Cholesky decomposition
Moderator: aileen
Forum rules
Please ask questions in the other sub-forums
Please ask questions in the other sub-forums
Cholesky decomposition
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?
Re: Cholesky decomposition
There are multiple ways to calculate the Cholesky factorization in GAMS: implement the algorithm in GAMS, use a model and solver, or use $libInclude linalg (see the header of the file in <sysdir>\inclib\linalg.gms for details). These three methods are demonstrated in the model below. Thanks to Oyvind Hoveid and Arne Drud for contributing to this model.
Code: Select all
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) >= ord(jj)) = yes;
Parameters
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);
b1(j,jj) = sum(n, (jn1(j,n) - b1(j,j))*(jn1(jj,n) - b1(jj,jj)))/card(n);
* Symbols for the Cholesky decomposition algorithm
Sets
l(j,j) 'defines the lower triangle'
r(j) 'the pivot row' / #j /
npr(j) 'rows not yet used for pivot' / #j /
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 /;
* Cholesky decomposition algorithm in GAMS
ww(j,jj) = b1(j,jj);
loop(r,
big = ww(r,r);
if (big >= tol,
piv = sqrt(big);
l(npr,r) = yes;
npr(r) = no;
npc(r) = no;
w(r ,r) = piv;
w(npr,r) = ww(npr,r) / piv;
ww(npr,npc) = ww(npr,npc) - ww(npr,r)*ww(npc,r)/big;
ww(npr,r) = 0;
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;
* Use $libInclude linalg to calculate the Cholesky factors
$libInclude linalg cholesky j b1 w
display w;