以前，KL展開でKL展開(主成分分析)を作成しました． Jacobi法を使用した前回の方法では次元が少し増えるだけで処理回数が異常に増え， 実質的に使用不可となりますので今回QR法を用いて高速に演算が行えるKL展開を作成してみます． 但し，QR法とLU分解を用い，高速化のため行列演算を減らしたため，固有値が重複した際に， 固有ベクトルを列成分とする行列が正規直交行列にならくなり，固有値が重複する場合は 計算できないという制限がつきました．
今回高速化された代償として固有値が重複すると演算が出来なくなりました． それでは固有値が重複するということはどういう事なのか考えてみます．
固有値と固有ベクトルの意味はJacobi法で 書きました．KL展開ではベクトルの集合から共分散行列を考え，その固有値と固有ベクトルを 考えています．KL展開の導出のときに，分散は固有値の和であるということがわかっています． つまり，固有ベクトルの方向の分散の大きさが固有値で表されているという事になります．
例えば三次元空間において，上の図のように5つのベクトルがあるとします． このとき，x1，x2，x3の三次元空間のベクトルにも関わらず，赤の点線で示す通り，y1，y2の二次元平面 で完全に表現することができます．ここで固有値はy1，y2，y3方向の分散を表すため， y1の固有ベクトルが一番大きく，y2の固有ベクトルがその次，y3の固有ベクトルはゼロになります．
また，上の2つの図はKL展開をすると同じ結果になります． 上の図はx1方向に値を持っているものの， 分散という観点では上も下もy2の固有値(分散)はゼロになるからです． それではここで固有値が重複する場合を考えてみます．
固有値はそれに対応する固有ベクトル方向の分散ですので， KL展開後の基底において分散が等しいときに固有値が重解を持ちます． つまり次の図のようになったときだと考えられます．
すこし判り辛いですが，上の4つの点は上から見ると正方形を描き，x1-x2軸平面に対して並行 です．便宜上KL展開しても基底が変わらないものを選びましたが，変化しても同じです． ここでx1，x2，x3軸の変換後の基底をy1，y2，y3とすると，y1，y2方向の分散が 等しくなるのは自明です．これが固有値が重複するときの状況ですが， 例えばどれかの点が少しx1やx2方向にずれただけで固有値は重複しなくなります． このことから固有値が重複しないという条件は大してきつい条件ではなさそうです． 今回固有値が重複解を持つときを平面で考えましたが，3次元以上でも同じです．
以上より，固有値が重複する場合は計算できず，固有値がゼロのときには 重複する可能性は少し高いですが，そのときはそもそも固有ベクトルが必要ないので削ってしまいます． もっと言えば，仮に固有値が重複解を持っても，固有値が大きい部分であれば そもそもKL展開するようなデータではないですし，固有値が小さい部分であれば， まるごと削っても支障はないと思います．またデータが1つ増えたり減ったりするだけで， 分散のバランスも変わって来ると思いますので固有値は重解を持たないという前提で問題なさそうです．
QR法を利用すると高速に固有値を求める事が できますが，KL展開をするためには固有ベクトルを求めなくてはなりません． そこで基本に立ち帰って，
を満たすベクトルxを求めることで固有ベクトルを計算します．これは連立方程式を といているのと同じですので，LU分解を利用すれば 簡単に解くことができます．ただし，固有ベクトルの自由度は最低1以上になります． 自由度は固有値の重複度によって変わるわけですが，固有値が重複しなくても， 固有ベクトルは長さはいくらでもよく，方向のみが求められていますので， LU分解後は誤差を考えないとすると対角成分のどこかが必ずゼロになります． ちなみに固有値が重複すると自由度は重複した数だけ上り，対角成分のゼロの数も 増えて行きます．逆に考えると，対角成分のゼロが複数存在すれば 固有値が重複解を持ったと考えられます．
後退代入中に対角成分で割る事になりますが，結論から言うと対角成分がゼロであれば， それに対応するベクトルの要素の値はなんでも構いません．なので1とでもしておくことで， 残りの要素を決定することができます． 以上で固有ベクトルを求める事ができます．KL展開の実装については， 前回のJacobi法をQR法に差し替えただけなので特に問題は無いでしょう．
それではソースコードを見て行きます． 今回のソースコードはいくつも関数を使用していますが，その殆どは 以前作ったQR法，LU分解，Householder変換等の関数とその部品です． KLEncode2はKLEncodeがQR法に差し替えられただけです．引数が 増えていますが，これはdouble型がゼロかどうかの判断などに 許容誤差が必要になるので宣言しているだけです．そのため固有ベクトルの算出部分だけみてみます．
tmp = CreateMatrix(m->width, m->height);
for(i=0; i<m->height; i++){
for(j=0; j<m->width; j++){
if(i==j){
tmp->a[i*tmp->width + j] = m->a[i*m->width + j] -eigenvalue;
}
else{
tmp->a[i*tmp->width + j] = m->a[i*m->width + j];
}
}
}
order = LU(tmp);
min = GetMinvaluesIndex(tmp);
A-λIを計算し，それをLU分解しています．固有値が重複することは無いという前提で 計算していますので(もちろん後でチェックしています)，対角成分がゼロになるのは 一ヶ所です．そのため絶対値が一番小さいものがゼロであると考えて問題なさそうです．
dotproduct = 0;
for(i=m->height-1; i>=0; i--){
double sum=0;
for(j=m->height-1; j>i; j--){
sum += tmp->a[i*tmp->width + j]*eigvec->a[j*eigvec->width + index];
}
if(i==min){
dotproduct += SQR(eigvec->a[i*eigvec->width + index] = 1);
}
else{
if(fabs(tmp->a[i*tmp->width + i]) < eps){
fprintf(stderr, "Error: Include Multiple Root\n");
exit(1);
}
dotproduct += SQR(eigvec->a[i*eigvec->width + index] =
-sum/tmp->a[i*tmp->width + i]);
}
}
for(i=0; i<eigvec->height; i++)
eigvec->a[i*eigvec->width + index] /= sqrt(dotproduct);
この部分は後退代入を行っているのですが，iがminのとき，値はいくらでもいいので， 1を代入しています．しかし，それ以外のときに対角成分の絶対値がeps以下であると， 固有値に重複解が存在したということで終了します．最後に内積のルートで割っています． これは固有ベクトルを正規化しています．
以上で固有ベクトルを求めることができますが，前身消去は必要ないのかという疑問が湧くかと 思います．LU分解の際に下三角行列の対角成分は全て1にしていますので， 計算してみればわかりますが，解は全てゼロになります．そのためLUx=0はUx=0と同値です．
(*)記事中に何度か固有値が重複解を持つと計算出来ないと書きましたが， これは計算できないというよりも，同じ固有値を持つ固有ベクトル同士が直交する ベクトルを求められないという意味です．固有ベクトルの特徴として， 固有値が違えばそれぞれは必ず直交します．今回の場合，実対称行列ですので， 固有値が同じでも直交する固有ベクトルは存在し，実際Jacobi法では求まっています． QR法においても三重対角化のときからずっと直交行列を掛けることで変換してきていますので それをJacobi法のように単位行列に掛けながら保持しておけばその列ベクトルが 固有ベクトルなので固有値が重複するかどうかを考える必要はないのですが， 速度の面で行列演算を減らしたかったためにこのようになりました．
