C++の行列ライブラリ Eigenの紹介

Tag

# ライブラリ

Daisuke Okanohara

Executive Vice President

C++で行列計算をする場合に便利なライブラリEigenを紹介したいと思います。

ベクトル・行列演算は知っているからEigenの使い方だけを教えてくれというかたは最初の章は読み飛ばしてください。

—

多くの統計処理がベクトル・行列演算を用いるとコンパクトに表すことが知られています。ちょっと複雑そうにみえる問題も整理してみるとベクトル・行列演算で書ける場合が多いです。（ベクトル・行列という言葉に抵抗がある方はそれぞれを単に配列、配列の配列とでも思ってもらえればいいでしょう）。ベクトルの内積は\(u^T v = u_1 v_1 + u_2 v_2 + \ldots +\)として求められ、ベクトルのノルムは自分自身のベクトルとの内積の平方根、\(|u| = \sqrt{ u^T u} \)として求められます（以降ベクトルは全て列ベクトルを指すとします）。

例えば、あるユーザーの商品の購買履歴は、\(i\)番目の商品を買っていたら\(u_i=1\)、そうでなければ\(u_i=0\)であるようなベクトル\(u\)で表せます。二人のユーザー\(u, v\)の購買履歴が似ているかどうかは共通の商品を多く買っているかどうかで判断できるので、二つのベクトルの間の角度がどれだけ小さいか、つまりコサイン距離\(cos(u, v) = u^Tv/|u||v|\)で測ることが多いです。

また機械学習の多くも行列演算で表すことができます。入力\(x \in R^m\)から出力\(y \in \{1,2,\ldots, k\}\)を推定する多クラス分類器\(f(x)\)をパーセプトロンを使って学習する場合を考えてみます。モデルパラメータとしてm行k列の重み行列\(A\)を用意します。そしてi行目の値のみ1で他が全て0であるようなベクトルをe_iとしたとき、入力\(x\)に対する推定結果は\(f(x) = \arg \max_y x^TAe_y\)と表せます。また訓練データ\((x, y)\)が与えられた時の学習則は\(A := A + x(e_y – e_{f(x)})^T\)となります。

—

Eigenは行列・ベクトル演算をC++から簡単に使えるようにしたライブラリでありLGPL version3またはGPL version2のdual licenseで公開されています。最初、二人で始められたプロジェクトですが現在では数十人の開発者によって進められており、年内に大幅な性能向上がされたEigen3.0の正式版がリリースされる予定になっています。

Eigenの特徴について、公式ホームページから文言を借りて説明します。

多用途
固定、可変サイズの行列、密行列、疎行列のサポート、列/行方向のデータ格納、線形方程式solver、幾何、クオータニオンのサポート、複素数など様々な成分値のサポート
高速
Expression templatesを利用することにより、遅延評価を行い全体最適化を行うと同時に中間データ生成を極力押さえる。SSEなどベクトル演算が利用可能な場合は積極的に利用。不要な動的メモリ確保の除去、コスト評価した上でのループアンローリング、キャッシュ向けに最適化。
洗練されたAPI
擬似コードをそのままコードとして書けるよう直感的なAPI
多くのコンパイラをサポート
GCC , MSVC, ICC, MinGW, Sun Compilerなどで利用可能

Eigenの性能は有名な行列ライブラリBLASやその亜種に匹敵するほど優秀です（ベンチマーク）。また、わざわざ行列で表現する必要がないように思える、簡単なベクトルの足し算、内積などの演算でもCで直に書くより、Eigenを使ったほうがSSEなどのベクトル化をしてくれるので数倍高速になる場合が多く、利用する価値があります。

Eigenを実際に利用している分野はロボット、CG、ゲームなど多岐にわたり、Googleも機械学習やコンピュータビジョンで利用しいるそうです（公式ページより）

ではさっそく使ってみましょう。Eigenはテンプレートライブラリ、つまりヘッダファイルだけからなるので、ダウンロードしてきたヘッダファイルをソースからインクルードするだけで利用することができます。以降の例はEigen 3.0-beta2を利用した場合を想定しています。

まず、Eigenを利用するにはEigen/Coreを読み込んでおけば大抵の場合OKです。様々な種類の行列が利用できますが、ここでは動的サイズでfloat成分値である行列MatrixXfを使った例を示してみます。

以下のように行列の演算は大体そのまま思った通りにかけます（サイズが合わない演算の場合はコンパイルエラーになります。）

#inlcude <Eigen/Core>
using namespace Eigen;
....
MatrixXf X(100, 200);             // 100行200列の行列
MatrixXf Y(200, 50);              // 200行50列の行列
X(11,22) = 333;                   // 11行22列目に333を代入
...色々処理
MatrixXf A = X * Y;               // X*Yの結果をZに代入

MatrixXf Z = MatrixXf::Zero(100, 50); // 零行列で初期化
...色々処理
MatrixXf B = 1.5 * X * Y + 3.0 * Z; 
// X*Yなどの中間要素をつくらず直接Bに計算結果が代入

MatrixXfを使う上でちょっと注意しなければならないのは、std:vectorとは違って成分値が0初期化されないということです。0初期化したい場合は例にあるようにMatrixXf::Zero(rows, cols)を利用して初期化します。

また、行列の様々な部分行列を取り出す方法も充実しています。行、列を指定するrow(i), col(j)はもちろん、ベクトルを対角行列に変換するasDiagonal()、成分ごとの操作を行うcwise()などもあり、ほしいと思った機能は大抵揃っています。もちろんこれらの操作を組み合わせた後の全体最適化はEigenが勝手にやってくれます。

疎行列、つまり多くの成分が0であるような行列を扱う場合にはEigen/Sparseをインクルードし、SMatrixXfを利用します。疎行列は購買履歴や自然言語情報、画像解析（Bag of visual words)など、いろいろな場面ででてきます。

SMatrixXfでは、成分値が0でない場所と、その成分値のペアを格納するため、0でない成分数に比例した分だけしかメモリを使用せず、また演算も速いという特徴があります。このため数百万行、数百万列といった大規模行列も疎行列であれば高速に扱えます。しかし密行列と違って各成分に自由にアクセスできませんので通常のMatrixと比べて機能は制限されています。また、成分値の設定の仕方が前から順番に設定する方法のみが効率的にサポートされてます（これは疎行列を実装したことある人なら大体わかるとは思いますが）

#include <Eigen/Sparse>
using namespace Eigen;
...
SMatrixXf S(100, 200);
A.reserve(200);                 // 非零成分数を大体指定．
　　　　　　　　　　　　　 　　　　// サイズが少ない場合は適当に確保してくれる
for (int i = 0; i < 100; ++i){
  A.startVec(i);                 // i行目の成分をこれから設定する宣言
  A.insertBack(i, 10) = 5;       // i行10列目の成分値を5に設定
  A.insertBack(i, i*2) = 7;      // i行i*2列目の成分に7を設定
}
A.finalize();
MatrixXf X (300, 100);
MatrixXf Y = X * A;              // 行列演算は大体普通にできる

より詳しい使い方は
チュートリアル（英語）が充実していますのでそちらを参考にしてみてください。

また、私が先日公開した大規模行列向け特異値分解ライブラリredsvdもEigenを利用していますので参考にしてみてください。redsvd src

Eigenはまだ発展途上です。最善の手法をまだ使ってない場合も多いですし（すごい勢いで次々実装されていますが）、APIも安定していないです。が、これまで紹介してきたメリットは非常に大きく、試してみる価値はあると思います。是非使ってみてください。

Tag

# ライブラリ