Boost.uBLAS でQR分解 - yanoの日記

前回の特異値分解に続いて，今回は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;
}

VC9
g++ (GCC) 3.2.2
g++ (GCC) 3.4.4 (cygming special)
g++ (GCC) 4.1.3

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

次回は対称行列の固有値分解を実装する予定です．
それでは．