この記事は，Boost Advent Calendar 2011の14日目の記事です．
(2時間近くも遅刻してすみません＞＜)

スキルはまだまだですが，それでももうboostが無いと生きていけない体になってしまいました．
そういえば魔導書Vol2のid:kikairoya先生の章(一番最後ですね)あたりは，ホントに「boostを使い倒す！」かんじがして最高に読み応えがありました．ぜひお読みください．

さて，本来はこのアドベでは無難に「あまり日本語情報の出回ってないライブラリを(ほとんど翻訳で)紹介する」つもりだったんですが，意外とネタ探しにつまづいたので，方針転換して僕が引っかかってしまった落とし穴を紹介することに．
筋違いなことを言ってたら，たぶん言ってますが，ツッコミをお願いします．
ホントに不当なdisがあったら全力で謝罪記事を書かせて頂きます…．

(2011.12.15 12:30) 編集上のミスについて指摘をいただいたので訂正

※正確にはboost::numeric::ublasですが，便宜的にboost::ublasと表記します．

boost::ublasとは

特に細かい説明の必要は無いとは思いますが．
1.29.0で正式にboost入りした，線形代数ライブラリです．
BLASというメジャーな線形代数APIがあって(ライブラリというよりは規格と言ったほうが正確かもしれません)，このサブセットとしてuBLASがあります．読みは「ゆーぶらす」．

#include <iostream>
#include <algorithm>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

using namespace boost::numeric::ublas;

int main()
{
	vector<double> vec(3);
	vec[0] = 1;
	vec[1] = 2;
	vec[2] = 3;

	matrix<double> mat(3, 2);
	mat(0, 0) = 1; mat(0, 1) = 2;
	mat(1, 0) = 3; mat(1, 1) = 4;
	mat(2, 0) = 5; mat(2, 1) = 6;

	std::cout << vec << std::endl;
	std::cout << mat << std::endl;
	std::cout << (3 * vec) << std::endl;
	std::cout << prod(vec, mat) << std::endl;
}

に対して

 [3](1,2,3)
 [3,2]((1,2),(3,4),(5,6))
 [3](3,6,9)
 [2](22,28)

なかんじ．

uBLASを紹介することは本筋じゃないので，このあたりはさらっと．

ublas::bounded_vectorについて

さてさて．
みなさんは自分のプログラムでちょっとした線形代数の数値計算をするとき，だいたい何次元の行列なりベクトルを利用されますか？
データ処理や機械学習の分野では，数万から数百万次元というレベルがわりとザラです．天気予報やX線CTなんかもあれは中で超高次元な一次連立方程式を解いてるだけです．
世に出回ってる線形代数ライブラリは基本的にはこのようなタスクに最適化されています(特にC系のものは)．

でもたぶん普通の(普通のって何よって話ですが)プログラム書く上で使うとしたら，平面ベクトルとしての2次元から空間のアフィン変換としての4×4次元というごく低次元かつ固定次元な場合がほとんどでしょう．こう低次元ともなると，データ構造やアルゴリズムも高次元用のものでは非常に非効率になります．ちっさい領域が内部でうじゃうじゃと動的メモリ確保なんかされるとなんかもうね．
そんなとき性能面や扱いやすさの面で便利なのが，固定次元ベクトル型ublas::bounded_vectorです(次元ごとの最適な数値計算アルゴリズムの話は置いとくとして)．

• ublas::bounded_vector
• 次元はテンプレート引数Nで与える(コンパイル時に固定！)
• 内部で固定長配列を用いる(動的メモリ確保を行わない！)

という特徴を持っていて，

• 低次元かつ固定次元という目的にぴったり
• 十分次元が低ければ，性能および扱いの面でプリミティブ型の感覚で使える

という利点があり，かなりガンガン使えます．
もちろん，通常のublas::vectorと混在して使うことができます．

実際に少しだけ使ってみた例．

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

using namespace boost::numeric::ublas;

int main()
{
	//内部で固定長配列なので，使う領域は次元数×sizeof(T)+4(次元数も保管するため)
	std::cout << sizeof(bounded_vector<int, 4>) << std::endl;
	std::cout << sizeof(bounded_vector<float, 4>) << std::endl;
	std::cout << sizeof(bounded_vector<double, 4>) << std::endl;
	std::cout << sizeof(bounded_vector<double, 8>) << std::endl;

	//内部で固定長配列なので，
	//vector<char> huge1(16777216);			//ヒープでならなんでもない量(16MB)でも
	//bounded_vector<char, 16777216> huge2;		//スタックに配置されてるとsegfault※スタックサイズ依存

	//bounded_vectorと
	bounded_vector<double, 3> vec1;
	vec1[0] = 1; vec1[1] = 2; vec1[2] = 3;

	//ふつうのvectorの
	vector<double> vec2(3);
	vec2[0] = 1; vec2[1] = 2; vec2[2] = 3;

	//混在した演算も余裕
	std::cout << (vec1 + vec2) << std::endl;

	//というかクラス定義を見たらなーんだってなります
}

で

 20
 20
 36
 68
 [3](2,4,6)

いいかんじですね．低次元固定次元で使いやすそう．
でもね，使うときはちょっと気をつけて欲しいんです．

ちょっと設計が中途半端な気もする！

ublas::bounded_vectorは，次のコードのコンパイルに成功し，実行時に例外を吐きます．

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

using namespace boost::numeric::ublas;

int main()
{
	bounded_vector<double, 3> vec3;
	bounded_vector<double, 4> vec4;

	std::cout << (vec3 + vec4) << std::endl;
}

Check failed in file /usr/local/include/boost/numeric/ublas/vector_expression.hpp at line 548:
size1 == size2
terminate called after throwing an instance of 'boost::numeric::ublas::bad_argument'
  what():  bad argument
zsh: abort      ./a.out

おいちょっと待てと．
なんで．コンパイル時にエラーにさせればいいじゃん！できるじゃん！

まぁこれだけなら気をつけてればで片付けられるかなって気もします．

で，次のコードはコンパイルが通り例外も出ません．

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

using namespace boost::numeric::ublas;

int main()
{
	bounded_vector<double, 3> vec3;
	bounded_vector<double, 4> vec4;

	vec4 = vec3;

	std::cout << vec3 << std::endl;
	std::cout << vec4 << std::endl;
}

[3](-1.08785e-41,-4.66095e-42,-1.0874e-41)
[3](-1.08785e-41,-4.66095e-42,-1.0874e-41)

やーん＞＜
次数が違うのに代入できるの＞＜
vec4の次数が3になってる＞＜

一方で，

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

using namespace boost::numeric::ublas;

int main()
{
	bounded_vector<double, 3> vec3;
	bounded_vector<double, 4> vec4;

	vec3 = vec4;	//代入を逆にしてみた

	std::cout << vec3 << std::endl;
	std::cout << vec4 << std::endl;
}

Check failed in file /usr/local/include/boost/numeric/ublas/storage.hpp at line 343:
size <= N
terminate called after throwing an instance of 'boost::numeric::ublas::bad_size'
  what():  bad size
zsh: abort      ./a.out

これは落ちますか…

でも↓こう書いちゃえば，高次元→低次元もいけるし，異なる次元の演算もいけます．

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

using namespace boost::numeric::ublas;

int main()
{
	bounded_vector<double, 3> vec3;
	bounded_vector<double, 4> vec4;

	vec4.resize(3);

	vec3 = vec4;
	std::cout << (vec3 + vec4) << std::endl;
	std::cout << vec3 << std::endl;
	std::cout << vec4 << std::endl;
}

[3](-2.38555e-41,-0.578328,4.37691e-314)
[3](-1.19277e-41,-0.289164,2.18846e-314)
[3](-1.19277e-41,-0.289164,2.18846e-314)

えーでもなんかなぁなんかなぁ＞＜；；

この辺りってのは設計思想からすると当然と言えば当然で，要はbounded_vectorは「ストレージを長さNの固定長配列とすることで，最大N次元のベクトルを表現できるようにした」ものであって，「次元を固定値Nにした」というわけではないのですねー．

うーん，でもやっぱり

bounded_vector<double, 4> vec;

て書いてあったら，vecは4次元の固定長ベクトルだしそれに対して3次元ベクトルを代入したりしたらコンパイル時でエラーになるんだろうと誰でも期待しちゃう気がします．
理解が浅いなら使うなっていう話で型付け片付けちゃってもいいんですが．

というか，一つの(実質プリミティブな)オブジェクトを3次元には使うわ4次元には使うわ，その時点で習慣としては良くないし，そのうえ(上の例では)5次元以上からは例外吐くのでアウトとか，やっぱり使いにくさの元になるんじゃないでしょうか．
「使い回しは悪しき習慣」てことでconst教の布教につなげることもできなくはないですね．

これtypedefとか関数渡しとかでちょっと複雑になると，ずいぶんとわかりにくいバグになります．
次のプログラムはいろんな意味で全く無意味ですが，ちょっとデバッガ無しでどこで落ちるか当てらんないんじゃないでしょうか．
ここではわざとらしく間違えてるのですが，代入に代入を経てそのオブジェクトは今何次元ベクトルなのってのが曖昧な状態でやらかしてしまいます(やらかしてしまいました)．

//どこで落ちるでしょう
#include <iostream>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
using namespace boost::numeric::ublas;

typedef bounded_vector<double, 3> vec3d;
typedef bounded_vector<double, 4> vec4d;

vec4d unit(vec4d v)
{
	double norm = norm_2(v);
	return v / norm;
}

double angle(vec4d v1, vec4d v2)
{
	v1 = unit(v1);
	v2 = unit(v2);
	return acos(inner_prod(v1, v2));
}

int main()
{
	vec3d v1;
	v1[0] = 1; v1[1] = 1; v1[2] = 0;

	vec4d v2;
	v2[0] = 0; v2[1] = 1; v2[2] = 1;

	std::cout << angle(v1, v2) << std::endl;
}

そう，アホなんですよ，こんなしょーもないことでつまづいた僕がアホなんです．
でもなんとなくそれで済ませたくなかった…

夢を見たいならEigenでどうぞ

ここで，線形代数ライブラリEigenの紹介（ちょ
uBLASと同じくバリバリのC++で，ヘッダオンリーなライブラリです．
行列やベクトルの性質(次数やsparseかdenseか，などいろいろ)に併せたデータ構造とアルゴリズムを持つこと，低次元・固定次元のためにプリミティブな型の感覚で使えるベクトル型をもつことはuBLASと同じです．
さらに

• 真に型レベルで固定長
  • 違う次数への代入や演算なんかは全部型レベルでハジく
• 行列も固定次元にできる
  • というかベクトル自体が1行または1列の行列
  • つまり行ベクトルも列ベクトルも表現できる
• 次数を型で管理するのでフィールドに保持する必要がなく，ublas::bounded_vectorより常に4バイト分おトク
• 特に次元が2〜4のときのベクトルの場合の扱いやすさに
  • 添字じゃなくx,y,z,wというメソッド(参照返しなので代入も可)でアクセスできる
    • ただ例えば2次元のときzやwにアクセスするとコンパイルエラーじゃなく実行時エラーになるけどこれうまくやってなんとかできないだろうか
  • 基本的な型がtypedefされている(Vector2dとかMatrix3dとか)
  • コンストラクタに次元数だけの値を渡すことでベクトルの初期値を設定できる

こんな雰囲気．

#include <iostream>
#include <Eigen/Core>

#include <boost/format.hpp>

int main()
{
	std::cout << sizeof(Eigen::Matrix<int, 4, 1>) << std::endl;
	std::cout << sizeof(Eigen::Matrix<double, 4, 1>) << std::endl;
	std::cout << sizeof(Eigen::Matrix<double, 8, 1>) << std::endl;
	std::cout << sizeof(Eigen::Matrix<double, 8, 8>) << std::endl;

	//Eigen::Vector
	Eigen::Matrix3d m;
	std::cout << m << std::endl;

	Eigen::Vector2d v2d(0, 1);
	//Eigen::Vector3d v3d = v2d;	//ちゃんとコンパイルエラー．(コピコンじゃなく)代入でももちろん
	Eigen::Vector3d v3d(0, 1, 2);

	//演算も型が違うとエラー
	//std::cout << (v2d + v3d) << std::endl;
	//Eigen::Matrix<double, 3, 3>() * Eigen::Matrix<double, 4, 4>();

	//単に型が違うとエラーなんじゃなくて，算術的に有効な型同士なら問題ない
	Eigen::Matrix<double, 3, 1>() * Eigen::Matrix<double, 1, 3>();

	//Eigenではベクトルは1行または1列の行列
	Eigen::Matrix<double, 3, 1> vec = Eigen::Vector3d(0, 0, 0);	//なので代入できる
	

	//4次元まではコンストラクタで初期値を指定でき，x,y,z,wというメンバで参照(左辺に置いて代入することも)できる
	Eigen::Vector2f v2f(0, 1);
	std::cout << boost::format("{%1%,%2%}") % v2f.x() % v2f.y() << std::endl;

	Eigen::Vector3f v3f(0, 1, 2);
	std::cout << boost::format("{%1%,%2%,%3%}") % v3f.x() % v3f.y() % v3f.z() << std::endl;

	Eigen::Vector4f v4f(0, 1, 2, 3);
	std::cout << boost::format("{%1%,%2%,%3%,%4%}") % v4f.x() % v4f.y() % v4f.z() % v4f.w() << std::endl;
}

 16
 32
 64
 512
 -7.30135e-42 -1.27833e-41 -1.37662e-41
 -1.27833e-41 -7.30002e-42 3.64203e-314
 9.38725e-323 -6.62786e-42 -1.27992e-41
 {0,1}
 {0,1,2}
 {0,1,2,3}

まとめ

あれ，結局何が言いたかったんだろう．
あ，思い出した．

• uBLASには低次元で扱いやすいbounded_vectorがあるよ！
• でも次数周りでちょっとハマりやすいかもだから気をつけてね！
• Eigenも便利だよ！
• 便利だからって，次数を数百とかにして乱用することないようにね！

さて，遅刻にぐだぐだな内容と徹底的にハードルを下げてしまったので，お次の@yutoppさん，よろしくお願いします！！