Boost.uBLAS で対称行列の固有値分解(Eigenvalue Decomposition)

boost C++ lapack uBLAS

前回のQR分解に続き，今回は対称行列の固有値分解（スペクトル分解）を実装します．
QR分解のときと違い，lapack.dsyev を用いるだけなので，特異値分解のときと同様に非常に簡単に実装することができます．
以下に私の実装例を示します．

math.hpp

#ifndef MATH_HPP_20080924
#define MATH_HPP_20080924

#if defined(_MSC_VER) && (_MSC_VER >= 1020)
# pragma once
#endif


#include <vector>
#include <boost/numeric/ublas/fwd.hpp>


namespace math {
    void eigen(const boost::numeric::ublas::symmetric_matrix<double>& A, boost::numeric::ublas::matrix<double>& Q, std::vector<double>& d);
}


#endif

math.cpp

#include "math.hpp"
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#define BOOST_NUMERIC_BINDINGS_USE_CLAPACK
#include <boost/numeric/bindings/lapack/syev.hpp>
#include <boost/numeric/bindings/traits/std_vector.hpp>
#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
#undef  BOOST_NUMERIC_BINDINGS_USE_CLAPACK


void math::eigen(const boost::numeric::ublas::symmetric_matrix<double>& A, boost::numeric::ublas::matrix<double>& Q, std::vector<double>& d)
{
    namespace ublas = boost::numeric::ublas;

    std::vector<double>                        cd(A.size1());
    ublas::matrix<double, ublas::column_major> CQ(A.size1(), A.size2());
    int                                        info;

    for (std::size_t i = 0; i < A.size1(); ++i) {
        for (std::size_t j = 0; j <= i; ++j) {
            CQ(i, j) = A(i, j);
        }
    }

    info = boost::numeric::bindings::lapack::syev('V', 'L', CQ, cd);
    BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic());

    Q = CQ;
    d.swap(cd);
}

math.cpp に固有値分解を計算する関数 eigen を実装しました．math.cpp をプロジェクトに加え（他のソースコードと一緒にコンパイル），math.hpp を include することにより eigen が使用可能になります．
関数 eigen は N x N 対称行列 A を受け取り，これを N x N 正規直交行列 Q と N x N 対角行列 D と Q の転置の積に因数分解します．N x N 正規直交行列が関数の引数 Q に，N x N 対角行列 D の対角成分が関数の引数 d に保存されます．
N x N 対角行列 D の対角成分は行列 A の固有値になっており，固有値 D(i, i) は 固有ベクトル column(Q, i) に対応します．
N x N 対角行列 D の対角成分は昇順になっており，順番に関数の引数 d に保存されます．
コードの簡単な説明ですが，まず入力行列 A の下三角行列を CQ に保存します．
bindings.lapack.syev を呼び出すと CQ は A = Q x D x trans(Q) の Q に，cd は D の対角成分になるため，これを結果として関数の引数に保存して終了です．
以下のテストコードで eigen の動作を確認してみます．

test.cpp

#include "math.hpp"
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/symmetric.hpp>


#ifdef _MSC_VER
# pragma comment(lib, "libf2c.lib")
# pragma comment(lib, "BLAS.lib")
# pragma comment(lib, "clapack.lib")
#endif


int main()
{
    namespace ublas = boost::numeric::ublas;

    ublas::symmetric_matrix<double> A(2, 2);
    ublas::matrix<double>           Q;
    std::vector<double>             d;

    A(0, 0) = 2;
    A(1, 0) = 1; A(1, 1) = 2;

    math::eigen(A, Q, d);

    // 以下結果の表示
    // 計算結果を http://en.wikipedia.org/wiki/Eigenvalue,_eigenvector_and_eigenspace#Example と比較してみましょう
    for (std::size_t i = 0; i < d.size(); ++i) {
        std::cout << d[i] << ' ' << ublas::column(Q, i) << std::endl;
    }

    return 0;
}
  • VC9
  • g++ (GCC) 3.4.4 (cygming special)

でのコンパイルを確認しています．
実行してみると，正しく固有値分解できていることが分かります．
エラーが出てコンパイルできない場合は，前回の記事を見て正しくインストールできているか確認してみてください．


次回は分割統治法を用いたより高速な固有値分解を実装する予定です．
とはいっても，syev を syevd に機械的に置換するだけでほぼ終了ですが．
それでは．