Help us understand the problem. What is going on with this article?

eigen3を使ってみた備忘録

More than 1 year has passed since last update.

概要

任意精度対角化(密行列,疎行列)を追い求め,色々試してます.昨日eigenというライブラリーに出会い,今日インストールして使ってみたのでログとして記録.

任意精度対角化としてMpackは密行列のsymmetric, Hermiteクラスを対角化を実装していますが,疎行列はそもそもlapackで提供していないということもあり粗行列を効率的に解くことはできないようです.

色々調べていくとEigenというc++テンプレートライブラリーが密行列,粗行列の対角化を実行可能なようです.かつunsupportedのmpfrのライブラリーを使うとどうも任意精度で対角化できそうなので忘れないうちにeigenを使ったログを残す.

ログを残す主たる理由は私はc++ユーザーではなくすぐ忘れてしまいそうなので.

環境
ubuntu 14.04

install

メインページEigenのGet it からlatestをダウンロード,展開.

そもそもテンプレートライブラリーというものが何なのか理解していないが,Eigen getting startを読むとテンプレートライブラリーはmake installする必要はないようで,展開したディレクトリでプログラムを掛けばそれだけで動くようだ.

試しに展開したディレクトリ(EigenやREADME.mdがあるディレクトリー)でgetting startにあるサンプルプログラム

#include <iostream>
#include <Eigen/Dense>
using Eigen::MatrixXd;
int main()
{
  MatrixXd m(2,2);
  m(0,0) = 3;
  m(1,0) = 2.5;
  m(0,1) = -1;
  m(1,1) = m(1,0) + m(0,1);
  std::cout << m << std::endl;
}

をコンパイルしてみたが

$ g++ test.cpp -o a.out
test.cpp:3:23: fatal error: Eigen/Dense: No such file or directory
compilation terminated.

とコンパイルが通らなかったので

#include "Eigen/Core"

としなければならなかった(Cでも同様な現象があった気がするが忘れてしまった).

さらに任意精度計算をするために eigen の unsupported モジュールMPFRC++ Support module
のサンプル

#include <iostream>
#include <Eigen/MPRealSupport>
#include <Eigen/LU>
using namespace mpfr;
using namespace Eigen;
int main()
{
  // set precision to 256 bits (double has only 53 bits)
  mpreal::set_default_prec(256);
  // Declare matrix and vector types with multi-precision scalar type
  typedef Matrix<mpreal,Dynamic,Dynamic>  MatrixXmp;
  typedef Matrix<mpreal,Dynamic,1>        VectorXmp;
  MatrixXmp A = MatrixXmp::Random(100,100);
  VectorXmp b = VectorXmp::Random(100);
  // Solve Ax=b using LU
  VectorXmp x = A.lu().solve(b);
  std::cout << "relative error: " << (A*x - b).norm() / b.norm() << std::endl;
  return 0;
}

のinclude 下記のように修正

#include "unsupported/Eigen/MPRealSupport"
#include "Eigen/LU"

しコンパイルするも

$ g++ test.cpp -o a.out -lgmp -lmpfr
In file included from test.cpp:2:0:
unsupported/Eigen/MPRealSupport:15:22: fatal error: Eigen/Core: No such file or directory
compilation terminated.

してもエラーを出す.しょうが無いので, Eigen/ と unsupported/ を /usr/local/include/ 以下に sudo コピーすることで,サンプルプログラムを実行可能になった.

MPRealSupport

コンパイルする際にgmpとmpfrのライブラリーをリンクする必要があるのでMPRealSupport
のサンプルは下記のようにしてコンパイルする必要がある

g++ test.cpp -o a -lgmp -lmpfr

任意精度対角化の準備

行列
$$
A = \begin{pmatrix}
0 & 1 & 0\\
2 & 0 & 2\\
0 & 1 & 0
\end{pmatrix}
$$

の固有値は$\pm2 ,0$である. double 精度,例えば numpy (python)で解けば

>>> import numpy as np
>>> np.set_printoptions(precision=18)
array([ -2.000000000000000000e+00,  -1.732555071631836086e-16,
         1.999999999999999556e+00])

とな0固有値に誤差が伝播する.Mathematicaで機械精度の行列を対角化した場合も同様で

In[15]:= Eigenvalues@N@{{0,1,0},{2,0,2},{0,1,0}}//FullForm
Out[15]//FullForm= List[-2.`,1.9999999999999996`,-1.732555071631836`*^-16]

である. もちろんnMathematicaは任意の精度で($\le\infty$)で対角化可能である(が速度が気になる).
eigen3 と mpfrc++を使って対角化するプログラムは

#include <iostream>
#include <Eigen/Eigenvalues>
#include <unsupported/Eigen/MPRealSupport>
#include <Eigen/Dense>
using namespace Eigen;

using namespace mpfr;

int main()
{
  // set precision to 256 bits (double has only 53 bits)
  mpreal::set_default_prec(53);
  // Declare matrix and vector types with multi-precision scalar type
  typedef Matrix<mpreal,Dynamic,Dynamic>  MatrixXmp;
  MatrixXmp A(3,3); 

  A(0,0) = 0;
  A(0,1) = 1;
  A(0,2) = 0;

  A(1,0) = 2;
  A(1,1) = 0;
  A(1,2) = 2;

  A(2,0) = 0;
  A(2,1) = 1;
  A(2,2) = 0;

  std::cout << A << std::endl;

  EigenSolver<MatrixXmp> es(A);
  std::cout << es.eigenvalues() << std::endl;

  return 0;
}

みたいな感じで書けば良いようで,実行結果は

$ g++ eigenvalues-mpfrc-min.cpp -o a -lgmp -lmpfr 
$ ./a 
0 1 0
2 0 2
0 1 0
          (-2,0)
           (2,0)
(-3.28269e-77,0)

と高精度で評価されていることが分かる.説明にあるように

  // set precision to 256 bits (double has only 53 bits)
  mpreal::set_default_prec(53);

とすると実行結果は

$ ./a 
0 1 0
2 0 2
0 1 0
         (-2,0)
          (2,0)
(1.77435e-16,0)

となりdouble程度の精度しかない.

興味深いことに 0 固有値のlapack (numpy)の結果と, 機械精度(double)のMathematicaの結果は精度範囲で一致している.
私の記憶が正しければMathematicaはlapack方で実装されていると思うので,同じ結果なのであろう.

一方eigen3は表示桁数を増大させる方法が分からないが見れる範囲で異なる値を持っている.
これは採用しているアルゴリズムが違うことに起因するのであろうか?

ちなみに,mpmath (python)でデフォルト精度(=仮数部16桁)で対角化すると

>>> import mpmath
>>> mpmath.eig(mpmath.matrix([[0,1,0],[2,0,2],[0,1,0]]))
([mpf('-1.9999999999999993'), mpf('1.6358461123585456e-16'), mpf('2.0')],

この結果も異なるので,これはこれで興味深い.

to do

MPRealSupportは名前から実数のみのサポートなのだろうか?複素数もサポートされているのだろうか?
そもそもc++を使い始めて半日.使い方がよくわかっていない.

また速度比較

(追記)

eigen, mpack, mathematica で対角化した際の速度を比較する.
mpackは解ける行列のクラスに制限があるので次の対称行列

$$
M=
\begin{pmatrix}
1&1&0&\cdots &0\\
1&1&1&\cdots &0\\
0&1&1&\cdots &0\\
&\vdots&&\ddots&\\
0&0&0&\cdots &1\\
\end{pmatrix}
$$
を80桁の精度で密行列として対角化する.

diag-time.png

やはりMpack早い.mpackのプログラムを実行中cpuが400になるのでblasで並列化されているのだろう.
Eigenもblasが使えるようだがどうすれば良いのだろうか.

凡例の ditto (Anorldi 20%)のはmathematicaの対角化ルーチンであるEigensystemでAnorldi法(Lanczos法)を指定して固有値を大きい方を20%分求めたもの.

小さい方から求めた際の,比較はこれ
diag-time-arnoldi.png
Arnoldi法やLanczos法であれば$10^5\times10^5$の対角化が現実味を帯びてくる.

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした