User Tools

Site Tools


gams:cholesky_decomposition

This is an old revision of the document!


Cholesky Decomposition

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?

Contributed by Oyvind Hoveid: It is not that bad. The program enclosed seems to work. Basically, you need only one additional equation and one variable.

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;

A note contributed by Arne Drud: Oyvind Hoveid has a solution based on solving a model. You can include this inside a larger model so it is very flexible. 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 Oyvinds 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.

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;
IMPRESSUM / LEGAL NOTICEPRIVACY POLICY gams/cholesky_decomposition.1190861595.txt.gz · Last modified: 2007/10/20 07:19 (external edit)