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?
There are multiple ways to calculate the Cholesky factorization in GAMS: implement the algorithm in GAMS, use a model and solver, or use the cholesky
tool. These three methods are demonstrated in the model below. Thanks to Oyvind Hoveid and Arne Drud for contributing to this model.
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) >= ord(jj)) = yes; Parameters 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); b1(j,jj) = sum(n, (jn1(j,n) - b1(j,j))*(jn1(jj,n) - b1(jj,jj)))/card(n); * Symbols for the Cholesky decomposition algorithm Sets l(j,j) 'defines the lower triangle' r(j) 'the pivot row' / #j / npr(j) 'rows not yet used for pivot' / #j / npc(j) 'columns not yet used for pivot' / #j / parameters w(j,j) 'factor of b' ww(j,j) 'part of b not yet factored' big 'the diagonal pivot at each point in the factorization' piv 'the square root of the pivot' tol 'pivot tolerance' / 1e-10 /; * Cholesky decomposition algorithm in GAMS ww(j,jj) = b1(j,jj); loop(r, big = ww(r,r); if (big >= 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; * Using a CNS model to determine the Cholesky factors Variables BB(j,jj) 'Cholesky factors' Equations b1_eq(j,jj) 'Cholesky factorization'; b1_eq(b_l(j,jj)).. b1(j,jj) =e= sum(jjj, (BB(j,jjj)$b_l(j,jjj))*(BB(jj,jjj)$b_l(jj,jjj))); Model cholesky / b1_eq / ; BB.l(j,j) = 1; option limrow=0, limcol=0, solprint=silent; Solve cholesky using cns; display bb.l; * Using a cholesky tool to determine the Cholesky factors execute_unload 'a.gdx', j, b1; execute 'cholesky a.gdx j b1 b.gdx w'; execute_load 'b.gdx', w; display w;