読者です 読者をやめる 読者になる 読者になる

Boost.uBLAS で逆行列

Boost.uBLASを使い始めて,一番最初に疑問に思うことの多くは
「uBLASには逆行列を計算する関数はないのか」
だと思います.
残念ながらuBLASにはそのような関数はありません.
しかし,多くの先人が既にこの問題を解決してくれています.
具体的には次のサイトを参考にすればよいと思います.

以下に私が実装した例を紹介します.といっても,上記のサイトのコードに変更を少し加えただけですが.

math.hpp

#ifndef MATH_HPP_20080914
#define MATH_HPP_20080914

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


#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>


namespace math {
    template <class M, class MI> void invert(const M& m, MI& mi);
}


template <class M, class MI>
void math::invert(const M& m, MI& mi)
{
    namespace ublas = boost::numeric::ublas;

    BOOST_UBLAS_CHECK(m.size1() == m.size2(), ublas::external_logic());

    ublas::matrix<double>       lhs(m);
    ublas::matrix<double>       rhs(ublas::identity_matrix<double>(m.size1()));
    ublas::permutation_matrix<> pm(m.size1());

    ublas::lu_factorize(lhs, pm);
    ublas::lu_substitute(lhs, pm, rhs);

    BOOST_UBLAS_CHECK(rhs.size1() == m.size1(), ublas::internal_logic());
    BOOST_UBLAS_CHECK(rhs.size2() == m.size2(), ublas::internal_logic());

#if BOOST_UBLAS_TYPE_CHECK
    BOOST_UBLAS_CHECK(
        ublas::detail::expression_type_check(
            ublas::prod(m, rhs),
            ublas::identity_matrix<typename M::value_type>(m.size1())
        ),
        ublas::internal_logic()
    );
#endif

    mi.resize(rhs.size1(), rhs.size2(), false);
    mi.assign(rhs);
    // mi.assign_temporary(rhs);
}


#endif

math.hpp に逆行列を求める関数 invert を実装しました.
math.hpp を include するだけで invert が使用可能になります.
デバッグ時には実際に逆行列を乗算し,単位行列になるかをチェックするコードを入れています.
何故 mi = rhs とせずに assign を呼び出すかは

で noalias や assign が出てくる付近に詳しく書かれています(効率の問題です).
MI の型が ublas::matrix である場合は,assign の代わりに assign_temporary を使うことができます.
assign_temporary は実際には単なる swap です.よって,assign よりも高速になる可能性があります.
また,assign_temporary を使用する場合は,resize する必要もありません.

以下のテストコードで invert の動作を確認してみます.

test.cpp

#include "math.hpp"
#include <boost/numeric/ublas/io.hpp>


int main()
{
    boost::numeric::ublas::matrix<double> m(3, 3), mi;

    m(0, 0) = 1; m(0, 1) =  2; m(0, 2) =  3;
    m(1, 0) = 2; m(1, 1) =  1; m(1, 2) = -3;
    m(2, 0) = 4; m(2, 1) = -3; m(2, 2) =  2;

    math::invert(m, mi); // mi => m^{-1}

    std::cout << mi << std::endl;

    std::cout << boost::numeric::ublas::prod(m, mi) << std::endl; // 単位行列になるはず
    std::cout << boost::numeric::ublas::prod(mi, m) << std::endl; // 単位行列になるはず

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

でのコンパイルを確認しています.
実行してみると,正しく逆行列が求められていることが分かります.


余談ですが,DirectXでuBLASを用いるには少し注意が必要です.
というのも,D3DのCreateDeviceは行儀が悪く,FPUの設定を変えてしまうため,
BOOST_UBLAS_CHECKが通らない恐れがあるからです.
もしそうなったら,D3DCREATE_FPU_PRESERVE を検討してみてください.


さて今後についてですが,今のところ以下を予定しています.

  • math.hpp
    • determinant の実装
  • math.cpp
  • numeric.hpp
    • Expression Template による遅延評価な
      • 外積(cross product, vector product)
      • トレース(trace)
      • abs, log, sqrt などをベクトル,行列の全要素に適用
      • sum(ベクトルには定義されているが,行列には何故か定義されていないため)
      • norm_2_sq(ユークリッドノルムの2乗)
      • norm_frobenius_sq(フロベニウスノルムの2乗)

それでは.