Boost.uBLAS で対称行列の一般化固有値分解(Generalized eigenvalue decomposition)

boost C++ lapack uBLAS

前回の固有値分解に続き，今回は対称行列の一般化固有値問題を解きます．
一般化固有値問題には lapack.dsygv を用いるだけで，コードも固有値分解のときとほとんど同じで，簡単に実装することができます．
以下に私の実装例を示します．

math.hpp

#ifndef MATH_HPP_20080926
#define MATH_HPP_20080926

#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, const boost::numeric::ublas::symmetric_matrix<double>& B, boost::numeric::ublas::matrix<double>& Z, 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/sygv.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, const boost::numeric::ublas::symmetric_matrix<double>& B, boost::numeric::ublas::matrix<double>& Z, std::vector<double>& d)
{
    namespace ublas  = boost::numeric::ublas;
    namespace lapack = boost::numeric::bindings::lapack;

    BOOST_UBLAS_CHECK(A.size1() == B.size1(), ublas::external_logic());

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

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

    info = lapack::sygv(1, 'V', 'L', CA, CB, cd, lapack::optimal_workspace());
    BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic());

    Z = CA;
    d.swap(cd);
}

math.cpp に一般化固有値問題を解く関数 eigen を実装しました．math.cpp をプロジェクトに加え（他のソースコードと一緒にコンパイル），math.hpp を include することにより eigen が使用可能になります．
関数 eigen は N x N 対称行列 A と N x N 正値対称行列 B を受け取り，Ax = (lambda)Bx となるすべて（N 個）の N 次元固有ベクトル x と固有値 lambda を求めます．N 個の固有値 lambda は昇順に関数の引数 d に，その固有値 lambda に対応する N 次元固有ベクトル x は関数の引数 Z の各列（column(Z, 0) から column(Z, N - 1)）に順番に保存されます．Z が N x N 行列になることに注意してください．
また，各固有ベクトルは trans(Z) x B x Z が単位行列になるように正規化されています．
コードの簡単な説明ですが，まず入力行列 A と B の下三角行列を CA と CB に保存します．
bindings.lapack.sygv を呼び出すと CA の各列は固有ベクトルに，cd は対応する固有値になるため，これを結果として関数の引数に保存して終了です．
また，行列 CB の上（下）三角成分は Cholesky 分解 B = trans(U) x U あるいは B = L x trans(L) の U あるいは L になるのですが，これはどこにも保存せずに破棄しています．
以下のテストコードで 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(3, 3), B(3, 3);
    ublas::matrix<double>           Z;
    std::vector<double>             d;

    A(0, 0) = 1;
    A(1, 0) = 1;   A(1, 1) = 1;
    A(2, 0) = 0.5; A(2, 1) = 0.25; A(2, 2) = 2;
    
    B(0, 0) = 2;
    B(1, 0) = 2; B(1, 1) = 5;
    B(2, 0) = 2; B(2, 1) = 5; B(2, 2) = 11;

    math::eigen(A, B, Z, d);

    // 以下結果の確認
    for (std::size_t i = 0; i < d.size(); ++i) {
        std::cout << d[i] << ' ' << ublas::column(Z, i) << ' ' << std::endl;

        // Ax = (lambda)Bx ?
        BOOST_UBLAS_CHECK(
            ublas::detail::expression_type_check(
                ublas::prod(A, ublas::column(Z, i)),       // Ax
                d[i] * ublas::prod(B, ublas::column(Z, i)) // (lambda)Bx
            ),
            ublas::internal_logic()
        );
    }

    std::cout << ublas::prod(ublas::trans(Z), ublas::matrix<double>(ublas::prod(B, Z))) << std::endl; // 単位行列になるはず

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

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


次回は uBLAS の遅延評価について説明する予定です．
それでは．