固有値分解 - akihiko’s tech note

固有値分解 (Eigendecomposition) の用途：主成分分析など多数．正方行列に対してのみ適用可能（非正方の場合は特異値分解）．

行列 (Matrix) A に対する固有値分解は，

  EIG  eig (const Matrix& A, bool calc_eigenvectors=true);
  EIG  eig (const ComplexMatrix& A, bool calc_eigenvectors=true);

によって得られる．ここで EIG クラスのオブジェクト eig には固有値分解の結果，すなわち固有値 (eigenvalues) と固有ベクトル (eigenvectors) が格納される． calc_eigenvectors は固有ベクトルを計算するかどうかを表すフラグである（省略してもよい）．

EIG オブジェクトは以下の public メンバ関数を持つ:

  // 固有値を返す(ComplexColumnVector) :
  ComplexColumnVector  eigenvalues (void) const;
  // 固有ベクトルを返す(ComplexMatrix) :
  ComplexMatrix  eigenvectors (void) const;

個々の固有値には eig.eigenvalues()(i), i:int のようにしてアクセスできる（iの範囲はAの行数と同じ）． eig.eigenvalues() は固有値を縦に並べた列ベクトル (ComplexColumnVector) である．個々の固有ベクトルには eig.eigenvectors().column(i), i:int のようにしてアクセスできる． eig.eigenvectors() は固有ベクトル（列ベクトル）を横に並べた行列 (ComplexMatrix) である．

固有値・固有ベクトルのいずれも，要素の一部または全部が複素数になる場合がある．実部を取り出す必要がある場合は， real 関数を使う． real(eig.eigenvalues()) とすると実数の ColumnVector ， real(eig.eigenvalues()(i)) とすると実数 (double)， real(eig.eigenvectors()) とすると実数の Matrix ， real(eig.eigenvectors().column(i)) とすると実数の ColumnVector がそれぞれ得られる．

octave のマニュアルにあるように， The eigenvalues returned by eig are not ordered. ，つまり eig.eigenvalues() はソートされていない．サンプルデータの共分散行列に EIG を適用すると昇順にソートされているように見えるが，保証はない．

主成分分析などで固有値の大きいものから順に取得したい場合， eig.eigenvalues() をソートすればよいのだが，固有ベクトルのソートも同時に行わなければならない．この対処のひとつとして， eig.eigenvalues() を直接ソートせずに， {std::vector idx, idx[n] はn番目に大きい固有値のインデックス} となるようなインデックスベクトル idx を定義する方法が考えられる．具体的な実装は主成分分析のサンプルプログラムを参照．ただし，主成分分析をする場合は特異値分解を用いることが多いようで（精度がいいらしい），この場合， liboctave では固有値と固有ベクトルはソートされたものが得られるから，特に気にする必要はない．

サンプル:

  Matrix A(4,4);
  stringstream ss (
    "  4.0  2.0  3.0  1.0"
    "  0.0 -1.0  1.0  2.0"
    "  1.0  3.0  2.0  4.0"
    "  0.0 -1.0 -5.0 -4.0"); ss >> A;
  int info;
  cout<<"A="<<endl<<A;
  // EIG オブジェクトの生成 :
  EIG  eig (A, info);
  cout<<"info= "<<info<<endl;
  cout<<"eig.eigenvalues()= "<<endl<< eig.eigenvalues();
  cout<<"eig.eigenvectors()= "<<endl<< eig.eigenvectors();
  // A を対称にする :
  A= 0.5*(A+A.transpose());
  cout<<"A="<<endl<<A;
  // 再度固有値分解 :
  eig= EIG(A,info);
  cout<<"info= "<<info<<endl;
  cout<<"eig.eigenvalues()= "<<endl<< eig.eigenvalues();
  cout<<"eig.eigenvectors()= "<<endl<< eig.eigenvectors();
  // 固有値，固有ベクトルを実数値で抽出する :
  cout<<"real(eig.eigenvalues()(2))= "<< real(eig.eigenvalues()(2))<<endl;
  cout<<"real(eig.eigenvectors().column(2))= "<<endl<< real(eig.eigenvectors().column(2));

結果::

  A=
    4 2 3 1
    0 -1 1 2
    1 3 2 4
    0 -1 -5 -4
  info= 0
  eig.eigenvalues()=
    (4.4782,0)
    (-0.524237,3.27929)
    (-0.524237,-3.27929)
    (-2.42973,0)
  eig.eigenvectors()=
    (0.9733,0) (0.0591816,-0.335765) (0.0591816,0.335765) (-0.168139,0)
    (-0.00621851,0) (0.0758315,0.282186) (0.0758315,-0.282186) (0.796196,0)
    (0.197962,0) (0.450522,0.382927) (0.450522,-0.382927) (0.0231463,0)
    (-0.116014,0) (-0.669908,0) (-0.669908,-0) (-0.580745,0)
  A=
    4 1 2 0.5
    1 -1 2 0.5
    2 2 2 -0.5
    0.5 0.5 -0.5 -4
  info= 0
  eig.eigenvalues()=
    (-4.24974,0)
    (-1.80298,0)
    (1.25496,0)
    (5.79776,0)
  eig.eigenvectors()=
    (0.0718392,0) (-0.0357097,0) (0.618216,0) (-0.781904,0)
    (0.2315,0) (0.870549,0) (-0.331435,0) (-0.280539,0)
    (-0.17344,0) (-0.402036,0) (-0.70642,0) (-0.556108,0)
    (-0.95455,0) (0.28149,0) (0.0945012,0) (-0.0258393,0)
  real(eig.eigenvalues()(2))= 1.25496
  real(eig.eigenvectors().column(2))=
    0.618216
    -0.331435
    -0.70642
    0.0945012

eig.eigenvalues() などの結果が (r,i) のように書かれているのは，複素数だからだ． r が実部で， i が虚部である．