Boost.uBLAS で正値対称係数行列の線形方程式系(System of linear equations)の解の計算

boost C++ lapack uBLAS

今回は AX = B の解行列 X あるいは Ax = b の解ベクトル x を計算する関数を実装します．ただし，A は正値対称行列であるとします．
lapack には dposv という関数があり，これが内部で係数行列をコレスキー分解して解を計算してくれます．つまり，前回と同様 bindings.lapack.posv に渡す column_major な行列を用意するだけでほとんどの作業は終わりです．以下に私の実装例を示します．

math.hpp

#ifndef MATH_HPP_20081107
#define MATH_HPP_20081107

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


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


namespace math {
    void solve(const boost::numeric::ublas::symmetric_matrix<double>& A, const boost::numeric::ublas::vector<double>& b, boost::numeric::ublas::vector<double>& x);
    void solve(const boost::numeric::ublas::symmetric_matrix<double>& A, const boost::numeric::ublas::matrix<double>& B, boost::numeric::ublas::matrix<double>& X);
}


#endif

math.cpp

#include "math.hpp"
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/vector.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::solve(const boost::numeric::ublas::symmetric_matrix<double>& A, const boost::numeric::ublas::vector<double>& b, boost::numeric::ublas::vector<double>& x)
{
    BOOST_UBLAS_CHECK(A.size1() == b.size(), boost::numeric::ublas::external_logic());

    boost::numeric::ublas::matrix<double> B(b.size(), 1), X;

    boost::numeric::ublas::column(B, 0).assign(b);
    solve(A, B, X);

    x = boost::numeric::ublas::column(X, 0);
}


void math::solve(const boost::numeric::ublas::symmetric_matrix<double>& A, const boost::numeric::ublas::matrix<double>& B, boost::numeric::ublas::matrix<double>& X)
{
    namespace ublas = boost::numeric::ublas;

    BOOST_UBLAS_CHECK(A.size1() == B.size1(), ublas::external_logic());

    ublas::matrix<double, ublas::column_major> CA(A.size1(), A.size2()), CX(B);
    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::posv('L', CA, CX);
    BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic());

#if BOOST_UBLAS_TYPE_CHECK
    BOOST_UBLAS_CHECK(ublas::detail::expression_type_check(ublas::prod(A, CX), B), ublas::internal_logic());
#endif

    X = CX;
}

math.cpp に実正値対称係数行列の線形方程式系の解を計算する関数 solve を実装しました．デバッグ時には実際に計算した X を用いて，AX = B になることを確認するコードを追加しています．math.cpp をプロジェクトに加え（他のソースコードと一緒にコンパイル），math.hpp を include することにより solve が使用可能になります．
solve は AX = B の解行列 X あるいは Ax = b の解ベクトル x を計算します．ただし，A は N 次正値対称行列，X と B は N x NRHS 行列（NRHS は 1 以上の任意の数），x と b は N 次元ベクトルを表します．
NRHS が 1 のとき，X と B は N x 1 の行列になりますが，これは N 次元ベクトルに他なりません（NRHS が 1 のときは，AX = B という式はベクトルを用いて Ax = b になるのです）．uBLAS ではベクトルを表す型もちゃんと用意されているため，引数にベクトルを取る関数(Ax = b を解く関数)も実装しました．コードを見ると分かりますが，この関数は AX = B を解く関数を利用しているだけです．


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

test.cpp

#include "math.hpp"
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/vector.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(4, 4);
    ublas::vector<double>           b(4), x;

    A(0, 0) = 4.16;
    A(1, 0) = -3.12; A(1, 1) = 5.03;
    A(2, 0) = 0.56;  A(2, 1) = -0.83; A(2, 2) = 0.76;
    A(3, 0) = -0.1;  A(3, 1) = 1.18;  A(3, 2) = 0.34; A(3, 3) = 1.18;

    b(0) = 8.7; b(1) = -13.35; b(2) = 1.89; b(3) = -4.14;

    math::solve(A, b, x);

    // 以下結果の表示
    // 計算結果を http://www.nag.co.uk/lapack-ex/node11.html と比較してみましょう
    std::cout << x << std::endl;

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

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


次回は最小二乗法を実装する予定です．
それでは．