福州网站建设设计公司百度浏览器入口

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))

结束