特異値分解

C++ liboctave

特異値分解 (Singular value decomposition, SVD) の用途： 主成分分析，疑似逆行列の計算，機械学習，パタン認識，数値解析など．任意の形の行列を分解できる．

行列 (Matrix/ComplexMatrix) A に対する特異値分解は

  SVD  svd (const Matrix& A, SVD::type svd_type=SVD::std);
  SVD  svd (const Matrix& A, int& info, SVD::type svd_type=SVD::std);
  ComplexSVD  svd (const ComplexMatrix& A, SVD::type svd_type=SVD::std);
  ComplexSVD  svd (const ComplexMatrix& A, int& info, SVD::type svd_type=SVD::std);

によって得られる．

octave では，行列Aは

  A = USV*

のように分解される．ここで U, V はユニタリ行列， V* は V の共役転置を表す． S は対角成分が A の特異値（対角成分以外はゼロ）である． S は大きい順にソートされている （この情報は octave のマニュアルに基づくものではなく， Singular Value Decomposition -- LAPACK Users' Guide に記述されていることからの判断である）．

ここで SVD クラス, ComplexSVD クラス のオブジェクトには特異値分解の結果が代入される． info は計算後にゼロ以外なら何らかのエラーが発生したことを示すフラグである． SVD::type は列挙型で，

  00032 class
  00033 OCTAVE_API
  00034 SVD
  00035 {
  00036 public:
  00038   enum type
  00039     {
  00040       std,
  00041       economy,
  00042       sigma_only
  00043     };

のように定義されている．デフォルト値 SVD::std は通常の特異値分解（U,S,Vすべて計算）， SVD::sigma_only は S のみ計算， SVD::economy はエコノミーサイズ（UとVの不要な列や行を消去）で計算する．

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

  // S を返す :
  DiagMatrix   singular_values (void) const;
  // U を返す :
  Matrix   left_singular_matrix (void) const;
  // V を返す :
  Matrix   right_singular_matrix (void) const;

サンプル:

  Matrix A(3,6);
  stringstream ss (
    "  4   0   4   6   5   4"
    "  9   6   2   4   1   0"
    "  1   6   9   0   7   5"); ss >> A;
  int info;
  cout<<"A="<<endl<<A;
  // SVD オブジェクトの生成 :
  SVD  svd (A, info);
  cout<<"info= "<<info<<endl;
  cout<<"svd.singular_values()= "<<endl<< svd.singular_values();
  cout<<"svd.left_singular_matrix()= "<<endl<< svd.left_singular_matrix();
  cout<<"svd.right_singular_matrix()= "<<endl<< svd.right_singular_matrix();
  // USV* を計算（Aに一致するか?） :
  cout<<"svd.left_singular_matrix()*svd.singular_values()*svd.right_singular_matrix().transpose()= "<<endl
    << svd.left_singular_matrix()*svd.singular_values()*svd.right_singular_matrix().transpose();
  // エコノミーモードで特異値分解 :
  svd= SVD(A, info, SVD::economy);
  cout<<"info= "<<info<<endl;
  cout<<"svd.singular_values()= "<<endl<< svd.singular_values();
  cout<<"svd.left_singular_matrix()= "<<endl<< svd.left_singular_matrix();
  cout<<"svd.right_singular_matrix()= "<<endl<< svd.right_singular_matrix();
  // USV* を計算（Aに一致するか?） :
  cout<<"svd.left_singular_matrix()*svd.singular_values()*svd.right_singular_matrix().transpose()= "<<endl
    << svd.left_singular_matrix()*svd.singular_values()*svd.right_singular_matrix().transpose();

結果:

  A=
    4 0 4 6 5 4
    9 6 2 4 1 0
    1 6 9 0 7 5
  info= 0
  svd.singular_values()=
    17.6341 0 0 0 0 0
    0 9.51377 0 0 0 0
    0 0 6.12601 0 0 0
  svd.left_singular_matrix()=
    -0.509954 0.0830227 0.856186
    -0.498851 0.782325 -0.372982
    -0.700782 -0.617312 -0.357534
  svd.right_singular_matrix()=
    -0.410016 0.710098 -0.0472778 -0.556602 -0.109191 -0.060699
    -0.408176 0.104068 -0.715489 0.416853 0.235993 0.284919
    -0.529915 -0.384608 -0.0879892 0.0729675 -0.59143 -0.45652
    -0.286668 0.381283 0.595035 0.646778 0.00468758 -0.00584839
    -0.451064 -0.32834 0.229385 -0.224744 0.730871 -0.22679
    -0.314376 -0.289525 0.267234 -0.205602 -0.219999 0.80948
  svd.left_singular_matrix()*svd.singular_values()*svd.right_singular_matrix().transpose()=
    4 -8.88178e-16 4 6 5 4
    9 6 2 4 1 -1.77636e-15
    1 6 9 2.22045e-16 7 5
  info= 0
  svd.singular_values()=
    17.6341 0 0
    0 9.51377 0
    0 0 6.12601
  svd.left_singular_matrix()=
    -0.509954 0.0830227 0.856186
    -0.498851 0.782325 -0.372982
    -0.700782 -0.617312 -0.357534
  svd.right_singular_matrix()=
    -0.410016 0.710098 -0.0472778
    -0.408176 0.104068 -0.715489
    -0.529915 -0.384608 -0.0879892
    -0.286668 0.381283 0.595035
    -0.451064 -0.32834 0.229385
    -0.314376 -0.289525 0.267234
  svd.left_singular_matrix()*svd.singular_values()*svd.right_singular_matrix().transpose()=
    4 -8.88178e-16 4 6 5 4
    9 6 2 4 1 -1.77636e-15
    1 6 9 2.22045e-16 7 5