KL展開2 をテンプレートにして作成

開始行:
RIGHT:寄稿：東條遼平 * KL展開の高速化 [#fdb7742e] 以前、KL展開でKL展開(主成分分析)を作成しました。 Jacobi法を使用した前回の方法では次元が少し増えるだけで処... 実質的に使用不可となりますので今回QR法を用いて高速に演算... 但し、QR法とLU分解を用い、高速化のため行列演算を減らした... 固有ベクトルを列成分とする行列が正規直交行列にならくなり... 計算できないという制限がつきました。 * 固有値の意味 [#rc480819] 今回高速化された代償として固有値が重複すると演算が出来な... それでは固有値が重複するということはどういう事なのか考え... 固有値と固有ベクトルの意味はJacobi法で 書きました。KL展開ではベクトルの集合から共分散行列を考え... 考えています。KL展開の導出のときに、分散は固有値の和であ... つまり、固有ベクトルの方向の分散の大きさが固有値で表され... #ref(klt2_1.png,nolink) 例えば三次元空間において、上の図のように5つのベクトルがあ... このとき、x1、x2、x3の三次元空間のベクトルにも関わらず、... で完全に表現することができます。ここで固有値はy1、y2、y3... y1の固有ベクトルが一番大きく、y2の固有ベクトルがその次、y... #ref(klt2_2.png,nolink) #ref(klt2_3.png,nolink) また、上の2つの図はKL展開をすると同じ結果になります。 上... 分散という観点では上も下もy2の固有値(分散)はゼロになるか... それではここで固有値が重複する場合を考えてみます。 固有値はそれに対応する固有ベクトル方向の分散ですので、 KL展開後の基底において分散が等しいときに固有値が重解を持... つまり次の図のようになったときだと考えられます。 #ref(klt2_4.png,nolink) すこし判り辛いですが、上の4つの点は上から見ると正方形を描... です。便宜上KL展開しても基底が変わらないものを選びました... ここでx1、x2、x3軸の変換後の基底をy1、y2、y3とすると、y1... 等しくなるのは自明です。これが固有値が重複するときの状況... 例えばどれかの点が少しx1やx2方向にずれただけで固有値は重... このことから固有値が重複しないという条件は大してきつい条... 今回固有値が重複解を持つときを平面で考えましたが、3次元以... 以上より、固有値が重複する場合は計算できず、固有値がゼロ... 重複する可能性は少し高いですが、そのときはそもそも固有ベ... もっと言えば、仮に固有値が重複解を持っても、固有値が大き... そもそもKL展開するようなデータではないですし、固有値が小... まるごと削っても支障はないと思います。またデータが1つ増え... 分散のバランスも変わって来ると思いますので固有値は重解を... * 固有ベクトルの算出 [#d33e00c1] QR法を利用すると高速に固有値を求める事が できますが、KL展開をするためには固有ベクトルを求めなくて... そこで基本に立ち帰って、 #ref(img190.png,nolink) を満たすベクトルxを求めることで固有ベクトルを計算します。... といているのと同じですので、LU分解を利用すれば 簡単に解くことができます。ただし、固有ベクトルの自由度は... 自由度は固有値の重複度によって変わるわけですが、固有値が... 固有ベクトルは長さはいくらでもよく、方向のみが求められて... LU分解後は誤差を考えないとすると対角成分のどこかが必ずゼ... ちなみに固有値が重複すると自由度は重複した数だけ上り、対... 増えて行きます。逆に考えると、対角成分のゼロが複数存在す... 固有値が重複解を持ったと考えられます。 後退代入中に対角成分で割る事になりますが、結論から言うと... それに対応するベクトルの要素の値はなんでも構いません。な... 残りの要素を決定することができます。 以上で固有ベクトルを求める事ができます。KL展開の実装につ... 前回のJacobi法をQR法に差し替えただけなので特に問題は無い... * ソースコード [#v9f9a23a] それではソースコードを見て行きます。 今回のソースコードはいくつも関数を使用していますが、その... 以前作った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->w... } else{ tmp->a[i*tmp->width + j] = m->a[i*m->w... } } } 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... } if(i==min){ dotproduct += SQR(eigvec->a[i*eigvec->width + ... } 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 + ... -sum/tmp->a[i*tmp->width + ... } } for(i=0; i<eigvec->height; i++) eigvec->a[i*eigvec->width + index] /= sqrt(dotpr... この部分は後退代入を行っているのですが、iがminのとき、値... 1を代入しています。しかし、それ以外のときに対角成分の絶対... 固有値に重複解が存在したということで終了します。最後に内... これは固有ベクトルを正規化しています。 以上で固有ベクトルを求めることができますが、前身消去は必... 思います。LU分解の際に下三角行列の対角成分は全て1にしてい... 計算してみればわかりますが、解は全てゼロになります。その... (*)記事中に何度か固有値が重複解を持つと計算出来ないと書き... これは計算できないというよりも、同じ固有値を持つ固有ベク... ベクトルを求められないという意味です。固有ベクトルの特徴... 固有値が違えばそれぞれは必ず直交します。今回の場合、実対... 固有値が同じでも直交する固有ベクトルは存在し、実際Jacobi... QR法においても三重対角化のときからずっと直交行列を掛ける... それをJacobi法のように単位行列に掛けながら保持しておけば... 固有ベクトルなので固有値が重複するかどうかを考える必要は... 速度の面で行列演算を減らしたかったためにこのようになりま... - &ref(main.c); - &ref(calculation.c); - &ref(calculation.h);

終了行:
RIGHT:寄稿：東條遼平 * KL展開の高速化 [#fdb7742e] 以前、KL展開でKL展開(主成分分析)を作成しました。 Jacobi法を使用した前回の方法では次元が少し増えるだけで処... 実質的に使用不可となりますので今回QR法を用いて高速に演算... 但し、QR法とLU分解を用い、高速化のため行列演算を減らした... 固有ベクトルを列成分とする行列が正規直交行列にならくなり... 計算できないという制限がつきました。 * 固有値の意味 [#rc480819] 今回高速化された代償として固有値が重複すると演算が出来な... それでは固有値が重複するということはどういう事なのか考え... 固有値と固有ベクトルの意味はJacobi法で 書きました。KL展開ではベクトルの集合から共分散行列を考え... 考えています。KL展開の導出のときに、分散は固有値の和であ... つまり、固有ベクトルの方向の分散の大きさが固有値で表され... #ref(klt2_1.png,nolink) 例えば三次元空間において、上の図のように5つのベクトルがあ... このとき、x1、x2、x3の三次元空間のベクトルにも関わらず、... で完全に表現することができます。ここで固有値はy1、y2、y3... y1の固有ベクトルが一番大きく、y2の固有ベクトルがその次、y... #ref(klt2_2.png,nolink) #ref(klt2_3.png,nolink) また、上の2つの図はKL展開をすると同じ結果になります。 上... 分散という観点では上も下もy2の固有値(分散)はゼロになるか... それではここで固有値が重複する場合を考えてみます。 固有値はそれに対応する固有ベクトル方向の分散ですので、 KL展開後の基底において分散が等しいときに固有値が重解を持... つまり次の図のようになったときだと考えられます。 #ref(klt2_4.png,nolink) すこし判り辛いですが、上の4つの点は上から見ると正方形を描... です。便宜上KL展開しても基底が変わらないものを選びました... ここでx1、x2、x3軸の変換後の基底をy1、y2、y3とすると、y1... 等しくなるのは自明です。これが固有値が重複するときの状況... 例えばどれかの点が少しx1やx2方向にずれただけで固有値は重... このことから固有値が重複しないという条件は大してきつい条... 今回固有値が重複解を持つときを平面で考えましたが、3次元以... 以上より、固有値が重複する場合は計算できず、固有値がゼロ... 重複する可能性は少し高いですが、そのときはそもそも固有ベ... もっと言えば、仮に固有値が重複解を持っても、固有値が大き... そもそもKL展開するようなデータではないですし、固有値が小... まるごと削っても支障はないと思います。またデータが1つ増え... 分散のバランスも変わって来ると思いますので固有値は重解を... * 固有ベクトルの算出 [#d33e00c1] QR法を利用すると高速に固有値を求める事が できますが、KL展開をするためには固有ベクトルを求めなくて... そこで基本に立ち帰って、 #ref(img190.png,nolink) を満たすベクトルxを求めることで固有ベクトルを計算します。... といているのと同じですので、LU分解を利用すれば 簡単に解くことができます。ただし、固有ベクトルの自由度は... 自由度は固有値の重複度によって変わるわけですが、固有値が... 固有ベクトルは長さはいくらでもよく、方向のみが求められて... LU分解後は誤差を考えないとすると対角成分のどこかが必ずゼ... ちなみに固有値が重複すると自由度は重複した数だけ上り、対... 増えて行きます。逆に考えると、対角成分のゼロが複数存在す... 固有値が重複解を持ったと考えられます。 後退代入中に対角成分で割る事になりますが、結論から言うと... それに対応するベクトルの要素の値はなんでも構いません。な... 残りの要素を決定することができます。 以上で固有ベクトルを求める事ができます。KL展開の実装につ... 前回のJacobi法をQR法に差し替えただけなので特に問題は無い... * ソースコード [#v9f9a23a] それではソースコードを見て行きます。 今回のソースコードはいくつも関数を使用していますが、その... 以前作った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->w... } else{ tmp->a[i*tmp->width + j] = m->a[i*m->w... } } } 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... } if(i==min){ dotproduct += SQR(eigvec->a[i*eigvec->width + ... } 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 + ... -sum/tmp->a[i*tmp->width + ... } } for(i=0; i<eigvec->height; i++) eigvec->a[i*eigvec->width + index] /= sqrt(dotpr... この部分は後退代入を行っているのですが、iがminのとき、値... 1を代入しています。しかし、それ以外のときに対角成分の絶対... 固有値に重複解が存在したということで終了します。最後に内... これは固有ベクトルを正規化しています。 以上で固有ベクトルを求めることができますが、前身消去は必... 思います。LU分解の際に下三角行列の対角成分は全て1にしてい... 計算してみればわかりますが、解は全てゼロになります。その... (*)記事中に何度か固有値が重複解を持つと計算出来ないと書き... これは計算できないというよりも、同じ固有値を持つ固有ベク... ベクトルを求められないという意味です。固有ベクトルの特徴... 固有値が違えばそれぞれは必ず直交します。今回の場合、実対... 固有値が同じでも直交する固有ベクトルは存在し、実際Jacobi... QR法においても三重対角化のときからずっと直交行列を掛ける... それをJacobi法のように単位行列に掛けながら保持しておけば... 固有ベクトルなので固有値が重複するかどうかを考える必要は... 速度の面で行列演算を減らしたかったためにこのようになりま... - &ref(main.c); - &ref(calculation.c); - &ref(calculation.h);

ページ名: