Boost.uBLAS でLU分解 - lapack版
今回は前回実装したLU分解を lapack を用いて実装してみます．
前回は正方行列に限定してLU分解しましたが，今回は PA = LU の形に分解するため，任意の行列をLU分解できます．
その代償としてLU分解後の行列の切り出しの処理が必要になり，コード量が増えてしまいますが．
以下に私の実装例を示します．
math.hpp
#ifndef MATH_HPP_20081013 #define MATH_HPP_20081013 #if defined(_MSC_VER) && (_MSC_VER >= 1020) # pragma once #endif #include <boost/numeric/ublas/fwd.hpp> #include <boost/numeric/ublas/lu.hpp> namespace math { void lu_factorize(const boost::numeric::ublas::matrix<double>& A, boost::numeric::ublas::permutation_matrix<>& P, boost::numeric::ublas::triangular_matrix<double>& L, 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/matrix_proxy.hpp> #include <boost/numeric/ublas/triangular.hpp> #define BOOST_NUMERIC_BINDINGS_USE_CLAPACK #include <boost/numeric/bindings/lapack/gesv.hpp> #include <boost/numeric/bindings/traits/ublas_matrix.hpp> #include <boost/numeric/bindings/traits/ublas_vector.hpp> #undef BOOST_NUMERIC_BINDINGS_USE_CLAPACK void math::lu_factorize(const boost::numeric::ublas::matrix<double>& A, boost::numeric::ublas::permutation_matrix<>& P, boost::numeric::ublas::triangular_matrix<double>& L, boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::upper>& U) { namespace ublas = boost::numeric::ublas; ublas::permutation_matrix<int> CP((std::min)(A.size1(), A.size2())); ublas::matrix<double, ublas::column_major> CA(A); int info; info = boost::numeric::bindings::lapack::getrf(CA, CP); BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic()); ublas::permutation_matrix<> CCP(CP.size()); for (std::size_t i = 0; i < CCP.size(); ++i) { CCP(i) = CP(i) - 1; } ublas::triangular_matrix<double> CL = ublas::triangular_adaptor<const ublas::matrix_range<ublas::matrix<double, ublas::column_major> > >( ublas::project(CA, ublas::range(0, A.size1()), ublas::range(0, CP.size())) ); for (std::size_t i = 0; i < CL.size2(); ++i) { CL(i, i) = 1; } ublas::triangular_matrix<double, ublas::upper> CU = ublas::triangular_adaptor<const ublas::matrix_range<ublas::matrix<double, ublas::column_major> >, ublas::upper>( ublas::project(CA, ublas::range(0, CP.size()), ublas::range(0, A.size2())) ); P.assign_temporary(CCP); L.assign_temporary(CL); U.assign_temporary(CU); }
math.cpp にLU分解を計算する関数 lu_factorize を実装しました．math.cpp をプロジェクトに加え（他のソースコードと一緒にコンパイル），math.hpp を include することにより lu_factorize が使用可能になります．
関数 lu_factorize は M x N 行列 A を受け取り，これを PA = LU の形に因数分解します．ただし，P は M x M 置換行列(Permutation matrix)，L は M x min(M, N) 単位下三角行列，U は min(M, N) x N 上三角行列を表します．そのため，M > N の場合 L は下台形行列，M < N の場合 U は上台形行列になることに注意してください．
M x M 置換行列が関数の引数 P に，M x min(M, N) 単位下三角行列が関数の引数 L に，min(M, N) x N 上三角行列が関数の引数 U に保存されます．
コードの簡単な説明ですが，bindings.lapack.getrf を呼び出すと，第一引数に下三角行列と上三角行列が，第二引数に置換行列が格納されます．
lapack が返す置換行列はインデックスが 1 から始まっているため，これを uBLAS の permutation_matrix に合わせるため 1 を引いています．
その後，M x min(M, N) 下三角行列と min(M, N) x N 上三角行列を project と triangular_adaptor で切り出します．そして下三角行列の対角成分を 1 にします．
これらの結果を関数の引数 P, L, U に保存して終了です．
今回のプログラムはただ lapack に処理を投げるだけでなく，for で回したり size1 や size2 が入り乱れていて，このコードが本当に正しいのか疑問をもたれるかもしれませんが，今回の実装例に限らず本ブログに載せているコードはすべて実際にコンパイルし様々な入力を与え出力を確認していますので，どうか安心して使ってください．同様に，説明文中の M や N といった間違い易そうな箇所についても，何回も確認し間違いがないように最新の注意を払っています．それでももし間違いじゃないかという部分を，あるいは間違いを発見された場合には，是非知らせてください．
以下のテストコードで lu_factorize の動作を確認できます．
test.cpp
#include "math.hpp" #include <boost/numeric/ublas/io.hpp> #include <boost/numeric/ublas/matrix.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::matrix<double> A(3, 5); A(0, 0) = 1; A(0, 1) = 2; A(0, 2) = 3; A(0, 3) = 4; A(0, 4) = 5; A(1, 0) = 2; A(1, 1) = 3; A(1, 2) = 4; A(1, 3) = 5; A(1, 4) = 6; A(2, 0) = 3; A(2, 1) = 4; A(2, 2) = 5; A(2, 3) = 6; A(2, 4) = 7; ublas::permutation_matrix<> pm(0); ublas::triangular_matrix<double> L; ublas::triangular_matrix<double, ublas::upper> U; math::lu_factorize(A, pm, L, U); // 以下結果の確認 ublas::matrix<double> P(ublas::identity_matrix<double>(A.size1())); ublas::swap_rows(pm, P); std::cout << P << std::endl; std::cout << L << std::endl; std::cout << U << std::endl; std::cout << ublas::prod(P, A) << std::endl; std::cout << ublas::prod(L, U) << std::endl; // PA = LU ? BOOST_UBLAS_CHECK(ublas::detail::expression_type_check(ublas::prod(P, A), ublas::prod(L, U)), ublas::internal_logic()); return 0; }
でのコンパイルを確認しています．
実行してみると，正しくLU分解できていることが分かります．
エラーが出てコンパイルできない場合は，前回の記事を見て正しくインストールできているか確認してみてください．
