読者です 読者をやめる 読者になる 読者になる

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

前回の固有値分解に続き,今回は対称行列の一般化固有値問題を解きます.
一般化固有値問題には 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 の遅延評価について説明する予定です.
それでは.