特異値分解
Neuromagのソフトウェアには固有値分解の親戚の特異値分解がよく出てきます。
Scilabを使って雰囲気を味わいましょう。
stacksize(50000000);
loadmatfile(xgetfile('.../xxx.mat'));
getf('CalcSSP2.sce');
time=ave_0001(:,$);MEG=CalcSSP2(ave_0001(:,1:306));
scf();plot2d(time,MEG,rect=[-50,-600,500,600]);xset('use color',0);
体性感覚誘発磁場のoriginal波形です。
次にsvdを使って特異値分解します。
[Usvd,Ssvd,Dsvd]=svd(MEG');
Z=zeros(size(Ssvd,1),size(Ssvd,2)-size(Ssvd,1));
Ssvd=diag(Ssvd);
SS=Ssvd;SS(2:$)=0;SS=Usvd*[diag(SS),Z]*Dsvd';
scf();plot2d(time,SS',rect=[-50,-600,500,600]);xset('use color',0);
Scilabの関数svdはMatlabのsvdと同じものです。誤差はなさそうです。
行列Aを[U,S,D]=Aと分解します。
Sは降順に並んだ特異値の対角行列を含む行列です。
図は最大特異値に対応する成分だけで作図したものです。
SS=Ssvd;SS(4:$)=0;SS=Usvd*[diag(SS),Z]*Dsvd';
scf();plot2d(time,SS',rect=[-50,-600,500,600]);xset('use color',0);
図は特異値の上位1〜3に対応する成分だけで作図したものです。
SS=Ssvd;SS(8:$)=0;SS=Usvd*[diag(SS),Z]*Dsvd';
scf();plot2d(time,SS',rect=[-50,-600,500,600]);xset('use color',0);
図は特異値の上位1〜7に対応する成分だけで作図したものです。
originalの波形に近づいてきました。
sumSsvd=sum(abs(Ssvd).^2);
for n=1:length(Ssvd);if sum(abs(Ssvd(1:n)).^2)/sumSsvd>0.99;break;end;end;
n
n=40
SS=Ssvd;SS(41:$)=0;SS=Usvd*[diag(SS),Z]*Dsvd';
scf();plot2d(time,SS',rect=[-50,-600,500,600]);xset('use color',0);
図は特異値の上位1〜40に対応する成分だけで作図したものです。
固有値分解と同じようなことをしている気がします。
私は特異値の幾何学的な意味を理解していません。
特異値は固有値の√で、特異値分解は固有値分解よりも効率よく分解できる程度にしか理解していません。