手抜きcholesky分解
結局、ublasのlu_factorize関数周りの疑問はイマイチすっきりしないのでcholesky分解アルゴリズムを手で書いちゃうことに。あんまり難しいこと考えなければ、結構簡単なんだな。(トートロジーです)
/*! g++ main.cpp */ #include <boost/numeric/ublas/matrix.hpp> #include <boost/numeric/ublas/triangular.hpp> #include <boost/numeric/ublas/banded.hpp> #include <boost/numeric/ublas/io.hpp> #include <string> #include <iostream> #include <sstream> using namespace std; using namespace boost::numeric; /** * @brief Cholesky factorize m. It is assumed that m is positive definite symmetric square matrix and * non diagonal elements m(i, j), |i - j| > w are zero. No error check is done. */ template<class M> void cholesky(boost::numeric::ublas::matrix_expression<M>& m, const size_t w){ const size_t n = m().size1(); for (size_t i = 0; i < n; i++) { for (size_t j = i + 1; j < std::min(n, i + w + 1); j++) { m()(i, j) = m()(j, i) = m()(i, j) / m()(i, i); } for (size_t j = i + 1; j < std::min(n, i + w + 1); j++) { for (size_t k = i + 1; k < std::min(n, i + w + 1); k++) { m()(j, k) -= m()(j, i) * m()(i, k) * m()(i, i); } } //std::cout << m << std::endl; } } int main (int argc, char* argv[]) { try { ublas::matrix<double> m(3, 3); string temp = "[3, 3]((10, 1, 0), (1, 10, 1), (0, 1, 10))"; stringstream sst(temp); sst >> m; ublas::matrix<double> fact(m); cholesky(fact, 1); ublas::matrix<double> l = ublas::triangular_adaptor<ublas::matrix<double>, ublas::unit_lower>(fact); ublas::matrix<double> d = ublas::banded_adaptor<ublas::matrix<double> >(fact); ublas::matrix<double> hoge = ublas::prod(d, ublas::trans(l)); cout << ublas::prod(l, hoge) << endl; //これはコンパイルエラー //cout << ublas::prod(l, ublas::prod(d, ublas::trans(l))) << endl; cout << m << endl; } catch (std::exception& e) { cout << e.what() << endl; } return 0; }
相変わらず、ublasのtemplate expression周りが分かってないので、よく分からないエラーメッセージに悩まされる。。。*