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