密行列 | Fluid Dynamics Lab. Kyoto Univ.

Eigenは, まあ手っ取り早く行列計算をするパッケージです．とにかく簡単にプログラムを書くのが目的らしいですよ．でも，もう一つの特徴は，通常あり得ないほど計算が高速であることです．

なんかこのページ，中身が古いというか要点外しすぎで萎えるんだけど?

Matrix

#include <Eigen/Dense> ←じゅげーむ！
Eigen::MatrixXd myArray(5,3);
Eigen::MatrixXd myBrray(3,5);
Eigen::MatrixXi myCrray(5,5);

これで5x3と3x5の蜜行列ができると. 要素のアクセスは，配列名(行，列)ばい．ちなみにXdのdが倍精度, Xi のiが整数じゃき．

myArray(0,0)=2.5;//0から始まるのがお茶目
std::cout << "Rows= " << myArray.rows() << " Cols= " << myArray.cols() << std::endl;

でもこれ，めんどくさかろ？下の方がぱっと見がええよ？

myArray << 2.5, 3.1, 4.2,
             5.2, 3.3, 1.0,
             3.1, 5.2, 5.2,
             -1.2,4.2, std::sqrt(2.5),
             5.3, 2.2, 5.6;

画面に書くのもちょろいっちゃん：

std::cout << myArray << std::endl;
std::cout << myArray.row(3) << std::endl; //0から数えて3(数学では4列目というやつ)
std::cout << myArray.col(1) << std::endl;   //律儀に縦ベクトルで出力するので注意

もちろん途中でresize可能：

Eigen::MatrixXd myDrrays;
myDrrays=myArrray;
Eigen::MatrixXd myErrays;
myErrays.resize(100,100);

行列ば足したり掛けたりしてよか:

auto D=myArray*myBrray;

ばってん下のはできんと. 倍精度なのか整数なのか，わかりよらん. 掛け算というか代入も無理．

auto E=D*myCrray;

こういう時には, cast<typename T>()関数ば使うったい:

Eigen::MatrixXd myFrrays=myCrray.cast<double>();
auto G=D*myCrray.cast<double>();

みゃるほどみゃーーー

他にも, 転置行列を得る transpose() とか trace(), diagonal()とか，まあ想像できるわな．

対角行列には特別な奴がおるったい：

using DiagonalXd=Eigen::DiagonalMatrix<double,Eigen::Dynamic>
DiagonalXd KAKU;
KAKU.resize(12);//引数は1つ
//KAKU(2)=1.5;     //なぜかこれはできん.
KAKU.diagonal()(2)=1.5;
//std::cout << KAKU    //これも，できにゃあ
Eigen::MatrixXd pKAKU=KAKU;  std::cout << pKAKU;

不思議な奴やね

Vector

Matrixを扱うにはベクトルが必要だが, これは VectorXd とかである. まあ普通に計算できるが, 内積は

Eigen::VectorXd A(3), B(3); A << 0,1,2; B << 1,2,3;
double innter=A.dot(B);  //これが内積

3次元に特化した Vector3d では，外積もある：

Eigen::Vector3d A, B; A << 0,1,2; B << 1,2,3;
auto C=A.cross(B);  //これが内積

これがあれば, 行列の列を入れ替えるとかが簡単なのよ：

Eigen::MatrixXd VV(4,5);
		VV << 1,2,3,4,5,
		11,12,13,14,15,
		21,22,23,24,25,
		31,32,33,34,35;
Eigen::VectorXd XX(4);
XX << -1.2,-1.3,-1.4,-1.5;
VV.col(2)=XX;

お前がstd::vector<double>とか持っているならしょきかはかんたんじゃき(C++17)

std::vector<double> Boke {1.2, 1.3, -2.4};
Eigen::Vector3d Kasu(Boke.data()); // Bokeの成分でKasuを作成

row(), col()の実態

実はrow(), col()で得られるのはブロック要素への参照というものである. いいですか, 参照であるし, ブロック要素なんですよ．ベクトルではないのです. 例えば 1xN 行列のj列を取り出そうとして

double X=N.col(j); // これはできない.結果は行列なんで
double X=N.col(j).value();  //これはできる

なーんと．ああそれと，Eigenではめっちゃ遅延評価するんで，普通に書いたら式のままでビルドされてしまう．実行時に値が確定する．だからauto使う時は注意ね.

auto B=A.transpose(); // B変数は行列型ではなく数式型. 行列がよかなら auto B=A.transpose().value() と書いて遅延評価ばなしにせんとね

ブロック行列

行列の一部を取り出せる.

Eigen::Matrix<double,4,4> MM;
MM <<   1,2,3,4,
...
Eigen::Matrix<double,1,3> m1;
m1=MM.block<1,3>(0,0);//1x3行列, 0,0起点

Matrix Reloaded

じつはMatrixは線形演算しかサポートしていない．もっと色々できるのがArrayである. Matrix MをArrayにするには M.array() を，Array AをMatrixにするには A.matrix() を使う：

auto X=myArray.arrat(); //XはArray型
auto Y=X.matrix();          //YはMatrix型

線形演算

いよいよ本領発揮だ．

QR分解

m x n 行列 Aは, m x m 直交行列 Q と m x n 上三角行列に分解できる. ハウスホルダー変換で実現するのが A.colPivHouseholderQr() であり, ハウスホルダー変換結果を表す HouseholderQrクラスのオブジェクトの参照を戻す．というわけで, A x = b を解きたい場合には, A.colPivHouseholderQr().solve(b) でOK. 詳細を知りたければ戻ってきたオブジェクトに尋ねれば良い.

分解は他にも色々準備されている.

固有値固有ベクトル

もうマニュアル見た方が

行列式と逆行列

正方行列では, A.determinant() と A.inverse() が利用できると．ばってん#include <Eigen/LU> をわすれんよう気いつけるとよ．

計算速度

計算速度は実に謎めいている. 準備されている関数は, ヘッダーライブラリーであるにも関わらず，速い！なんだかよくわからないがコンパイル済みのBLASよりも速い! わけがわからん. ところが，自分で書くと遅い．おかしい．なんだろうこれ

MatrixXd G(N,M)=...定義...;
std::cout << G.colwise().norm()<< std::endl

はひどく高速だが

MatrixXd G(N,M)=...定義...;
for(int i=0;i<G.cols();i++) std::cout << G.col(i).norm()<< std::endl

はひどく低速である. なんでやねんこれ？だから

MatrixXd G(N,M)=...定義...;
auto NORM=G.colwise().norm();
std::cout << NORM.col(12) << std::endl;

が高速で便利なのよ. 要素アクセスで，何か最適化がかかるのだろうけど・・・

丸田先生の解釈では，テンプレートだし，オブジェクトの実態が決まるまでコンパイルできないし，コンパイルするときには行列のサイズやら全部決定しているから，クソほど最適化がかかるのではないか？

だってさー．そんなもんなんですかね．まあええわ，クソほど早いのはいいことだし