近年C++でよく利用されている線型代数ライブラリであるEigenの使い方について簡単にまとめました。私もたまに利用するのですが、使用頻度が低いため使い方を忘れてしまいます。忘れる度にネット検索して調べるのですが、いいかげん面倒になってきたので、自分用の備忘録として使い方を書き留めておくことにしました。
はじめに
線型代数ライブラリとはベクトルや行列に関する色々な数値計算を実行するための関数（サブルーチン）をまとめたものです。行列同士の掛け算などの基本的な計算であれば簡単に自分でプログラミングすることも可能かと思います。プログラミング初心者には多重ループを理解するための良い演習問題になるかもしれません。しかしながら実際の計算機シミュレーションにおいては、プログラミング初心者が作成するような行列同士の掛け算プログラムでは、実行速度が遅くて使い物になりません。したがって行列やベクトルの計算をする場合は自分でプログラミングせずに、線型代数ライブラリを利用することが普通です。
線型代数ライブラリで一番有名なのはBLAS/LAPACKです。これはどこのスーパーコンピューターにでもインストールされています。しかもそれぞれのシステムに最適化されているので、メチャメチャ計算が速いです。そんなBLAS/LAPACKでも現代的なプログラミング言語からすると古臭く、使い勝手が良くないです。（あくまでも私の感想ですが）もともとFORTRANという古い言語用に開発されたものなので仕方ないのでしょう。特にC++のような常に進化し続けている言語で利用すると、明らかに使い勝手が悪く感じます。というかFORTRAN用のものをC言語で利用しようとする時点でもう面倒くさいです。そのような問題を解決してくれるのがEigenというライブラリです。
Eigenはヘッダーファイルのみで構成されるテンプレートライブラリとよばれるもので、事前にコンパイルしてビルドする必要がありません。STLのような設計を選択している時点で神ライブラリです。もちろんEigenの計算速度は高速で、CPUのSIMD命令を使いベクトル化しくれますし、一部の計算ではスレッド並列化もしてくれます。私自身は利用したことは無いですが、コンパイルオプションによってBLAS/LAPACKを利用するように切り替える事も出来るようです。そんな神ライブラリも常に利用しているわけではないので、Eigenを利用する場合はいつもWEBで調べながらプログラミングしています。「あーそうだった！」とか独り言しながらプログラミングするのも飽きてきたので、自分用にまとめてみようと思いました。
インストール
Linuxにインストールする場合はパッケージ管理コマンドを利用すると簡単です。Debian系では
$ sudo apt install -y libeigen3-devとなり、RedHat系では
$ dnf -y install eigen3-develとなります。最新版を利用したい場合はEigenのWebページ(http://eigen.tuxfamily.org/)から、
アーカイブファイルをダウンロードして、適当なディレクトリで展開すれば良い。この記事を書いている時点の最新版はバージョン3.4.0で、Eigenのページからeigen-3.4.0.tar.gzというファイルがダウンロード出来ます。例えばこれを/usr/local/incluceの下で展開すると、/usr/local/include/eigen-3.4.0というディレクトリが作成されるので、/usr/local/incluce/eigenにリンクを作成する。
# ln -s /usr/local/include/eigen-3.4.0 /usr/local/incluce/eigen簡単な使い方
Eigenを利用するプログラムを記述するときは、<Eigen/Dense>をインクルードする。using namespace Eigen;と暗黙的に利用する名前空間を設定しておくと便利です。簡単なサンプルプログラムは次の通りです。
1 #include <Eigen/Dense>
2 #include <iostream>
3
4 using namespace std;
5 using namespace Eigen;
6
7 int main() {
8 Matrix3d m = Matrix3d::Random();
9 m = (m + Matrix3d::Constant(1.2)) * 50;
10 cout << "m =" << endl << m << endl;
11 Vector3d v(1, 2, 3);
12
13 cout << "m * v =" << endl << m * v << endl;
14 }このプログラムをsample1.cppとして保存して、次のようなコマンドでコンパイルする。
$ g++ -I/usr/include/eigen3 sample1.cppここで-IオプションでEigenがインストールされているディレクトリを指定する。aptやdnfでインストールした場合は/usr/include/eigen3となる。
ベクトルや行列の定義と初期化
3次元ベクトルの定義
Vector3d v; 実数
Vector3cd v; 複素数3次正方行列
Matrix3d A; 実数
Matrix3cd A; 複素数N次元ベクトルの定義
VectorXd v(N);
VectorXcd v(N);N次正方行列
MatrixXd A(N,N); 実数
MatrixXcd A(N,N); 複素数ベクトルのゼロ初期化
V=VectorXd::Zero(N);ベクトルの乱数で初期化
V=VectorXd::Random(N);行列のゼロ初期化（N行M列）
A=MatrixXd::Zero(N,M);行列の乱数で初期化（N行M列）
A=MatrixXd::Random(N,M);単位行列として初期化（N次正方行列）
A=MatrixXd::Identity(N,N);演算
ベクトルの基本演算
VectorXd a = VectorXd::Random(5);
VectorXd b = VectorXd::Random(5);
double c = 10.0;
VectorXd x = a + c * b;ベクトルの内積
VectorXd a = VectorXd::Random(5);
VectorXd b = VectorXd::Random(5);
double x = a.dot(b);行列とベクトルの基本演算
MatrixXd A = MatrixXd::Random(5, 5);
VectorXd b = VectorXd::Random(5);
VectorXd x = A * b; 行列とベクトルの積
VectorXd y = b.transpose() * A; ベクトルと行列の積
double c = b.transpose() * A * b; 二次形式の値行列同士の基本演算
MatrixXd A = MatrixXd::Random(5, 5);
MatrixXd B = MatrixXd::Random(5, 5);
MatrixXd C = A * B;操作
ベクトルと行列の転置
VectorXd a;
cout << a << endl << endl << a.transpose() << endl;
MatrixXd A;
cout << A << endl << endl << A.transpose() << endl;ベクトルからの切り出し
VectorXd a;
a.head(n); 先頭からn要素
a.tail(n); 末尾のn要素
a.segment(i,n); i+1番目からn要素を切り出す行列からの切り出し
MatrixXd A;
A.row(i); i+1番目の行
A.col(j); j+1番目の列
A.topLeftCorner(p,q); 左上のp行q列
A.bottomLeftCorner(p,q); 左下のp行q列
A.topRightCorner(p,q); 右上のp行q列
A.bottomRightCorner(p,q); 右下のp行q列
A.topRows(p); 上からp行
A.bottomRows(p); 下からp行
A.leftCols(q); 左からq列
A.rightCols(q); 右からq列
A.middleCols(i,q); i+1番目の列からq個の列を切り出し
A.middleRows(i,p); i+1番目の行からp個の行を切り出し
A.block(i,j,p,q); i+1番目の行、j+1番目の列からp行q列を切り出し平方根などの初等関数の計算
ベクトル(1,2,3,4,5)のそれぞれの要素の平方根の計算は、arratyメソッドを利用する。
VectorXd v(5);
v << 1, 2, 3, 4, 5;
VectorXd u=v.array().sqrt();行列についても同様にできる。
MatrixXd A(3, 3);
A << 1, 2, 3, 4, 5, 6, 7, 8, 9;
MatrixXd B = A.array().sqrt();使える関数は、abs(), log(), exp(), sin(), cos(), tan(), square(), cube()などが使える。
詳しくはEigenのドキュメントを参照してください。
総和や平均などの計算
ベクトルの要素の総和、平均、最大値、最小値の計算
VectorXd v = VectorXd::Random(10);
cout << v.sum() << endl;
cout << v.mean() << endl;
cout << v.maxCoeff() << endl;
cout << v.minCoeff() << endl;ベクトルのノルムの計算
VectorXd v = VectorXd::Random(10);
cout << v.norm() << endl; L2ノルム
cout << v.squaredNorm() << endl; L2ノルムの2乗
cout << v.lpNorm<1>() << endl; L1ノルム行列の列毎の総和、平均、最大値、最小値の計算
MatrixXd A = MatrixXd::Random(5, 5);
cout << A.colwise().sum() << endl;
cout << A.colwise().mean() << endl;
cout << A.colwise().maxCoeff() << endl;
cout << A.colwise().minCoeff() << endl;行列の行毎の計算はrowwise()を利用する。
自分が利用するであろうベクトルと行列の基本的な計算については大体記述できたかな。連立一次方程式や固有値問題などについては次回のブログでまとめたいと思います。
Category: プログラミング関連