boostのbanded_matrixの使い方が良く判らん

boost/numeric/ublasには、帯行列用のbanded_matrixクラスがある。数学的には、が帯行列ならばの計算はなはず。という訳で、実験してみたんだけど、どうも、banded_matrixを使ってもベクトルとの掛け算に掛かっている気がする。以下ソース。

#include <boost/numeric/ublas/banded.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>

#include <ctime>
#include <iostream>

using namespace boost::numeric;
using namespace boost;

using namespace std;



void makeBigMat(ublas::matrix<double>& m, const size_t n){
  m.resize(n, n);
}

void makeBigBandMat(ublas::banded_matrix<double>& bm, const size_t n,
		    const size_t l, const size_t u){
  bm.resize(n, n, l, u);
}

void makeBigVec(ublas::vector<double>& v, const size_t n){
  v.resize(n);
}

int main(int argc, char* argv[]){
  const int numSize  = atoi(argv[1]);
  const int sizeBy = atoi(argv[2]);
  const int numIter = atoi(argv[3]);

  

  vector<int> n(numSize);
  n[0] = 10;
  for(size_t i = 1; i < n.size(); i++){
    n[i] = sizeBy + n[i - 1];
  }


  
  for(size_t i = 0; i < n.size(); i++){
    ublas::matrix<double> mat;
    makeBigMat(mat, n[i]);

    ublas::banded_adaptor<ublas::matrix<double> > adp(mat, 1, 1);

    ublas::banded_matrix<double> band;
    makeBigBandMat(band, n[i], 1, 1);

    ublas::vector<double> v, ret(n[i]);
    makeBigVec(v, n[i]);

    clock_t start, end, normal, adaptor, banded;

    start = clock();
    for(size_t j = 0; j < numIter; j++){
      ublas::noalias(ret) = ublas::prod(mat, v);
    }
    end = clock();
    normal = end - start;



    start = clock();
    for(size_t j = 0; j < numIter; j++){
      ublas::noalias(ret) = ublas::prod(adp, v);
    }
    end = clock();
    adaptor = end - start;



    start = clock();
    for(size_t j = 0; j < numIter; j++){
      ublas::noalias(ret) = ublas::prod(band, v);
    }
    end = clock();
    banded = end - start;

    cout << n[i] << "\t"
	 << normal << "\t"
	 << adaptor << "\t"
	 << banded << endl;
  }
  

  return 0;
}

• ublas::noaliasは、一時オブジェクトを減らすためのおまじない。
• 最初、forの中身をublas::prod(mat, v); etc だけにしてたら、行列のサイズを変えても計算時間が全く変わらないのでびびる。これは、template expressionという技術で、中身が本当に必要になるまで計算を遅らせているからみたいである。boost恐るべし。
• デバッガで簡単に追ってみただけだと、banded_matrixだからと言って特別に変わっているような部分が見当たらないので、単にprod関数を使うだけじゃだめなのかも。

追記 : 解決した。http://d.hatena.ne.jp/hotoku/20110120/1295526130

最近勉強したこと、とか、最近の生活風景

• ノンパラ本終わり

Nonparametric Regression and Generalized Linear Models: A roughness penalty approach (Chapman & Hall/CRC Monographs on Statistics & Applied Probability)

一応最後まで読んだ。ひたすらスプラインな本。スプラインなんてたかが多項式と思っていたけど、こんなに出来るやつだったとは。で、読み終わったので、自分のデータをブチ込むべく中身のアルゴリズムを実装中。

• 行列代数

統計のための行列代数 上 (1)

一般逆行列とか特異値分解とか聞いたこと有るけどよく知らねーなと思っていたので理工系の嗜みとして購入。「行列代数」のタイトル通り、行列の線形写像の表示としての側面や、数字が2次元的に並んだ物という側面をなるべく表に出さないで、なるべく「行列」そのものを操作する代数でいろんな性質を語る本。代数だから、証明の表面は、抽象的な記号がずらずら並んだ等式の連鎖が多いんだけど、かなり細かく証明のステップが提示されているので、抽象的な代数の背後にある意味・幾何的な直観を補う作業が、通勤電車読書に丁度良い。

• xcode

職場では、windows + visual studio一本槍。だがしかし、大学院の研究ではunixツールが使えた方が何かと楽かと思って、mac book proを購入(今は、airにすればよかったと激しく後悔)。で、開発はxcodeでやろうと思ったんだけど、これがまた、やっぱり使い慣れたvisual studioからの垣根が滅茶苦茶高い。全然使いこなせてない。まぁそのうち慣れるだろうけど。