統計のための線形代数 in C++ (Boost uBLASの参考サイトとtipsまとめ)

最近の僕はより大規模な実験を行うべく、これまで作ってきたtoy programをC++に移植する毎日です。
その過程でどうしてもベクトル・行列演算が必要になったので、時期標準の呼び声高いBoost内のベクトル・行列ライブラリ、uBLASを用いてみました。
今回の記事は、特に僕が必要であった演算を実現すべく、参考にしたサイト及びそれらをまとめたテンプレート関数群をご紹介します。
参考にした主なサイトは以下のとおりです。

1. Boost Basic Linear Algebra（本家）
2. Boost 数学関連ライブラリの使い方
3. yanoの日記

2で概要を掴み、1で詳細を詰め、3で必要なものを実装... がオススメです。
特に3の矢野さんの日記には大変助けられました。ほとんど全ての実装がblogを元にしているといっても過言ではありません（この場を借りて御礼申し上げます）。

先人達の知恵を集結し、まとめた関数群がこちら。

/*
 matrix.hpp is template functions using uBLAS
 References:
   boost Basic Linear Algebra http://www.boost.org/doc/libs/1_39_0/libs/numeric/ublas/doc/index.htm
  yanoの日記 http://d.hatena.ne.jp/blono/
  cholesky.hpp http://www.guwi17.de/ublas/examples/cholesky.hpp
  
 ver.0.1
 13/10/09
 Takaaki
*/

#ifndef _MATRIX_ 
#define _MATRIX_

#include <boost/numeric/ublas/lu.hpp>   
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <algorithm>

using namespace boost::numeric::ublas;

namespace math{
 namespace ublas = boost::numeric::ublas;

 /*mean*/
 template <class V> double mean(const V& v);

 /*variance*/
 template <class V> double var(const V& v);

 /*standard deviation*/
 template <class V> double sd(const V& v);

 /*Invert*/
 template <class M, class MI> void invert(const M& m, MI& mi); 

 /*Trace*/
 template <class M> double trace(const M& m);

 /*Determinant*/
 template <class M> double determinant(const M& m);
 template <class M> double det(const M& m);

 /*Cholesky decomposition*/
 template <class M> size_t cholesky(const M& A, M& L);
}

template <class V>
double math::mean(const V& v)
{
 return ublas::sum(v) / v.size();
}

template <class V>
double math::var(const V& v)
{
 /*Create mean vector*/
 ublas::vector<double> mv(v.size());
 std::fill(mv.begin(), mv.end(), math::mean(v));

 return sum(element_prod(v - mv, v-mv)) / (v.size() - 1.0);
}

template <class V>
double math::sd(const V& v)
{
 return std::sqrt(var(v));
}

template <class M, class MI>
void math::invert(const M& m, MI& mi)
{

 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);

 mi.resize(rhs.size1(), rhs.size2(), false);
 mi.assign(rhs);
}

template <class M>
double math::trace(const M& m)
{
 BOOST_UBLAS_CHECK(m.size1() == m.size2(), ublas::external_logic());

 M cm(m);
 ublas::matrix_vector_range<M> diag(cm, ublas::range(0,m.size1()), ublas::range(0,m.size2()));
 
 return ublas::sum(diag);
}

template <class M>
double math::determinant(const M& m)
{
 BOOST_UBLAS_CHECK(m.size1() == m.size2(), ublas::external_logic());

 ublas::matrix<double> lu(m);
 ublas::permutation_matrix<> pm(m.size1());

 ublas::lu_factorize(lu,pm);

 double det(1);

 typedef ublas::permutation_matrix<>::size_type size_type;
 
 for (size_type i = 0; i < pm.size(); ++i) {
  det *= (i == pm(i)) ? +lu(i, i) : -lu(i, i);
 }

 return det;
}

template <class M>
double math::det(const M& m)
{
 return math::determinant(m);
}

template <class M>
size_t math::cholesky(const M& A, M& L)
{

 typedef typename M::value_type T;

 assert( A.size1() == A.size2() );
 assert( A.size1() == L.size1() );
 assert( A.size2() == L.size2() );

 const size_t n = A.size1();
         
 for (size_t k=0 ; k < n; k++) {
                 
   double qL_kk = A(k,k) - inner_prod( project( row(L, k), range(0, k) ),project( row(L, k), range(0, k) ) );

  if (qL_kk <= 0) {
   return 1 + k;
  } else {
   double L_kk = sqrt( qL_kk );
   L(k,k) = L_kk;
   matrix_column<M> cLk(L, k);
   project( cLk, range(k+1, n) )
    = ( project( column(A, k), range(k+1, n) )
     - prod( project(L, range(k+1, n), range(0, k)), 
      project(row(L, k), range(0, k) ) ) ) / L_kk;
  }
 }
 return 0;      
}

#endif

これをmath.hppとして保存し、以下のmain.cppで動作を確認してみます。

#include <cstdlib>
#include "math.hpp"

/*function of random variables*/
double runif();
double rnorm(double, double);

const static double PI = 3.14159265358979323846;

/*main*/
int main(void)
{
 srand((unsigned int)time(NULL));

 int dim = 3;
 double mean = 0.0;
 double sd = 1.0;

 /*create vector & matrix*/
 vector<double> v(dim);
 matrix<double> m(dim,dim); /*synmetric matrix*/
 matrix<double> im(dim,dim); /*inverse of m*/
 matrix<double> L(dim,dim); /*triangular matrix*/

  /*fill the vector & matrix as diagnal matrix*/
 for(int i=0;i<dim;++i){
  v(i) = static_cast<double>(i + 3);
  for(int j=0;j<dim;++j){
   if(j >= i){
    i == j ? m(i,j) = 1.0 : m(i,j) = rnorm(mean,sd);
   }else{
    m(i,j) = m(j,i);
   }
  }
 }

 /*print*/
 std::cout<<"v : "<<v<<"\n";
 std::cout<<"m : "<<m<<"\n";

 /*mean, variance, standard deviation*/
 std::cout<<"mean(v) : "<<math::mean(v)<<"\n";
 std::cout<<"var(v) : "<<math::var(v)<<"\n";
 std::cout<<"sd(v) : "<<math::sd(v)<<"\n";

 /*inverse, trace, determinant, cholesky */
 math::invert(m,im);
 std::cout<<"m * m^-1 : "<<prod(m,im)<<"\n"; /*check the inverse*/
 std::cout<<"trace(m) : "<<math::trace(m)<<"\n";
 std::cout<<"determinant(m) : "<<math::determinant(m)<<"\n";
 math::cholesky(m,L);
 std::cout<<"cholesky decompsition of m\n";
 std::cout<<"L : "<<L<<"\n";
 std::cout<<"L' : "<<trans(L)<<"\n";
 std::cout<<"L * L' : "<<prod(L,trans(L))<<"\n";

 return 0;
}


double runif()
{
 return static_cast<double>(rand()) / RAND_MAX;
}

double rnorm(double mean, double sd)
{
 double U=0.0;
 double V=0.0;
 while(true){
  U = runif();
   V = runif();
  if(U>0.000001 && V > 0.000001){
   break;
  }
 }

 double R = sqrt(-2.0*log(U));
 double T = 2.0*PI*V;

 if(U>0.5){
  return sd * R * cos(T) + mean;
 }else{
  return sd * R * sin(T) + mean;
 }
}

確認環境はg++(gcc) 4.1.2、boost 1_39_0です。

$g++ main.cpp
$./a.out

の出力結果がこちら。

v : [3](3,4,5)
m : [3,3]((1,0.30186,0.0902471),(0.30186,1,0.0612414),(0.0902471,0.0612414,1))
mean(v) : 4
var(v) : 1
sd(v) : 1
m * m^-1 : [3,3]((1,-8.23994e-18,0),(-3.72966e-17,1,-1.38778e-17),(0,6.93889e-18,1))
trace(m) : 3
determinant(m) : 0.900322
cholesky decompsition of m
L : [3,3]((1,0,0),(0.30186,0.953352,0),(0.0902471,0.035663,0.995281))
L' : [3,3]((1,0.30186,0.0902471),(0,0.953352,0.035663),(0,0,0.995281))
L * L' : [3,3]((1,0.30186,0.0902471),(0.30186,1,0.0612414),(0.0902471,0.0612414,1))

Awesome!