C++中的数学计算库Eigen
Eigen是C++中常用的多功能数学库，支持矩阵运算，固定大小矩阵/任意大小的密集矩阵/稀疏矩阵都能支持。同时支持各种数值类型，包括复数整数等，可扩展到自定义的数值类型。 支持常用矩阵运算，包括各种矩阵分解。Eigen库还有一些非官方模块生态系统，比如非线性优化、矩阵函数、多项式求解器、FFT（快速傅里叶变换）。
用起来也挺友好，Eigen库没有源文件，只有头文件，只需要在编译的时候-I对应头文件目录就行了。代码中直接include进来，比如：
#include <Eigen/Dense>
#include <Eigen/Core>
常用的数据类型
Eigen使用中会见到各种各样的矩阵数据类型，其实质都是
Matrix<Type, Size, Size>
各种诸如MatrixXd/Matrix2Xd/MatrixXi/Vector2d/Vector4f都是typedef所定义的别名。
Matrix/Vector/RowVector数据类型后面会有后缀。其中数字表示某个固定size，X表示Dynamic动态大小，d表示double，同理i表示int，在预编译时，宏会在typedef中替换掉定义。
- 比如
Matrix2Xd等价于Matrix<double, 2, Dynamic>； Vector表示列向量，本质是多行单列的矩阵，Vector4f就是4个float元素构成的向量；RowVector表示行向量，本质是单行多列的矩阵；
另外还有Array数据类型，不能当矩阵用，用来做些运算的广播挺方便。跟Matrix相比，各种运算的定义上会有区别，比如乘法，对于Array来讲是pointwise乘积,相当于numpy中的A.*B；对于Matrix，是正经的矩阵乘法，相当于numpy中的A@B。
例如，下面的点态乘法结果是[[0,0],[0,0]].
Eigen::Array<double, 2, 2> A(2, 2), B(2, 2);
A << 1, 0, 0, 1;
B << 0, 1, 1, 0;
std::cout << "A*B=\n" << A * B << std::endl;
而矩阵乘法的结果是[[0,1],[1,0]]
Eigen::Matrix<double, 2, 2> A(2, 2), B(2, 2);
A << 1, 0, 0, 1;
B << 0, 1, 1, 0;
std::cout << "A*B=\n" << A * B << std::endl;
Array和Matrix可以相互转换，用.array()与.matrix()就能实现，比如下面又变成Array的点态乘法了
Eigen::Matrix<double, 2, 2> A(2, 2), B(2, 2);
A << 1, 0, 0, 1;
B << 0, 1, 1, 0;
std::cout << "A*B=\n" << A.array() * B.array() << std::endl;
这样转化有时也比较麻烦，在矩阵中进行点乘的话也可以用cwiseProduct,还有很多的cwiseXXXX方法可以使用，比如abs/sqrt/sign等等。
更多参考：https://eigen.tuxfamily.org/dox/group__TutorialArrayClass.html
常用初始化与赋值
矩阵必须要初始化才能使用，不然会看见一些奇怪的值。下面就是一些常用的初始化
- 随机初始化（标准正态分布），
numpy中的np.randn(2,3)
Eigen::Matrix<double, 2, 3>::Random(2, 3);
mat.setRandom(); //或者这个 - 常数值初始化，numpy中的
np.ones(2,3)*2.0
Eigen::Matrix<double, 2, 3> A = Eigen::Matrix<double, 2, 3>::Constant(2.0); - 单位阵初始化，numpy中的
np.identity(3)
Eigen::Matrix<double, 3, 3> identityMatrix = Eigen::Matrix<double, 3, 3>::Identity();
或者下面这些更方便的：
a.setRandom()a.setOnes()a.setZero()a.setIdentity()
然后是赋值：
- 单个元素给矩阵赋值很简单，用
<<运算符可以实现，依行优先顺序填入即可，注意是列优先的，跟numpy行优先不一样。 - 特定某个元素设置的话，可以用
C(2, 2) = 0.，这里用的圆括号，跟Numpy的方括号有区别。 - 也有类似于numpy中
D[1:3, 2:4]=1.的块操作，等价于下面的赋值
D.block<2, 2>(1, 2) = Eigen::Matrix<double, 2, 2>::Constant(1.0); - 特殊的块，选取列
.col(j)或者行.row(i)，比如 numpy中E[:, 1]=E[:,2], 等价于
E.col(1)=E.col(2);
也可以是
E(Eigen::all, 1) = E(Eigen::all, 2); - 切片操作，跟block操作差不多，numpy中的
F[:, 1:4]=F[:, 0:3]等价于
E(Eigen::all, Eigen::seqN(1,3))=E(Eigen::all, Eigen::seq(0,2))
这里需要注意的是seqN表示从2开始，数2个元素，seq表示0开始，2结束的序列，跟python的通常用法不同，2是被包含在内的，类似matlab。
常用运算
- 获取矩阵的行与列，（numpy中的
mat.shape)：
int mat_height = mat.rows();
int mat_weight = mat.cols(); - 运算比较
| numpy | Eigen | |
|---|---|---|
| F范数 | np.linalg.norm(mat) | mat.norm() |
| 1范数 | np.linalg.norm(mat, 1) | mat.norm<1>() |
| 无穷范数 | np.linalg.norm(mat, 'inf') | mat.norm<> |
| 转置 | mat.T | mat.transpose() |
| 等差间断 | np.linspace(0,1,10) | VectorXd::LinSpaced(10,0,1) |
| 条件矩阵 | mat[mat>1]=1 | mat=(mat.array()>1).select(1,mat) |
resize()操作，numpy默认行优先，Eigen中是列优先。resize会执行析构函数，先将原对象析构，再在新空间上填充，使用太频繁对性能有一点影响。resize()并不会检查元素数量，当resize前后元素数量不一致时，会得到没有初始化的矩阵。conservativeResize()操作，动态改变矩阵大小，保留之前的元素，形状不匹配的地方都截断了，会出现未初始化的值。跟resize做个比较：
Eigen::MatrixXi m(4, 3);
m << 1, 2, 3,
4, 5, 6,
7, 8, 9,
10, 11, 12;
m.conservativeResize(6, 2);
会得到：
1 2
4 5
7 8
10 11
-842150451 -842150451
-842150451 -842150451
而m.resize(6,2)会对元素进行重排：
1 8
4 11
7 3
10 6
2 9
5 12
规约运算
mat.sum(),mat.mean()跟numpy没有区别。
numpy中的mat.sum(1),mat.mean(0) 分别差不多等价于
mat.colwise(1).sum(), 跟mat.rowwise(1).sum()
另外numpy中的min差不多等价于minCoeff()
性能优化
Map
Eigen::Map的使用场景，通常是有一个已经分配好的内存块，想用Eigen库来操作这块内存，但又不想复制数据时，map可以将data指向这片内存。
Map is Eigen is used to "convert" or interface data types like arrays, std::vector, and std::array into Eigen matrices without copying them.
比较接近的例子是torch中的view方法，B=A.view(4,3)并不会分配新的内存存数据，对B元素的修改也会反映到A上。比如下面的代码：
Eigen::MatrixXd a(1,12);
a<<0,1,2,3,4,5,6,7,8,9,10,11;
Eigen::Map<Eigen::Matrix<double, 4, 3>> c(&a(0), 4, 3);
Eigen::Map<Eigen::Matrix<double, 6, 2>> d(&a(0), 2, 6);
c(1, 1) = 0; //会改变 c，d矩阵中的值
Ref
Ref跟Map作用类似，通常用在函数参数的传递上。比如，我们想传递矩阵中子块的值作为参数并可能更改值时，如果直接使用引用数据类型，就会报错，a.col(0)是右值，无法传递其引用。
void doubleVector(Eigen::Vector3d& a) {a = a*2;
}
doubleVector(a.col(0));
使用Ref可以解决这一问题：
void doubleVector(Eigen::Ref<Eigen::Vector3d> a) {a = a*2;
}
doubleVector(a.col(0));
如果只是为了避免增加新的临时变量、不更改值的话，用const Ref<const VectorXd>& x写形参。
此部分参考：https://eigen.tuxfamily.org/dox/classEigen_1_1Ref.html
例子：求解二维Poisson方程
Eigen有诸多高级功能，但运用本文所介绍的内容，已经可以凑合着写一个二维Poisson方程求解器了。如从前一样，我们考虑一个
−Δu=32π2sin(4πx)sin(4πy)(x,y)∈(0,1)×(0,1)u|∂Ω=0
其真值为u=sin(4πx)sin(4πy). 采用重排序的Gauss-Seidel迭代求解, 直观的Python与C++代码对比：
左边python，右边C++
C++代码如下：
# define M_PI 3.14159265358979323846
void solvePoisson() {int N;int k = 4;N = 100;float h = 1. / N;float h2 = h * h;float tol_error = 1e-3;Eigen::RowVectorXd x(N+1);Eigen::VectorXd y(N+1);x.setLinSpaced(N + 1, 0, 1);y.setLinSpaced(N + 1, 0, 1);Eigen::MatrixXd xx = x.replicate(N+1, 1);Eigen::MatrixXd yy = y.replicate(1, N+1);Eigen::MatrixXd source = 2*k*k*M_PI*M_PI*((xx * k * M_PI).array().sin())*((yy * k * M_PI).array().sin());Eigen::MatrixXd u_truth = ((xx * k * M_PI).array().sin()) * ((yy * k * M_PI).array().sin());Eigen::MatrixXd u(N + 1, N + 1);Eigen::MatrixXd u_old(N + 1, N + 1);u.setZero();u_old.setOnes();while((u-u_old).norm()>tol_error){u_old = u;u(Eigen::seqN(1, N / 2, 2), Eigen::seqN(1, N / 2, 2)) = \0.25 * (h2 * source(Eigen::seqN(1, N / 2, 2), Eigen::seqN(1, N / 2, 2)) + \u(Eigen::seqN(0, N / 2, 2), Eigen::seqN(1, N / 2, 2)) + \u(Eigen::seqN(2, N / 2, 2), Eigen::seqN(1, N / 2, 2)) + \u(Eigen::seqN(1, N / 2, 2), Eigen::seqN(0, N / 2, 2)) + \u(Eigen::seqN(1, N / 2, 2), Eigen::seqN(2, N / 2, 2)));u(Eigen::seqN(2, (N - 1) / 2, 2), Eigen::seqN(1, N / 2, 2)) = \0.25 * (h2 * source(Eigen::seqN(2, (N - 1) / 2, 2), Eigen::seqN(1, N / 2, 2)) + \u(Eigen::seqN(1, (N - 1) / 2, 2), Eigen::seqN(1, N / 2, 2)) + \u(Eigen::seqN(3, (N - 1) / 2, 2), Eigen::seqN(1, N / 2, 2)) + \u(Eigen::seqN(2, (N - 1) / 2, 2), Eigen::seqN(0, N / 2, 2)) + \u(Eigen::seqN(2, (N - 1) / 2, 2), Eigen::seqN(2, N / 2, 2)));u(Eigen::seqN(1, N / 2, 2), Eigen::seqN(2, (N - 1) / 2,2)) = \0.25 * (h2 * source(Eigen::seqN(1, N / 2, 2), Eigen::seqN(2, (N - 1) / 2, 2)) + \u(Eigen::seqN(1, N / 2, 2), Eigen::seqN(1, (N - 1) / 2, 2)) + \u(Eigen::seqN(1, N / 2, 2), Eigen::seqN(3, (N - 1) / 2, 2)) + \u(Eigen::seqN(0, N / 2, 2), Eigen::seqN(2, (N - 1) / 2, 2)) + \u(Eigen::seqN(2, N / 2, 2), Eigen::seqN(2, (N - 1) / 2, 2)));u(Eigen::seqN(2, (N - 1) / 2, 2), Eigen::seqN(2, (N - 1) / 2, 2)) = \0.25 * (h2 * source(Eigen::seqN(2, (N - 1) / 2, 2), Eigen::seqN(2, (N - 1) / 2, 2)) + \u(Eigen::seqN(1, (N - 1) / 2, 2), Eigen::seqN(2, (N - 1) / 2, 2)) + \u(Eigen::seqN(3, (N - 1) / 2, 2), Eigen::seqN(2, (N - 1) / 2, 2)) + \u(Eigen::seqN(2, (N - 1) / 2, 2), Eigen::seqN(1, (N - 1) / 2, 2)) + \u(Eigen::seqN(2, (N - 1) / 2, 2), Eigen::seqN(3, (N - 1) / 2, 2)));}std::cout << "Rel error: " << (u - u_truth).norm() / u.norm() << std::endl;
}
作为比对的Python代码：
import numpy.matlib
import numpy as np
N = 100
h = 1/N
h2 = h*h
tol_error=1e-3
x = np.linspace(0,1,N+1).reshape(1,-1)
y = np.linspace(0,1,N+1).reshape(-1,1)
xx = np.matlib.repmat(x, N+1, 1)
yy = np.matlib.repmat(y, 1, N+1)
source = 32*np.pi**2*np.sin(4*np.pi*xx)*np.sin(4*np.pi*yy)
u_truth = np.sin(4*np.pi*xx)*np.sin(4*np.pi*yy)
u = np.zeros((N+1,N+1))
u_old = np.ones((N+1, N+1))
while np.linalg.norm(u-u_old)>tol_error:u_old[:]=uu[1:-1:2, 1:-1:2]=0.25*(h2*source[1:-1:2, 1:-1:2]+u[0:-2:2, 1:-1:2]+u[2::2, 1:-1:2]+u[1:-1:2, 0:-2:2]+u[1:-1:2, 2::2])u[2:-1:2, 1:-1:2]=0.25*(h2*source[2:-1:2, 1:-1:2]+u[1:-2:2, 1:-1:2]+u[3::2, 1:-1:2]+u[2:-1:2, 0:-2:2]+u[2:-1:2, 2::2])u[1:-1:2, 2:-1:2]=0.25*(h2*source[1:-1:2, 2:-1:2]+u[1:-1:2, 1:-2:2]+u[1:-1:2, 3::2]+u[0:-2:2, 2:-1:2]+u[2::2, 2:-1:2])u[2:-1:2, 2:-1:2]=0.25*(h2*source[2:-1:2, 2:-1:2]+u[3::2, 2:-1:2]+u[1:-2:2, 2:-1:2]+u[2:-1:2, 1:-2:2]+u[2:-1:2, 3::2])
print('||u-u_truth||F/||u_truth||F = ', np.linalg.norm(u-u_truth)/np.linalg.norm(u_truth))