コレスキー分解

コレスキー分解　(Cholesky Decomposition)

正値対称行列Aを上三角行列Uとその転置行列Utの積
A=Ut*U
に分解することをコレスキー分解という．
与えられた分散共分散行列Aをコレスキー分解し得られるUを使うと，
独立の標準正規乱数から，Aの共分散をもった正規乱数に変換でき，
シミュレーションの乱数として利用される．

サンプルプログラム

****************************************;
* Cholesky Decomposition               *;
* 1999.2.5 翔                          *;
* input  A : Symmetric, all EigenVAL>0.*;
* output U : A=UtU                     *;
****************************************;


proc iml;reset log;

A={1  0.5 0.4 0.2,
  0.5 0.3 0.2 0.3,
  0.4 0.2 0.6 0.1,
  0.2 0.3 0.1 0.9};

eigenVAL=eigval(A);

n=ncol(A);
U=j(n,n,0);
do i=1 to n;
do j=1 to n;
  if i<=j then do;
    sUkiUkj=0;
    do k=1 to i-1;
      sUkiUkj=sUkiUkj+U[k,i]*U[k,j];
    end;
    if i=j then U[i,j]=sqrt(A[i,j]-sUkiUkj);
    else        U[i,j]=(A[i,j]-sUkiUkj)/U[i,i];
  end;
end;end;


UtU=t(U)*U;
print eigenVAL,A,UtU,U;

/*

 EIGENVAL

1.6354811
 0.793391
0.3684084
0.0027194
                   A

        1       0.5       0.4       0.2
      0.5       0.3       0.2       0.3
      0.4       0.2       0.6       0.1
      0.2       0.3       0.1       0.9
                  UTU

        1       0.5       0.4       0.2
      0.5       0.3       0.2       0.3
      0.4       0.2       0.6       0.1
      0.2       0.3       0.1       0.9
                   U

        1       0.5       0.4       0.2
        0 0.2236068         0 0.8944272
        0         0  0.663325 0.0301511
        0         0         0 0.2430862
*/


/*
もちろん関数一発でも求められる．
*/

half=half(a);
print half;

/*

                  HALF

        1       0.5       0.4       0.2
        0 0.2236068         0 0.8944272
        0         0  0.663325 0.0301511
        0         0         0 0.2430862
*/

サンプルプログラム その2

/*
行列演算は、一般にもちろんSAS/IMLソフトウェアで
行うことが推奨される。
しかし、コレスキー分解は、様々なプロシジャに
おいて内部で行われているため、それを応用すれば
SAS/IMLのライセンスがなくても実行できる。
以下のプログラムは、SAS/STATソフトウェアの
MIXEDプロシジャを利用したものである。
最終的に得られる行列は、上三角ではなく
下三角になっていることに注意。
*/

/*行列の準備。col1～や、rowという変数名は固定。*/
data a;
  input col1-col4;
  row+1;
  mean=0;
datalines;
1.0 0.5 0.4 0.2
0.5 0.3 0.2 0.3
0.4 0.2 0.6 0.1
0.2 0.3 0.1 0.9
;
run;

ods listing close;
ods output cholg=u(keep=col:);
proc mixed data=a;
  class row;
  parms /noiter;
  model mean=;
  random row*mean/type=un gdata=a gc;
run;
ods listing;

proc print data=u;
run;

/*下三角行列が返るが、本質的に同じ処理を行った結果である。*/
/*
OBS        Col1        Col2        Col3        Col4

 1       1.0000           0           0           0
 2       0.5000      0.2236           0           0
 3       0.4000           0      0.6633           0
 4       0.2000      0.8944     0.03015      0.2431
*/