特異値分解

特異値分解(SVD,Singular Value Decomposition)

mxnの行列Aに対して,Aとその転置行列t(A)の積A*t(A)の固有値の絶対値の平方根を特異値という．
Aの特異値を対角成分にもつ対角行列Lと，2つの直行行列U,Vにより，行列Aを
A=U*L*V
の形に分解することを特異値分解という．
小さい特異値のいくつかを0にすることによって，近似的にAを次元の小さい行列の積に分解することができ，
大きい行列構造を縮退させたい場合に応用される．

Lは，固有値が一意であるので，特異値の定義から一意に定められる．
Uは，A*t(A)の固有ベクトルからなる行列として得られる．Vは，A,L,Uから一意に定まる．
（注，直行行列Uとは，U*t(U)=I(単位行列)となる行列である）

サンプルプログラム

******************************************************
svd.sas
1998.10.9 翔
特異値分解による次元の縮退
******************************************************;


options nocenter;

proc iml;

******************************************************
a=u*l*vを求める
******************************************************;


a={    -1        -1        12
  ,     0        11        -3
  ,     0        -1        -4
  ,     2        14        -3
  ,     0        10         8};


aat=a*t(a);
u=eigvec(aat);
l=diag(sqrt(abs(eigval(aat))));
v=ginv(l)*t(u)*a;

print u,l,v;
/*
                        U

-0.065942 0.7699142 0.0328341 0.2298009 0.5907598
0.5380201 -0.171431 0.7555206 -0.146094 0.2983133
-0.043958 -0.257571 0.1750334 0.9480556 -0.047739
0.6906779  -0.17509 -0.625984 0.1149005 0.2953799
0.4766737 0.5299474 0.0749515 0.1176287 -0.687374
                        L

 20.52127         0         0         0         0
        0 15.595187         0         0         0
        0         0 1.2913615         0         0
        0         0         0 1.7225E-7         0
        0         0         0         0 2.0648E-7
              V

0.0705267 0.9972262 -0.023789
-0.071823 0.0288632 0.9969997
-0.994921 0.0686065 -0.073659
1.2806E-8         0 1.9791E-9
4.3306E-8 8.0618E-9 1.2806E-9
*/

******************************************************
固有値を３つとも残せば，aをuq*qvで再現できる
******************************************************;

q=sqrt(l)*diag({1 1 1 0 0});
q2=l*diag({1 1 1 0 0});
uq=u*q;
qv=q*v;
uqqv=uq*qv;
print  uqqv,uq,qv;
/*
             UQQV

       -1        -1        12
-2.44E-15        11        -3
-1.44E-15        -1        -4
        2        14        -3
5.767E-15        10         8
                        UQ

-0.298722 3.0404483 0.0373121         0         0
2.4372529 -0.676995 0.8585591         0         0
-0.199131 -1.017165 0.1989046         0         0
 3.128799 -0.691443 -0.711356         0         0
2.1593512 2.0928018 0.0851735         0         0
              QV

 0.319489 4.5174754 -0.107765
-0.283635  0.113983 3.9372257
-1.130609 0.0779631 -0.083705
        0         0         0
        0         0         0
*/


******************************************************
固有値を１つ落として，aをuq*qvで近似的に分解する
******************************************************;

q=sqrt(l)*diag({1 1 0 0 0});
q2=l*diag({1 1 0 0 0});
uq=u*q;
qv=q*v;
uqqv=uq*qv;
print  uqqv,uq,qv;

/*
             UQQV

-0.957815 -1.002909 12.003123
0.9706947 10.933064 -2.928134
0.2248833 -1.015507 -3.983351
 1.195734  14.05546 -3.059544
0.0962979 9.9933596 8.0071295
                        UQ

-0.298722 3.0404483         0         0         0
2.4372529 -0.676995         0         0         0
-0.199131 -1.017165         0         0         0
 3.128799 -0.691443         0         0         0
2.1593512 2.0928018         0         0         0
              QV

 0.319489 4.5174754 -0.107765
-0.283635  0.113983 3.9372257
        0         0         0
        0         0         0
        0         0         0
*/


******************************************************
[#p1305244]
コール関数svd()でも，簡単に求められる
******************************************************;

call svd(u,l,v,a);

print u,l,v;


/*

              U

-0.065942 0.7699142 -0.032834
0.5380201 -0.171431 -0.755521
-0.043958 -0.257571 -0.175033
0.6906779  -0.17509 0.6259842
0.4766737 0.5299474 -0.074952
    L

 20.52127
15.595187
1.2913615
              V

0.0705267 -0.071823 0.9949208
0.9972262 0.0288632 -0.068607
-0.023789 0.9969997 0.0736595

*/