Cholesky decomposition

FAQs GAMS
1

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?

2

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 \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.

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;