C++ - 行列ライブラリ Eigen 試用

Updated:


C++ 用の行列(線形代数)ライブラリである Eigen の環境を構築し、試用してみました。

0. 前提条件

  • Debian GNU/Linux 11.5 (64bit) での作業を想定。
  • GCC 12.1.0 (G++ 12.1.0) (C++17) でのコンパイルを想定。

1. Eigen のインストール

Eigen のページから最新のアーカイブファイルを取得し、適当な場所(/usr/local/include が良かろう)に展開するだけ。(当記事執筆時点では 3.4.0 が最新)

$ wget https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.tar.gz
$ sudo tar xvf eigen-3.4.0.tar.gz -C /usr/local/include

Debian の場合、 apt で以下のようにしてもインストールできるが、バージョンが 3.3.9 と少し古い。(当記事執筆時点)
よって、今回、当方は上記のアーカイブを取得して展開する方法を採用した。

$ sudo apt -y install libeigen3-dev

2. テスト用ソースコードの作成

  • 何をやっているのか等については、ソースコード内のコメントを参照。

File: eigen_test.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
/***********************************************************
  C++ 行列ライブラリ Eigen のテスト
  * 様々な定義、宣言、計算方法を試行
  * ベクトルの変数は英小文字、行列の変数は英大文字としている

    DATE        AUTHOR       VERSION
    2022.10.20  mk-mode.com  1.00 新規作成

  Copyright(C) 2022 mk-mode.com All Rights Reserved.
***********************************************************/
#include <cstdlib>   // for EXIT_XXXX
#include <iostream>
//#include <Eigen/Dense>      // 対応するヘッダを個別に include するのが面倒な場合
#include <Eigen/Core>         // 基本的な行列定義・行列演算
#include <Eigen/Eigenvalues>  // 固有値・固有ベクトル計算 etc. (外積)

int main(int argc, char* argv[]) {
  // 定義
  Eigen::Vector3d a;       // 3 次元ベクトル(double) a 用
  Eigen::VectorXd b(3);    // 3 次元ベクトル(double) b 用
  Eigen::Vector3d c;       // 3 次元ベクトル(double) c 用
  Eigen::Vector3d d;       // 3 次元ベクトル(double) d 用
  Eigen::Vector3d e;       // 3 次元ベクトル(double) e (単位ベクトル)用
  Eigen::Vector3d z;       // 3 次元ベクトル(double) z (全てゼロ)用
  Eigen::Matrix3d A;       // 3 x 3 行列(double) a 用
  Eigen::MatrixXd B(3, 3); // 3 x 3 行列(double) b 用
  Eigen::Matrix3d E;       // 3 x 3 行列(double) e (単位ベクトル)用
  Eigen::Matrix3d Z;       // 3 x 3 行列(double) z (全てゼロ)用

  try {
    // 初期化(英小文字の変数はベクトル、英大文字の変数は行列)
    a << 1.2, -2.3, 3.4;
    b(0) = -4.5;
    b(1) =  5.6;
    b(2) = -6.7;
    c << 8.0, 3.0, 4.0;
    d << 7.0, 5.0, 2.0;
    e = Eigen::Vector3d::Identity();
    z = Eigen::Vector3d::Zero();
    A << -1.2,  2.3, -3.4,
          4.5, -5.6,  6.7,
         -7.8,  8.9, -9.0;
    B(0, 0) =  9.0;
    B(0, 1) = -8.7;
    B(0, 2) =  7.6;
    B(1, 0) = -6.5;
    B(1, 1) =  5.4;
    B(1, 2) = -4.3;
    B(2, 0) =  3.2;
    B(2, 1) = -2.1;
    B(2, 2) =  1.0;
    E = Eigen::Matrix3d::Identity();
    Z = Eigen::Matrix3d::Zero();
    std::cout << "a =" << std::endl
              << a << std::endl;
    std::cout << "b =" << std::endl
              << b << std::endl;
    std::cout << "c =" << std::endl
              << c << std::endl;
    std::cout << "d =" << std::endl
              << d << std::endl;
    std::cout << "e =" << std::endl
              << e << std::endl;
    std::cout << "z =" << std::endl
              << z << std::endl;
    std::cout << "A =" << std::endl
              << A << std::endl;
    std::cout << "B =" << std::endl
              << B << std::endl;
    std::cout << "E =" << std::endl
              << E << std::endl;
    std::cout << "Z =" << std::endl
              << Z << std::endl;
    std::cout << "---" << std::endl;

    // 計算・出力 - ベクトル
    std::cout << "a^T  = " << a.transpose() << std::endl;    // 転置
    std::cout << "a    = " << std::endl
              << a.head(3) << std::endl;                     // 先頭から3要素
    std::cout << "     = " << std::endl
              << a.tail(3) << std::endl;                     // 末尾から3要素
    std::cout << "a(1) = " << a.segment(1, 1) << std::endl;  // 第2要素から1要素
    std::cout << "b(0) = " << b(0) << std::endl;             // 第1要素
    std::cout << "b(1) = " << b(1) << std::endl;             // 第2要素
    std::cout << "b(2) = " << b(2) << std::endl;             // 第3要素
    std::cout << "a.dot(b) =" << a.dot(b) << std::endl;      // 内積
    std::cout << "c.cross(d) =" << std::endl
              << c.cross(d) << std::endl;                    // 外積
    std::cout << "a.norm() = " << a.norm() << std::endl;     // ノルム
    std::cout << "sqrt(a) = " << std::endl
              << a.array().sqrt() << std::endl;              // √ベクトル
    std::cout << "sum(a) = " << std::endl
              << a.sum() << std::endl;                       // 総和
    std::cout << "mean(a) = " << std::endl
              << a.mean() << std::endl;                      // 平均値
    std::cout << "min(a) = " << std::endl
              << a.minCoeff() << std::endl;                  // 最小値
    std::cout << "max(a) = " << std::endl
              << a.maxCoeff() << std::endl;                  // 最大値
    std::cout << "a + b =" << std::endl
              << a + b << std::endl;                         // ベクトル + ベクトル
    std::cout << "a - b =" << std::endl
              << a - b << std::endl;                         // ベクトル - ベクトル
    std::cout << "a * 3 =" << std::endl
              << a * 3 << std::endl;                         // ベクトル * スカラ
    std::cout << "b / 2 =" << std::endl
              << b / 2 << std::endl;                         // ベクトル / スカラ
              // +, -, *, / は代入演算子 +=, -=, *=, /= も使用可
    std::cout << "---" << std::endl;

    // 計算・出力 - 行列
    std::cout << "A^T  = " << std::endl
              << A.transpose() << std::endl;        // 転置
    std::cout << "A(0, *) = " << std::endl
              << A.row(0) << std::endl;             // 1行目
    std::cout << "A(*, 1) = " << std::endl
              << A.col(1) << std::endl;             // 2列目
    std::cout << "A(0-1, *) = " << std::endl
              << A.topRows(2) << std::endl;         // 上から2行
    std::cout << "A(1-2, *) = " << std::endl
              << A.bottomRows(2) << std::endl;      // 下から2行
    std::cout << "A(*, 0-1) = " << std::endl
              << A.leftCols(2) << std::endl;        // 左から2列
    std::cout << "A(*, 1-2) = " << std::endl
              << A.rightCols(2) << std::endl;       // 右から2列
    std::cout << "A(1, *) = " << std::endl
              << A.middleRows(1, 1) << std::endl;   // 2行目から1行
    std::cout << "A(*, 1) = " << std::endl
              << A.middleCols(1, 1) << std::endl;   // 2列目から1列
    std::cout << "A(1, 1) = " << std::endl
              << A.block(2, 1, 1, 1) << std::endl;  // 3行目から1行、2列目から1列
    std::cout << "A.inverse = " << std::endl
              << A.inverse() << std::endl;          // 逆行列
    std::cout << "sqrt(A) = " << std::endl
              << A.array().sqrt() << std::endl;     // √行列
    std::cout << "sum(A) = " << std::endl
              << A.sum() << std::endl;              // 総和
    std::cout << "mean(A) = " << std::endl
              << A.mean() << std::endl;             // 平均値
    std::cout << "min(A) = " << std::endl
              << A.minCoeff() << std::endl;         // 最小値
    std::cout << "max(A) = " << std::endl
              << A.maxCoeff() << std::endl;         // 最大値
              // 行/列毎に総和、平均/最小/最大値を算出する場合は、
              // A.rowwise()..., A.colwise()... とする
    std::cout << "A + B =" << std::endl
              << A + B << std::endl;                // 行列 + 行列
    std::cout << "A - B =" << std::endl
              << A - B << std::endl;                // 行列 - 行列
    std::cout << "A * B =" << std::endl
              << A * B << std::endl;                // 行列 * 行列
    std::cout << "A * 3 =" << std::endl
              << A * 3 << std::endl;                // 行列 * スカラ
    std::cout << "B / 2 =" << std::endl
              << B / 2 << std::endl;                // 行列 / スカラ
              // +, -, *, / は代入演算子 +=, -=, *=, /= も使用可
    std::cout << "det(A) = " << A.determinant() << std::endl;  // 行列式

    // 計算・出力 - ベクトルと行列
    std::cout << "A * b =" << std::endl
              << A * b << std::endl;                // 行列 * ベクトル
    std::cout << "a^T * B =" << std::endl
              << a.transpose() * B << std::endl;    // ベクトル転置 * 行列 => ベクトル転置
    std::cout << "a * b^T =" << std::endl
              << a * b.transpose() << std::endl;    // ベクトル * ベクトル転置 => 行列
    std::cout << "A(*, 0) * B(1, *) =" << std::endl
              << A.col(0) * B.row(1) << std::endl;  // 行列1列目 * 行列2行目 => 行列
  } catch (...) {
      std::cerr << "EXCEPTION!" << std::endl;
      return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

3. ソースコードのコンパイル

$ g++ -std=c++17 -Wall -O2 --pedantic-errors \
  -I /usr/local/include/eigen-3.4.0 -o eigen_test eigen_test.cpp

4. 動作確認

$ ./eigen_test
a =
 1.2
-2.3
 3.4
b =
-4.5
 5.6
-6.7
c =
8
3
4
d =
7
5
2
e =
1
0
0
z =
0
0
0
A =
-1.2  2.3 -3.4
 4.5 -5.6  6.7
-7.8  8.9   -9
B =
   9 -8.7  7.6
-6.5  5.4 -4.3
 3.2 -2.1    1
E =
1 0 0
0 1 0
0 0 1
Z =
0 0 0
0 0 0
0 0 0
---
a^T  =  1.2 -2.3  3.4
a    =
 1.2
-2.3
 3.4
     =
 1.2
-2.3
 3.4
a(1) = -2.3
b(0) = -4.5
b(1) = 5.6
b(2) = -6.7
a.dot(b) =-41.06
c.cross(d) =
-14
 12
 19
a.norm() = 4.27668
sqrt(a) =
1.09545
   -nan
1.84391
sum(a) =
2.3
mean(a) =
0.766667
min(a) =
-2.3
max(a) =
3.4
a + b =
-3.3
 3.3
-3.3
a - b =
 5.7
-7.9
10.1
a * 3 =
 3.6
-6.9
10.2
b / 2 =
-2.25
  2.8
-3.35
---
A^T  =
-1.2  4.5 -7.8
 2.3 -5.6  8.9
-3.4  6.7   -9
A(0, *) =
-1.2  2.3 -3.4
A(*, 1) =
 2.3
-5.6
 8.9
A(0-1, *) =
-1.2  2.3 -3.4
 4.5 -5.6  6.7
A(1-2, *) =
 4.5 -5.6  6.7
-7.8  8.9   -9
A(*, 0-1) =
-1.2  2.3
 4.5 -5.6
-7.8  8.9
A(*, 1-2) =
 2.3 -3.4
-5.6  6.7
 8.9   -9
A(1, *) =
 4.5 -5.6  6.7
A(*, 1) =
 2.3
-5.6
 8.9
A(1, 1) =
8.9
A.inverse =
 2.5427 2.63361       1
3.23967 4.33058       2
      1       2       1
sqrt(A) =
   -nan 1.51658    -nan
2.12132    -nan 2.58844
   -nan 2.98329    -nan
sum(A) =
-4.6
mean(A) =
-0.511111
min(A) =
-9
max(A) =
8.9
A + B =
 7.8 -6.4  4.2
  -2 -0.2  2.4
-4.6  6.8   -8
A - B =
-10.2    11   -11
   11   -11    11
  -11    11   -10
A * B =
 -36.63      30  -22.41
  98.34  -83.46   64.98
-156.85  134.82 -106.55
A * 3 =
 -3.6   6.9 -10.2
 13.5 -16.8  20.1
-23.4  26.7   -27
B / 2 =
  4.5 -4.35   3.8
-3.25   2.7 -2.15
  1.6 -1.05   0.5
det(A) = -3.63
A * b =
 41.06
 -96.5
145.24
a^T * B =
36.63   -30 22.41
a * b^T =
  -5.4   6.72  -8.04
 10.35 -12.88  15.41
 -15.3  19.04 -22.78
A(*, 0) * B(1, *) =
   7.8  -6.48   5.16
-29.25   24.3 -19.35
  50.7 -42.12  33.54

5. 参照


これで、 C++ での行列計算が容易に行えます。
(自分で配列や vector 等を使用して行列計算するのも勉強になって良いですが)

以上。





 

Sponsored Link

 

Comments