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

boost C++ lapack uBLAS

前回は lapack.dsyev による固有値分解（スペクトル分解）を紹介しました．
lapack にはもう一つ固有値分解を行う関数があります．
それが今回紹介する lapack.dsyevd になります．
dsyev との違いは分割統治法(divide-and-conquer method)により処理を高速化したことです．
早速実装してみましょう（前回の math.hpp と test.cpp は変更を加える必要がありません）．

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/syevd.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;
    namespace lapack = boost::numeric::bindings::lapack;

    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 = lapack::syevd('V', 'L', CQ, cd, lapack::optimal_workspace());
    BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic());

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

include するヘッダファイルを syev.hpp から syevd.hpp に変更したのと，syev の呼び出しを syevd に変更しただけです．

  • VC9
  • g++ (GCC) 3.4.4 (cygming special)

でのコンパイルを確認しています．
syev と比べどの程度高速化されるかは与える行列にかなり影響されるため，はっきりとした数字はいえませんが，私の経験ではだいたい 0% から 70% くらいの高速化になりました．


次回は一般化固有値問題について説明する予定です．
それでは．