gams:cholesky_decomposition

**This is an old revision of the document!**

*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: … 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 NOTICE
PRIVACY POLICY
gams/cholesky_decomposition.1192896914.txt.gz · Last modified: 2007/10/20 18:15 by Franz Nelissen