# GAMS Support Wiki

### Site Tools

gams:cholesky_decomposition

# 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?

```\$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;```
gams/cholesky_decomposition.txt · Last modified: 2007/10/21 04:22 by franz