# GAMS Support Wiki

### 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 to the GAMS-User List: 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;```
gams/cholesky_decomposition.1192857575.txt.gz · Last modified: 2007/10/20 07:19 by Franz Nelissen