初出 https://qiita.com/kumpeishiraishi/items/f168f073b8433713935c

Eigenで固有値を計算するとstd::complexが返ってくる.md

やりたいこと

Eigenでdouble型の行列の固有値を計算し、結果を倍精度で標準出力する。

まずはPythonで正解を計算する

import numpy as np
A = np.array([[4,1],[2,-1]])
print(np.linalg.eig(A)[0][0])
print(np.linalg.eig(A)[0][1])

として、numpyで固有値を表示すると：

$ python3 hoge.py
4.372281323269014
-1.3722813232690143

となる。

コード例

# include <iostream>
# include <Eigen/Dense>
# include <complex>
# include <iomanip> //setprecision
# include <limits> //numeric_limits

using namespace Eigen;

int main() {
  Matrix<double, 2, 2> A;
  A << 4,1,
       2,-1;

  EigenSolver<Matrix<double,2,2> > es(A);
  
  std::cout << A << std::endl;
  std::cout << "----------------" << std::endl;
  std::cout << std::setprecision(std::numeric_limits<double>::digits10) << std::scientific << std::real(es.eigenvalues()[0]) << std::endl;
  std::cout << std::setprecision(std::numeric_limits<double>::digits10) << std::scientific << std::real(es.eigenvalues()[1]) << std::endl;
}

結果

$ g++ -I /EIGEN_DIRECTORY_FULL_PATH/ fuga.cpp
$ ./a.out
 4  1
 2 -1
----------------
4.372281323269014e+00
-1.372281323269014e+00

で正しい結果が返ってきた。

解説

公式ドキュメントなどで書かれている固有値の計算・出力方法は

EigenSolver<Matrix<double,2,2> > es(A);
std::cout << es.eigenvalues() << std::endl;

である。 このコードにより、

 (4.37228,0)
(-1.37228,0)

のように、行列Aの全ての固有値を表示させることができる。

上記のコード例にあるように、es.eigenvalues()[0]とすると、任意のindexの固有値のみを出力することができる。

行列の固有値は一般に複素数であるので、このとき出力されるのはstd::complex型である。 上の行列の場合は固有値が実数で、これでは具合が悪いので、std::real()とする。 もちろん、

std::cout << std::setprecision(std::numeric_limits<double>::digits10) << std::scientific << es.eigenvalues()[0] << std::endl;

とすれば、(4.372281323269014e+00,0.000000000000000e+00)のように実部・虚部の双方が、倍精度で出力できる。