手抜き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周りが分かってないので、よく分からないエラーメッセージに悩まされる。。。*