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?

$ontext
Cholesky Decomposition
Contributed by Oyvind Hoveid (oyvind.hoveid@nilf.no)
$offtext

SETS
  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) 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;

VARIABLES
 BB(j,jj) Cholesky factors
 OBJ objective
EQUATIONS
 B1_EQ(j,jj) Cholesky factorisation
 OBJ_DF ;

 B1_EQ(j,jj) $ b_l(j,jj)..
  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;

$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, big = ww(r,r);
   if ( big ge 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, l;