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

Boost.uBLAS でコレスキー分解

今回はコレスキー分解を実装します.コレスキー分解はLU分解の特別な場合であり,入力行列 A が正値対称行列なら,これを A = L x trans(L) あるいは A = trans(U) x U に分解することができます.コレスキー分解は lapack.dpotrf で計算することができます.
以下に私の実装例を示します.

math.hpp

#ifndef MATH_HPP_20081016
#define MATH_HPP_20081016

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


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


namespace math {
    void cholesky_factorize(const boost::numeric::ublas::symmetric_matrix<double>& A, boost::numeric::ublas::triangular_matrix<double>& L);
    void cholesky_factorize(const boost::numeric::ublas::symmetric_matrix<double>& A, boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::upper>& U);
}


#endif

math.cpp

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


void math::cholesky_factorize(const boost::numeric::ublas::symmetric_matrix<double>& A, boost::numeric::ublas::triangular_matrix<double>& L)
{
    namespace ublas = boost::numeric::ublas;

    ublas::matrix<double, ublas::column_major> CA(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) {
            CA(i, j) = A(i, j);
        }
    }

    info = boost::numeric::bindings::lapack::potrf('L', CA);
    BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic());

    L = ublas::triangular_adaptor<ublas::matrix<double, ublas::column_major> >(CA);
}


void math::cholesky_factorize(const boost::numeric::ublas::symmetric_matrix<double>& A, boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::upper>& U)
{
    namespace ublas = boost::numeric::ublas;

    ublas::matrix<double, ublas::column_major> CA(A.size1(), A.size2());
    int                                        info;

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

    info = boost::numeric::bindings::lapack::potrf('U', CA);
    BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic());

    U = ublas::triangular_adaptor<ublas::matrix<double, ublas::column_major>, ublas::upper>(CA);
}

math.cpp に正値対称行列のコレスキー分解を計算する関数 cholesky_factorize を実装しました.math.cpp をプロジェクトに加え(他のソースコードと一緒にコンパイル),math.hpp を include することにより cholesky_factorize が使用可能になります.
cholesky_factorize は N 次正値対称行列 A を受け取り,これを N 次下三角行列 L を用いて A = L x trans(L),あるいは N 次上三角行列 U を用いて A = trans(U) x U に分解します.どちらに分解するかはオーバーロードにより選ぶことができます.


以下のテストコードで cholesky_factorize の動作を確認できます.

test.cpp

#include "math.hpp"
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/triangular.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);

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

    ublas::triangular_matrix<double> L;
    ublas::triangular_matrix<double, ublas::upper> U;

    math::cholesky_factorize(A, L);
    math::cholesky_factorize(A, U);

    // 以下結果の表示
    std::cout << L << std::endl;
    std::cout << U << std::endl;

    std::cout << ublas::prod(L, ublas::trans(L)) << std::endl;
    std::cout << ublas::prod(ublas::trans(U), U) << std::endl;

    BOOST_UBLAS_CHECK(ublas::detail::expression_type_check(A, ublas::prod(L, ublas::trans(L))), ublas::internal_logic()); // A = L * L^T ?
    BOOST_UBLAS_CHECK(ublas::detail::expression_type_check(A, ublas::prod(ublas::trans(U), U)), ublas::internal_logic()); // A = U^T * U ?

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

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