Boost.uBLAS でQR分解
前回の特異値分解に続いて,今回はQR分解を実装してみます.
特異値分解はあっけないほど単純に実装できたのですが,QR分解はそれよりも少し長くなります.
なぜなら,特異値分解は lapack.dgesvd(lapack.dgesdd) を呼び出すだけで良かったのですが,QR分解は lapack.dgeqrf と lapack.dorgqr を呼び出さなければいけないからです.
以下に私の実装例を示します.
math.hpp
#ifndef MATH_HPP_20080923 #define MATH_HPP_20080923 #if defined(_MSC_VER) && (_MSC_VER >= 1020) # pragma once #endif #include <boost/numeric/ublas/fwd.hpp> namespace math { void qr_factorize(const boost::numeric::ublas::matrix<double>& A, boost::numeric::ublas::matrix<double>& Q, boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::upper>& R); } #endif
math.cpp
#include "math.hpp" #include <vector> #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/geqrf.hpp> #include <boost/numeric/bindings/lapack/orgqr.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::qr_factorize(const boost::numeric::ublas::matrix<double>& A, boost::numeric::ublas::matrix<double>& Q, boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::upper>& R) { namespace ublas = boost::numeric::ublas; BOOST_UBLAS_CHECK(A.size1() >= A.size2(), ublas::external_logic()); std::vector<double> tau(A.size2()); ublas::matrix<double, ublas::column_major> CQ(A); int info; info = boost::numeric::bindings::lapack::geqrf(CQ, tau); BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic()); ublas::triangular_matrix<double, ublas::upper> CR = ublas::triangular_adaptor<const ublas::matrix_range<ublas::matrix<double, ublas::column_major> >, ublas::upper>( ublas::project(CQ, ublas::range(0, A.size2()), ublas::range(0, A.size2())) ); info = boost::numeric::bindings::lapack::orgqr(CQ, tau); BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic()); #if BOOST_UBLAS_TYPE_CHECK BOOST_UBLAS_CHECK(ublas::detail::expression_type_check(ublas::prod(CQ, CR), A), ublas::internal_logic()); #endif ublas::matrix<double> CCQ(CQ); Q.assign_temporary(CCQ); R.assign_temporary(CR); }
math.cpp にQR分解を計算する関数 qr_factorize を実装しました.デバッグ時には実際に計算結果を用いて QR が A に等しくなるかを確認するコードを追加しています.math.cpp をプロジェクトに加え(他のソースコードと一緒にコンパイル),math.hpp を include することにより qr_factorize が使用可能になります.
関数 qr_factorize は M x N 行列 A (ただし,M >= N) を受け取り,これを M x N 正規直交行列 Q と N x N 上三角行列 R の積に因数分解します.
M x N 正規直交行列が関数の引数 Q に,N x N 上三角行列が関数の引数 R に保存されます.
コードの簡単な説明ですが,bindings.lapack.geqrf を呼び出すと,第一引数の右上の要素にQR分解のRとなる上三角行列が格納されます.
M > N のときがあるので,第一引数である CQ をまず ublas.project により N x N 正方行列に切り出し,これの右上成分を ublas.triangular_adaptor で取り出します.
次に,bindings.lapack.orgqr を呼び出すと,第一引数がQR分解の Q になります.
これらの結果を関数の引数 Q, R に保存して終了です.
以下のテストコードで qr_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> m(3, 3); ublas::matrix<double> Q; ublas::triangular_matrix<double, ublas::upper> R; m(0, 0) = 12; m(0, 1) = -51; m(0, 2) = 4; m(1, 0) = 6; m(1, 1) = 167; m(1, 2) = -68; m(2, 0) = -4; m(2, 1) = 24; m(2, 2) = -41; math::qr_factorize(m, Q, R); // 以下結果の表示 // 計算結果を http://en.wikipedia.org/wiki/QR_decomposition#Example と比較してみましょう std::cout << Q << std::endl; std::cout << R << std::endl; std::cout << ublas::prod(Q, R) << std::endl; // 元に戻る(m に等しくなる)はず return 0; }
でのコンパイルを確認しています.
実行してみると,正しくQR分解できていることが分かります.
エラーが出てコンパイルできない場合は,前回の記事を見て正しくインストールできているか確認してみてください.
次回は対称行列の固有値分解を実装する予定です.
それでは.