Eigen.org

私の研究分野では、実対称行列の固有値と固有ベクトルを求めることが多い。 固有値問題は最も基本的な数値計算であり、様々なパッケージが簡便に利用できる。 ここでは、C++のライブラリであるEigenを使って実対称行列の固有値を求める方法を、備忘のためにまとめる。

インストール

$ cd /home/shiraishi/.local
$ curl -O "https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.tar.gz"
$ tar xf eigen-3.4.0.tar.gz
$ mv ./eigen-3.4.0/Eigen ./include

C++コード

$3 × 3$ の対称行列の固有値と固有ベクトルを表示するプログラム

#include <iomanip>
#include <iostream>
#include <vector>
#include <random>
#include <Eigen/Dense>

int main() {
    // 適当な対称行列をstd::vectorとして作る
    std::mt19937 mt(12345);
    std::normal_distribution<double> dist(0.0, 1.0);
    std::vector<double> mat(3*3, 0.0);
    mat[3*0 + 0] = dist(mt);
    mat[3*0 + 1] = dist(mt);
    mat[3*0 + 2] = dist(mt);
    mat[3*1 + 1] = dist(mt);
    mat[3*1 + 2] = dist(mt);
    mat[3*2 + 2] = dist(mt);

    mat[3*1 + 0] = mat[3*0 + 1];
    mat[3*2 + 0] = mat[3*0 + 2];
    mat[3*2 + 1] = mat[3*1 + 2];

    // std::vectorへのマップとしてEigen::MatrixXdを作る
    Eigen::Map<Eigen::MatrixXd> p(mat.data(), 3, 3);

    // ここで固有値を計算
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(p);

    // 固有値・固有ベクトルを取得
    Eigen::MatrixXd val = es.eigenvalues().real();
    Eigen::MatrixXd vec = es.eigenvectors().real();

    // 取得した固有値・固有ベクトルを出力する
    for (int i=0; i<3; i++) {
        double lambda = val.data()[i];
        std::vector<double> veci(3);
        std::copy(vec.data()+3*i, vec.data()+3*i + 3, veci.begin());

        std::cout << lambda << std::endl;
        std::cout << veci[0] << ","
                  << veci[1] << ","
                  << veci[2] << std::endl;
        std::cout << "----------" << std::endl;
    }
}

コメント：

• 他で扱うのが楽なので、行列は std::vector で計算して、mapしてEigenに渡している
• 固有値、固有ベクトルの値を取り出すのも、扱いが気楽な std::vector に一旦コピーしている

CMake

CMakeはこんな感じで書ける。

cmake_minimum_required(VERSION 3.13)

project(test_cmake CXX)

add_executable(a.out main.cpp)

target_include_directories(a.out PRIVATE /home/shiraishi/.local/include)
target_link_directories(a.out PRIVATE /home/shiraishi/.local/lib)

target_compile_options(a.out PRIVATE -O3 -inline-forceinline -xCORE_AVX512)
target_compile_features(a.out PRIVATE cxx_std_17)

MKLをバックエンドで使う場合は、上記プログラムを

#define EIGEN_USE_MKL_ALL
#include <Eigen/Dense>

と書き換えて、CMakeLists.txtを以下の通りにする。

cmake_minimum_required(VERSION 3.13)

project(test_cmake CXX)

add_executable(a.out main.cpp)

target_include_directories(a.out PRIVATE /home/shiraishi/.local/include ${MKLROOT}/include)
target_link_directories(a.out PRIVATE /home/shiraishi/.local/lib ${MKLROOT}/lib/intel64)

target_compile_options(a.out PRIVATE -O3 -inline-forceinline -xCORE_AVX512 -DMKL_LP64)
target_compile_features(a.out PRIVATE cxx_std_17)

target_link_libraries(a.out mkl_intel_lp64 mkl_sequential mkl_core pthread m dl)

実行

$ mkdir build
$ cd build
$ cmake -DCMAKE_CXX_COMPILER=icpc ../
$ make
$ ./a.out
-1.48416
0.304764,-0.571674,-0.761779
----------
0.34374
0.677115,0.692535,-0.248818
----------
1.28243
-0.669801,0.439982,-0.598149
----------

ベンチマーク

東京大学情報基盤センターのスパコンOakbridge-CXを使って、MKLのバックエンドへの使用の有無で対角化の時間を測った。 使用したコードはリポジトリを参照。

行列サイズ $N$ が大きくなるとMKLが速いのだが、 $N = 1000$ 程度ならMKLの方が遅くなってしまうようだ。

test.png