Boost.uBLAS で特異値分解(SVD)

Boost.uBLAS と Boost.Bindings, clapack をインストールした今,特異値分解を簡単に実装することができます.
というのも,Boost.Bindings を通して bindings.lapack.gesvd(lapack.dgesvd) を呼び出すだけだからです.
以下に私の実装例を示します.

math.hpp

#ifndef MATH_HPP_20090501
#define MATH_HPP_20090501

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


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


namespace math {
    void svd(const boost::numeric::ublas::matrix<double>& A, boost::numeric::ublas::matrix<double>& U, boost::numeric::ublas::diagonal_matrix<double>& D, boost::numeric::ublas::matrix<double>& VT);
}


#endif

math.cpp

#include "math.hpp"
#include <algorithm>
#include <vector>
#include <boost/numeric/ublas/banded.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#define BOOST_NUMERIC_BINDINGS_USE_CLAPACK
#include <boost/numeric/bindings/lapack/gesvd.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::svd(const boost::numeric::ublas::matrix<double>& A, boost::numeric::ublas::matrix<double>& U, boost::numeric::ublas::diagonal_matrix<double>& D, boost::numeric::ublas::matrix<double>& VT)
{
    namespace ublas = boost::numeric::ublas;

    std::vector<double>                        s((std::min)(A.size1(), A.size2()));
    ublas::matrix<double, ublas::column_major> CA(A), CU(A.size1(), A.size1()), CVT(A.size2(), A.size2());
    int                                        info;

    info = boost::numeric::bindings::lapack::gesvd('A', 'A', CA, s, CU, CVT);
    BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic());

    ublas::matrix<double>          CCU(CU), CCVT(CVT);
    ublas::diagonal_matrix<double> CD(A.size1(), A.size2());

    for (std::size_t i = 0; i < s.size(); ++i) {
        CD(i, i) = s[i];
    }

#if BOOST_UBLAS_TYPE_CHECK
    BOOST_UBLAS_CHECK(
        ublas::detail::expression_type_check(ublas::prod(ublas::matrix<double>(ublas::prod(CCU, CD)), CCVT), A),
        ublas::internal_logic()
    );
#endif

    U.assign_temporary(CCU);
    D.assign_temporary(CD);
    VT.assign_temporary(CCVT);
}

math.cpp に特異値分解を計算する関数 svd を実装しました.デバッグ時には実際に計算結果を用いて U x D x VT が A に等しくなるかを確認するコードを追加しています.
関数 svd は M x N 行列 A を受け取り,これを M x M ユニタリ行列 U と M x N 対角行列 D と N x N ユニタリ行列 V の転置の積に因数分解します.M x M ユニタリ行列が関数の引数 U に,M x N 対角行列が関数の引数 D に,N x N ユニタリ行列 V の転置が関数の引数 VT に保存されます.その結果,関数の引数 A は A = U x D x VT と分解されます.
M x N 対角行列 D の対角成分は A の特異値になっており,非負の実数で,降順に格納されています.また,M x N 対角行列 D の対角成分は min(M, N) 個であることに注意してください.
math.cpp をプロジェクトに加え(他のソースコードと一緒にコンパイル),math.hpp を include することにより svd が使用可能になります.

今回は clapack を使用しているため,bindings のヘッダファイルを読む前に BOOST_NUMERIC_BINDINGS_USE_CLAPACK を定義し,bindings に clapack を使用していることを告げます.これにより bindings は正しいシンボルを見つけることができます.

さらに注意が必要なのは,lapackFORTRAN で書かれており,FORTRAN と C(C++) はメモリレイアウトが異なるということです.これは lapack を機械的に f2c した clapack にも当てはまります.
つまり,C の行列を FORTRAN(lapack) にそのまま渡すことはできません.lapack に行列を渡す前に,行列のメモリレイアウトを列優先に変換しなければなりません.
これは少し面倒な作業に聞こえますが,uBLAS には column_major というまさにこのための構造があります.
これを用いることにより,FORTRAN のメモリレイアウトに合わせることができます.

lapack.gesvd の後は column_major な行列を関数の row_major な引数に結果として保存するだけです.
ここで,U, D, VT を resize して assign するのではなく,CCU, CD, CCVT をコンストラクトして assign_temporary しているのは強い例外安全保証(Strong exception-safety guarantee)を考慮してのことです.つまり,関数 svd は強い例外安全を保証するため,処理に失敗してもリソースをリークすることはなく,かつ関数の引数の状態も呼び出し前から変更されることは決してありません.

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

test.cpp

#include "math.hpp"
#include <boost/numeric/ublas/banded.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix.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(4, 5);

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

    ublas::matrix<double> U, VT;
    ublas::diagonal_matrix<double> D;

    math::svd(m, U, D, VT); // m = U x D x VT

    // 以下 U, D, VT の表示
    // 計算結果を http://en.wikipedia.org/wiki/Singular_value_decomposition#Example と比較してみましょう
    std::cout << U << "\n";
    std::cout << D << "\n";
    std::cout << VT << std::endl;

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

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


次回は分割統治法を用いたより高速な特異値分解を実装する予定です.
とはいっても,bindings.lapack.gesvd の代わりに bindings.lapack.gesdd を使うだけですが.
それでは.