Differences

This shows you the differences between two versions of the page.

--- gams:cholesky_decomposition [2007/10/21 06:22]
Franz Nelissen
+++ gams:cholesky_decomposition [2020/05/27 11:30] (current)
Michael Bussieck
@@ Line 3: / Line 3: @@
 //Q: I have an estimate of a variance-covariance matrix and I need the Cholesky decomposition of it. Does anybody have any suggestions on how to program this in GAMS?//
+There are multiple ways to calculate the Cholesky factorization in GAMS: implement the algorithm in GAMS, use a model and solver, or use the ''cholesky'' tool. These three methods are demonstrated in the model below. Thanks to Oyvind Hoveid and Arne Drud for contributing to this model.
 <code>
-$ontext
+Sets
-Cholesky Decomposition
+     j          'variables'                /j1*j10/
-Contributed by Oyvind Hoveid (oyvind.hoveid@nilf.no)
+     n          'observations'             /n1*n15/
-$offtext
+     b_l(j,j)   'lower triangular marix';
-SETS
+Alias (j,jj,jjj);
-  j variables    /j1*j10/,
-  n observations /n1*n10/,
-  b_l(j,j) lower triangular marix;
-ALIAS (j,jj,jjj);
+b_l(j,jj) $ (ord(j) >= ord(jj)) = yes;
-b_l(j,jj) $ (ORD(j) GE ORD(jj)) = YES;
+Parameters
+ jn1(j,n)   'variable observations'
+ b1(j,jj)   'covariance matrix';
-PARAMETER
-  jn1(j,n) variable observations
-  b1(j,jj) covariance matrix;
 * forming a covariance matrix
-jn1(j,n) = NORMAL(0,1);
+jn1(j,n) = normal(0,1);
-b1(j,j) = SUM(n, jn1(j,n)) /CARD(n);
+b1(j,j)  = sum(n, jn1(j,n))/card(n);
-DISPLAY b1;
+b1(j,jj) = sum(n, (jn1(j,n) - b1(j,j))*(jn1(jj,n) - b1(jj,jj)))/card(n);
-b1(j,jj) = SUM(n, (jn1(j,n) - b1(j,j))*(jn1(jj,n) -b1(jj,jj)))/CARD(n);
-DISPLAY b1;
-VARIABLES
+* Symbols for the Cholesky decomposition algorithm
- BB(j,jj) Cholesky factors
+Sets
- OBJ objective
+    l(j,j)  'defines the lower triangle'
-EQUATIONS
+    r(j)    'the pivot row'                  / #j /
- B1_EQ(j,jj) Cholesky factorisation
+    npr(j)  'rows not yet used for pivot'    / #j /
- OBJ_DF ;
+    npc(j)  'columns not yet used for pivot' / #j /
+parameters
+    w(j,j)  'factor of b'
+    ww(j,j) 'part of b not yet factored'
+    big     'the diagonal pivot at each point in the factorization'
+    piv     'the square root of the pivot'
+    tol     'pivot tolerance' / 1e-10 /;
- B1_EQ(j,jj) $ b_l(j,jj)..
+* Cholesky decomposition algorithm in GAMS
-  b1(j,jj) =E= SUM(jjj, (BB(j,jjj) $ b_l(j,jjj))*(BB(jj,jjj) $ b_l(jj,jjj)));
-OBJ_DF.. OBJ =E= 1;
-MODEL CHOLESKY /ALL/ ;
-BB.L(j,j) = 1;
-SOLVE CHOLESKY MAXIMISING OBJ USING NLP
-display bb.l;
-</code>
-<code>
-$ontext
-Contributed by Arne Drud (adrud@arki.dk):
-If you just need the Cholesky decomposition then you can compute it in GAMS with a
-loop as shown below. It is an extension of Oyvind's program in which I increased the number of observations
-to ensure that the matrix was positive definite (10 observations gives a matrix that is very close to
-being singular). There are additional example discussing Cholesky decompositions and positive semidefinite
-matrices in the GAMS-L archive -- it is a popular area.
-$offtext
-SETS
- j           variables                            /j1*j10/
- n          observations                          /n1*n15/
- b_l(j,j)   lower triangular marix;
-ALIAS (j,jj,jjj);
-b_l(j,jj) $ (ORD(j) GE ORD(jj)) = YES;
-PARAMETER
- jn1(j,n)  variable observations
- b1(j,jj)   covariance matrix;
-* forming a covariance matrix
-jn1(j,n) = NORMAL(0,1);
-b1(j,j) = SUM(n, jn1(j,n)) /CARD(n);
-DISPLAY b1;
-b1(j,jj) = SUM(n, (jn1(j,n) - b1(j,j))*(jn1(jj,n) - b1(jj,jj)))/CARD(n);
-DISPLAY b1;
-set l(j,j) defines the lower triangle;
-parameters  w(j,j)    factor of b
-parameters  ww(j,j)   part of b not yet factored
-sets   r(j)   the pivot row
-       npr(j) rows not yet used for pivot
-       npc(j) columns not yet used for pivot
-scalar big    the diagonal pivot at each point in the factorization
-       piv    the square root of the pivot
-       tol    pivot tolerance;
-r(j)     = yes;
-npr(j)   = yes;
-npc(j)   = yes;
 ww(j,jj) = b1(j,jj);
-tol      = 1e-10;
+loop(r,
-loop(r, big = ww(r,r);
+   big = ww(r,r);
-   if ( big ge tol,
+   if (big >= tol,
       piv = sqrt(big);
       l(npr,r) = yes;
@@ Line 102: / Line 51: @@
       ww(r,npc) = 0;
       ww(r,r)   = 0;
-      );
    );
+);
+display w;
+* Using a CNS model to determine the Cholesky factors
+Variables
+   BB(j,jj) 'Cholesky factors'
+Equations
+   b1_eq(j,jj) 'Cholesky factorization';
+b1_eq(b_l(j,jj))..
+  b1(j,jj) =e= sum(jjj, (BB(j,jjj)$b_l(j,jjj))*(BB(jj,jjj)$b_l(jj,jjj)));
+Model cholesky / b1_eq / ;
+BB.l(j,j) = 1;
+option limrow=0, limcol=0, solprint=silent;
+Solve cholesky using cns;
+display bb.l;
-display w, l;
+* Using a cholesky tool to determine the Cholesky factors
+execute_unload 'a.gdx', j, b1;
+execute 'cholesky a.gdx j b1 b.gdx w';
+execute_load 'b.gdx', w;
+display w;
 </code>

GAMS Support Wiki

User Tools

Site Tools

Differences

Page Tools