Boost.uBLAS で特異値分解(SVD) - 高速版

boost C++ lapack uBLAS

前回は bindings.lapack.gesvd による特異値分解を紹介しました．
lapack にはもう一つ特異値分解を行う関数があります．
それが今回紹介する bindings.lapack.gesdd(lapack.dgesdd) になります．
gesvd との違いは分割統治法(divide-and-conquer method)により処理を高速化したことです．
早速実装してみましょう（前回の math.hpp と test.cpp は変更を加える必要がありません）．

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/gesdd.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::gesdd('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);
}
  • VC9
  • g++ (GCC) 3.4.4 (cygming special)
  • g++ (GCC) 4.3.3

でのコンパイルを確認しています．
include するヘッダファイルを gesvd.hpp から gesdd.hpp に変更したのと，gesvd の呼び出しを gesdd に変更しただけです．
gesvd と比べどの程度高速化されるかは与える行列にかなり影響されるため，はっきりとした数字はいえませんが，私の経験ではだいたい -10% から 70% くらいの高速化になりました（遅くなる場合もかなり速くなる場合もあったということです）．
たいていの場合は gesvd に比べ速くなるとは思いますが，精度が悪くなる場合もあるので注意してください．


次回はQR分解について説明する予定です．
それでは．