Boost.uBLAS で逆行列

boost C++ 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乗）

それでは．