C++矩阵库：Eigen 3.4.90 中文使用文档（一）

标签：Eigen 矩阵 C++ matrix 3.4 数组 include mat

写在前面：我在学习Eigen库时，没找到好的中文文档，因此萌发了汉化Eigen官网文档的想法。其中一些翻译可能不是特别准确，欢迎批评指正。感兴趣的同学可以跳转到官网查看原文：

Eigen: Main Pagehttps://eigen.tuxfamily.org/dox/index.html

Eigen库，是一个开源的C++模板库，专门用于线性代数计算。本文只讲运用，不讨论安装问题。如果还未安装Eigen，请先自行安装。共分为四大部分，这篇文章属于第一份部分，主要介绍：Eigen中常用的矩阵和数组操作。

一、模版类：Matrix

1、Matrix 模板类的前三个模板参数

Matrix 类总共有六个模板参数，但我们只需先了解前三个参数，剩下的三个参数有默认值，可以暂时忽略。

Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>

Scalar: 标量类型，即系数的类型。例如，如果你需要一个浮点数矩阵，可以选择 float 类型。
RowsAtCompileTime 和 ColsAtCompileTime: 在编译时已知的矩阵行数和列数。

2、便捷类型定义

Eigen提供了许多便捷的类型定义来简化常用矩阵和向量的定义。例如 Matrix4f 表示 4x4 的浮点数矩阵，定义如下：

typedef Matrix<float, 4, 4> Matrix4f;

3、向量

在Eigen中，向量其实可以看做是具有一行或者一列的矩阵。其中列向量是最常见的。例如 Vector3f 其实是一个利用 typedef 定义的3维浮点数列向量，定义如下：

typedef Matrix<float, 3, 1> Vector3f;

Eigen也提供了行向量的便捷类型定义，其本质上也是一个矩阵。例如一个2维的整形行向量其实本质上也是一个矩阵。定义如下：

typedef Matrix<int, 1, 2> RowVector2i;

4、动态尺寸

Eigen不仅限于在编译时创建已知尺寸的矩阵和向量。模板参数 RowsAtCompileTime 和 ColsAtCompileTime 可以取特殊值 Dynamic，表示尺寸在运行时确定。例如，MatrixXd 表示动态尺寸的双精度矩阵，定义如下：

typedef Matrix<double, Dynamic, Dynamic> MatrixXd;

5、构造函数

默认构造函数不进行动态内存分配，也不初始化矩阵系数。例如：

Matrix3f a;   // 一个 3x3 的未初始化浮点数矩阵
MatrixXf b;   // 一个动态尺寸的未初始化浮点数矩阵，当前尺寸为 0x0

还可以使用带尺寸的构造函数来分配系数数组，但不初始化系数：

MatrixXf a(10, 15);  // 一个 10x15 的未初始化浮点数矩阵
VectorXf b(30);      // 一个尺寸为 30 的未初始化浮点数向量

可以使用初始化列表来初始化矩阵和向量的系数：

Vector2d a(5.0, 6.0);   // 列向量
Vector3d b(5.0, 6.0, 7.0); // 列向量

MatrixXi m {
  {1, 2}, 
  {3, 4}
};

6、系数访问器

Eigen 提供了重载的括号操作符来访问和修改系数。对于矩阵，先行后列，如果你安装了Eigen你可以执行下面的代码。例如：

Eigen::MatrixXd m(2, 2);
m(0, 0) = 3;
m(1, 0) = 2.5;
m(0, 1) = -1;
m(1, 1) = m(1, 0) + m(0, 1);
std::cout << "Here is the matrix m:\n" << m << std::endl;

7、大小调整

可以使用 resize() 方法调整动态尺寸矩阵的大小。例如：

Eigen::MatrixXd m(2, 5);
m.resize(4, 3);
std::cout << "The matrix m is of size " << m.rows() << "x" << m.cols() << std::endl;
std::cout << "It has " << m.size() << " coefficients" << std::endl;

Eigen::VectorXd v(2);
v.resize(5);
std::cout << "The vector v is of size " << v.size() << std::endl;

但是对于固定大小的矩阵和向量，resize() 方法的操作是无效的。

Eigen::Matrix4d m;
m.resize(4, 4); // 无操作
std::cout << "The matrix m is of size " << m.rows() << "x" << m.cols() << std::endl;

8、固定尺寸与动态尺寸

对于使用固定尺寸还是动态尺寸，官方给的建议是：在尺寸非常小（如小于 16）且已知的情况下，使用固定尺寸有助于提高性能。对于较大的尺寸，或在编译时尺寸不确定的情况下，使用动态尺寸更为合适。例如：

Matrix4f mymatrix;  // 固定尺寸矩阵
MatrixXf mymatrix(rows, columns);  // 动态尺寸矩阵

9、可选模版参数

其实Matrix 类一共有六个模版参数，前三个参数是必选的，也是最常用，后面三个模版参数是可选的。Matrix 类的完整模板参数列表如下：

Matrix<typename Scalar,
       int RowsAtCompileTime,
       int ColsAtCompileTime,
       int Options = 0,
       int MaxRowsAtCompileTime = RowsAtCompileTime,
       int MaxColsAtCompileTime = ColsAtCompileTime>

Options 是一个位字段，例如 RowMajor 表示行优先存储。默认情况下，矩阵采用列优先存储。例如，以下定义表示行优先的 3x3 矩阵：

Matrix<float, 3, 3, RowMajor>

MaxRowsAtCompileTime 和 MaxColsAtCompileTime 用于指定矩阵的最大尺寸，可以避免动态内存分配。例如：

Matrix<float, Dynamic, Dynamic, 0, 3, 4>

10、便捷类型定义

Eigen 定义了以下便捷类型：

MatrixNt 表示 Matrix<type, N, N>，例如 MatrixXi 表示 Matrix<int, Dynamic, Dynamic>。
VectorNt 表示 Matrix<type, N, 1>，例如 Vector2f 表示 Matrix<float, 2, 1>。
RowVectorNt 表示 Matrix<type, 1, N>，例如 RowVector3d 表示 Matrix<double, 1, 3>。

这些类型定义使得使用 Eigen 更加方便和直观。

二、使用Eigen进行矩阵和向量运算

Eigen提供了一个全面的运算库，用于矩阵和向量的算术操作，使得矩阵和向量的操作和计算变得直观和简单。以下是基本操作的概述，包括加法、减法、标量乘法、转置、矩阵乘法和基本的归约操作。

1、加法和减法

矩阵和向量的加法与减法要求左右操作数具有相同的行数和列数，并且必须具有相同的标量类型，因为Eigen不会自动进行类型提升。可用的运算符包括：

二元运算符 +（如 a + b）
二元运算符 -（如 a - b）
一元运算符 -（如 -a）
复合运算符 +=（如 a += b）
复合运算符 -=（如 a -= b）

#include <iostream>
#include <Eigen/Dense>

int main()
{
  Eigen::Matrix2d a;
  a << 1, 2,
       3, 4;
  Eigen::MatrixXd b(2, 2);
  b << 2, 3,
       1, 4;
  std::cout << "a + b =\n" << a + b << std::endl;
  std::cout << "a - b =\n" << a - b << std::endl;
  std::cout << "执行 a += b;" << std::endl;
  a += b;
  std::cout << "现在 a =\n" << a << std::endl;
  Eigen::Vector3d v(1, 2, 3);
  Eigen::Vector3d w(1, 0, 0);
  std::cout << "-v + w - v =\n" << -v + w - v << std::endl;
}

输出：

a + b =
3 5
4 8
a - b =
-1 -1
 2  0
执行 a += b;
现在 a =
3 5
4 8
-v + w - v =
-1
-4
-6

2、标量乘法和除法

矩阵和向量的标量乘法和除法也非常简单。可用的运算符包括：

二元运算符 *（如 matrix * scalar 或 scalar * matrix）
二元运算符 /（如 matrix / scalar）
复合运算符 *=（如 matrix *= scalar）
复合运算符 /=（如 matrix /= scalar）

示例：

#include <iostream>
#include <Eigen/Dense>

int main()
{
  Eigen::Matrix2d a;
  a << 1, 2,
       3, 4;
  Eigen::Vector3d v(1, 2, 3);
  std::cout << "a * 2.5 =\n" << a * 2.5 << std::endl;
  std::cout << "0.1 * v =\n" << 0.1 * v << std::endl;
  std::cout << "执行 v *= 2;" << std::endl;
  v *= 2;
  std::cout << "现在 v =\n" << v << std::endl;
}

输出：

a * 2.5 =
2.5   5
7.5  10
0.1 * v =
0.1
0.2
0.3
执行 v *= 2;
现在 v =
2
4
6

3、转置和共轭

矩阵或向量的转置（a.transpose()）、共轭（a.conjugate()）和伴随矩阵（即共轭转置，a.adjoint()）通过成员函数分别获取。

示例：

#include <iostream>
#include <Eigen/Dense>

int main()
{
  Eigen::MatrixXcf a = Eigen::MatrixXcf::Random(2, 2);
  std::cout << "这是矩阵 a\n" << a << std::endl;

  std::cout << "这是矩阵 a^T\n" << a.transpose() << std::endl;

  std::cout << "这是 a 的共轭\n" << a.conjugate() << std::endl;

  std::cout << "这是矩阵 a^*\n" << a.adjoint() << std::endl;
}

输出：

这是矩阵 a
(-0.211,0.68) (-0.605,0.823)
 (0.597,0.566) (0.536,-0.33)
这是矩阵 a^T
(-0.211,0.68)  (0.597,0.566)
(-0.605,0.823) (0.536,-0.33)
这是 a 的共轭
(-0.211,-0.68) (-0.605,-0.823)
 (0.597,-0.566) (0.536,0.33)
这是矩阵 a^*
(-0.211,-0.68) (0.597,-0.566)
(-0.605,-0.823) (0.536,0.33)

对于实数矩阵，共轭操作是无效的，因此 adjoint() 等同于 transpose()。

需要注意的是，transpose() 和 adjoint() 返回一个代理对象而不进行实际转置。如果执行 a = a.transpose();Eigen会在转置操作完成之前开始写入 a，导致意外结果。解决方法是使用 transposeInPlace() 进行就地转置。

#include <iostream>
#include <Eigen/Dense>

int main()
{
  Eigen::MatrixXf a(2, 3);
  a << 1, 2, 3,
       4, 5, 6;
  std::cout << "这是初始矩阵 a:\n" << a << std::endl;

  a.transposeInPlace();
  std::cout << "转置后的矩阵 a:\n" << a << std::endl;
}

输出：

这是初始矩阵 a:
1 2 3
4 5 6
转置后的矩阵 a:
1 4
2 5
3 6

4、矩阵乘法和矩阵向量乘法

矩阵与矩阵的乘法、矩阵与向量的乘法通过运算符 * 实现。向量实际上是矩阵的一种特殊形式，因此矩阵向量乘法是矩阵乘法的一个特例。

示例：

#include <iostream>
#include <Eigen/Dense>

int main()
{
  Eigen::Matrix2d mat;
  mat << 1, 2,
         3, 4;
  Eigen::Vector2d u(-1, 1), v(2, 0);
  std::cout << "这是 mat*mat:\n" << mat * mat << std::endl;
  std::cout << "这是 mat*u:\n" << mat * u << std::endl;
  std::cout << "这是 u^T*mat:\n" << u.transpose() * mat << std::endl;
  std::cout << "这是 u^T*v:\n" << u.transpose() * v << std::endl;
  std::cout << "这是 u*v^T:\n" << u * v.transpose() << std::endl;
  std::cout << "让我们将 mat 自乘" << std::endl;
  mat = mat * mat;
  std::cout << "现在 mat 是:\n" << mat << std::endl;
}

输出：

这是 mat*mat:
 7 10
15 22
这是 mat*u:
1
1
这是 u^T*mat:
2 2
这是 u^T*v:
-2
这是 u*v^T:
-2 0
 2 0
让我们将 mat 自乘
现在 mat 是:
 7 10
15 22

5、点积和叉积

对于点积和叉积，需要使用 dot() 和 cross() 方法。当然，点积也可以通过 u.adjoint() * v 得到。

示例：

#include <iostream>
#include <Eigen/Dense>

int main()
{
  Eigen::Vector3d v(1, 2, 3);
  Eigen::Vector3d w(0, 1, 2);

  std::cout << "点积: " << v.dot(w) << std::endl;
  double dp = v.adjoint() * w; // 自动将内积转换为标量
  std::cout << "通过矩阵乘积获得的点积: " << dp << std::endl;
  std::cout << "叉积:\n" << v.cross(w) << std::endl;
}

输出：

点积: 8
通过矩阵乘积获得的点积: 8
叉积:
 1
-2
 1

注意，叉积仅适用于大小为3的向量。点积适用于任意大小的向量。对于复数，Eigen的点积在第一个变量上是共轭线性的，在第二个变量上是线性的。

6、基本算数归约操作

Eigen还提供了一些归约操作，将给定的矩阵或向量归约为单个值，如和（通过 sum() 计算）、积（prod()）、或所有系数的最大值（maxCoeff()）和最小值（minCoeff()）示例：

#include <iostream>
#include <Eigen/Dense>

using namespace std;

int main()
{
  Eigen::Matrix2d mat;
  mat << 1, 2,
         3, 4;
  cout << "这是 mat.sum(): " << mat.sum() << endl;
  cout << "这是 mat.prod(): " << mat.prod() << endl;
  cout << "这是 mat.mean(): " << mat.mean() << endl;
  cout << "这是 mat.minCoeff(): " << mat.minCoeff() << endl;
  cout << "这是 mat.maxCoeff(): " << mat.maxCoeff() << endl;
  cout << "这是 mat.trace(): " << mat.trace() << endl;
}

输出：

这是 mat.sum(): 10
这是 mat.prod(): 24
这是 mat.mean(): 2.5
这是 mat.minCoeff(): 1
这是 mat.maxCoeff(): 4
这是 mat.trace(): 5

矩阵的迹（trace）由函数 trace() 返回，是对角线系数的和，也可以使用 a.diagonal().sum() 高效计算。

此外，还有返回最小和最大系数坐标的 minCoeff 和 maxCoeff 函数变体：

示例：

#include <iostream>
#include <Eigen/Dense>

using namespace std;

int main()
{
  Eigen::Matrix3f m = Eigen::Matrix3f::Random();
  std::ptrdiff_t i, j;
  float minOfM = m.minCoeff(&i, &j);
  cout << "这是矩阵 m:\n" << m << endl;
  cout << "它的最小系数 (" << minOfM << ") 位于 (" << i << "," << j << ")\n\n";

  Eigen::RowVector4i v = Eigen::RowVector4i::Random();
  int maxOfV = v.maxCoeff(&i);
  cout << "这是向量 v: " << v << endl;
  cout << "它的最大系数 (" << maxOfV << ") 位于 " << i << endl;
}

输出：

这是矩阵 m:
  0.68  0.597  -0.33
-0.211  0.823  0.536
 0.566 -0.605 -0.444
它的最小系数 (-0.605) 位于 (2,1)

这是向量 v:  1  0  3 -3
它的最大系数 (3) 位于 2

7、运算的有效性

Eigen会检查所执行操作的有效性。如果可能，它会在编译时进行检查，产生编译错误。错误消息可能会很长，但Eigen会用大写字母突出重要信息。例如：

Matrix3f m;
Vector4f v;
v = m * v; // 编译时错误: YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES

当然，在许多情况下（例如检查动态大小），检查无法在编译时进行。此时，Eigen使用运行时断言。这意味着如果在“调试模式”下执行非法操作，程序会中止并显示错误消息，如果断言被关闭，程序可能会崩溃。

MatrixXf m(3, 3);
VectorXf v(4);
v = m * v; // 运行时断言失败：“无效的矩阵乘积”

通过这些基本操作，你可以高效地使用Eigen进行矩阵和向量的算术操作。

三、使用Eigen的Array类和系数操作

Eigen的Array类提供了通用数组，而Matrix类则用于线性代数。此外，Array类提供了一种简便的方法来执行系数操作，这些操作在数学上可能没有线性代数意义，例如为数组中的每个系数加上一个常数或按系数逐个相乘两个数组。

1、Array类型

Array是一个类模板，接受与Matrix相同的模板参数。与Matrix一样，前三个模板参数是必需的：

Array<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>

后三个模板参数是可选的。由于这些参数与Matrix完全相同，我们不再详细解释，详见Matrix类。

Eigen还提供了一些常见类型的typedef，这与Matrix的typedef类似，但有一些微小的区别，因为“array”用于一维和二维数组。我们采用约定，形式为ArrayNt的typedef表示一维数组，其中N和t分别表示大小和标量类型。对于二维数组，我们使用形式为ArrayNNt的typedef。以下是一些示例：

Array<float,Dynamic,1>

ArrayXf

Array<float,3,1>

Array3f

Array<double,Dynamic,Dynamic>

ArrayXXd

Array<double,3,3>

Array33d

2、访问Array中的值

括号操作符被重载，以提供对数组系数的读写访问，与矩阵类似。此外，可以使用<<操作符来初始化数组（通过逗号初始化器）或打印它们。

示例：

#include <Eigen/Dense>
#include <iostream>
 
int main()
{
  Eigen::ArrayXXf m(2,2);
  
  // 逐个系数赋值
  m(0,0) = 1.0; m(0,1) = 2.0;
  m(1,0) = 3.0; m(1,1) = m(0,1) + m(1,0);
  
  // 输出值到标准输出
  std::cout << m << std::endl << std::endl;
 
  // 使用逗号初始化器
  m << 1.0,2.0,
       3.0,4.0;
     
  // 输出值到标准输出
  std::cout << m << std::endl;
}

输出：

3、加法和减法

添加和减去两个数组与矩阵相同。如果两个数组具有相同的大小，操作是有效的，且按系数逐个进行加法或减法。数组还支持形式为array + scalar的表达式，将一个标量加到数组的每个系数上。这提供了矩阵对象无法直接实现的功能。

#include <Eigen/Dense>
#include <iostream>
 
int main()
{
  Eigen::ArrayXXf a(3,3);
  Eigen::ArrayXXf b(3,3);
  a << 1,2,3,
       4,5,6,
       7,8,9;
  b << 1,2,3,
       1,2,3,
       1,2,3;
       
  // 两个数组相加
  std::cout << "a + b = " << std::endl << a + b << std::endl << std::endl;
 
  // 从数组中减去标量
  std::cout << "a - 2 = " << std::endl << a - 2 << std::endl;
}

4、数组乘法

首先，您当然可以将数组与标量相乘，这与矩阵相同。数组与矩阵的本质区别在于两者相乘时的解释。矩阵将乘法解释为矩阵乘积，而数组将乘法解释为按系数逐个乘积。因此，只有在两个数组具有相同维度时，才能将它们相乘。

示例：

#include <Eigen/Dense>
#include <iostream>
 
int main()
{
  Eigen::ArrayXXf a(2,2);
  Eigen::ArrayXXf b(2,2);
  a << 1,2,
       3,4;
  b << 5,6,
       7,8;
  std::cout << "a * b = " << std::endl << a * b << std::endl;
}

5、其他系数操作

除了上面描述的加法、减法和乘法操作，Array类还定义了其他系数操作。例如，.abs()方法获取每个系数的绝对值，而.sqrt()计算系数的平方根。如果您有两个相同大小的数组，可以调用.min(.)来构造一个新数组，其系数是两个给定数组对应系数的最小值。这些操作在以下示例中进行了说明。

示例：

#include <Eigen/Dense>
#include <iostream>
 
int main()
{
  Eigen::ArrayXf a = Eigen::ArrayXf::Random(5);
  a *= 2;
  std::cout << "a =" << std::endl
            << a << std::endl;
  std::cout << "a.abs() =" << std::endl
            << a.abs() << std::endl;
  std::cout << "a.abs().sqrt() =" << std::endl
            << a.abs().sqrt() << std::endl;
  std::cout << "a.min(a.abs().sqrt()) =" << std::endl
            << a.min(a.abs().sqrt()) << std::endl;
}

输出：

a =
  1.36
-0.422
  1.13
  1.19
  1.65
a.abs() =
 1.36
0.422
 1.13
 1.19
 1.65
a.abs().sqrt() =
1.17
0.65
1.06
1.09
1.28
a.min(a.abs().sqrt()) =
  1.17
-0.422
  1.06
  1.09
  1.28

6、数组与矩阵的相互转换

什么时候应该使用Matrix类对象，什么时候应该使用Array类对象？您不能在数组上应用矩阵操作，或在矩阵上应用数组操作。因此，如果需要进行线性代数操作，如矩阵乘法，则应使用矩阵；如果需要进行系数操作，则应使用数组。然而，有时需要同时使用Matrix和Array操作。在这种情况下，您需要将矩阵转换为数组或反过来转换。这使您可以访问所有操作，而不管是将对象声明为数组还是矩阵。

矩阵表达式有一个.array()方法，将其“转换”为数组表达式，以便可以轻松应用系数操作。相反，数组表达式有一个.matrix()方法。与所有Eigen表达式抽象一样，这没有任何运行时成本（前提是让编译器进行优化）。当然.array()和.matrix()都可以作为右值和左值使用。

在表达式中混合矩阵和数组在Eigen中是禁止的。例如，不能直接将矩阵和数组相加；+运算符的操作数要么都是矩阵，要么都是数组。然而，可以使用.array()和.matrix()在两者之间轻松转换。此规则的例外是赋值操作：允许将矩阵表达式赋值给数组变量，或将数组表达式赋值给矩阵变量。

以下示例展示了如何使用.array()方法在Matrix对象上进行数组操作。例如，语句result = m.array() * n.array()将两个矩阵m和n转换为数组，按系数逐个相乘，并将结果赋值给矩阵变量result（这是合法的，因为Eigen允许将数组表达式赋值给矩阵变量）。

事实上，这种使用情况非常常见，Eigen为矩阵提供了.cwiseProduct(.)方法来计算系数乘积。这个方法也在示例程序中展示了。

示例：

#include <Eigen/Dense>
#include <iostream>
 
using Eigen::MatrixXf;
 
int main()
{
  MatrixXf m(2,2);
  MatrixXf n(2,2);
  MatrixXf result(2,2);
 
  m << 1,2,
       3,4;
  n << 5,6,
       7,8;
 
  result = m * n;
  std::cout << "-- Matrix m*n: --\n" << result << "\n\n";
  result = m.array() * n.array();
  std::cout << "-- Array m*n: --\n" << result << "\n\n";
  result = m.cwiseProduct(n);
  std::cout << "-- With cwiseProduct: --\n" << result << "\n\n";
  result = m.array() + 4;
  std::cout << "-- Array m + 4: --\n" << result << "\n\n";
}

输出：

-- Matrix m*n: --
19 22
43 50

-- Array m*n: --
 5 12
21 32

-- With cwiseProduct: --
 5 12
21 32

-- Array m + 4: --
5 6
7 8

同样，如果array1和array2是数组，则表达式array1.matrix() * array2.matrix()计算它们的矩阵乘积。

以下是一个更高级的示例。表达式(m.array() + 4).matrix() * m将4添加到矩阵m的每个系数，然后计算结果与m的矩阵乘积。类似地，表达式(m.array() * n.array()).matrix() * m计算矩阵m和n的系数乘积，然后计算结果与m的矩阵乘积。

示例：

#include <Eigen/Dense>
#include <iostream>
 
using Eigen::MatrixXf;
 
int main()
{
  MatrixXf m(2,2);
  MatrixXf n(2,2);
  MatrixXf result(2,2);
 
  m << 1,2,
       3,4;
  n << 5,6,
       7,8;
  
  result = (m.array() + 4).matrix() * m;
  std::cout << "-- Combination 1: --\n" << result << "\n\n";
  result = (m.array() * n.array()).matrix() * m;
  std::cout << "-- Combination 2: --\n" << result << "\n\n";
}

输出：

-- Combination 1: --
23 34
31 46

-- Combination 2: --
 41  58
117 170

四、块操作

1、矩阵和数组的块操作

块是矩阵或数组的矩形部分。块表达式可以用作右值，也可以用作左值。与Eigen表达式一样，只要让编译器进行优化，这种抽象就不会产生运行时成本。

Eigen中最通用的块操作叫做.block()。有两个版本，语法如下：

块操作	构造动态尺寸块表达式的版本	构造固定尺寸块表达式的版本
尺寸为(p,q)，起始位置为(i,j)的块	`matrix.block(i,j,p,q);`	`matrix.block<p,q>(i,j);`

如同在Eigen中一样，索引从0开始。

这两种版本都可以用于固定尺寸和动态尺寸的矩阵和数组。这两个表达式在语义上是等效的。唯一的区别是，如果块尺寸较小，固定尺寸版本通常会生成更快的代码，但要求在编译时已知该尺寸。

以下程序使用动态尺寸和固定尺寸版本打印矩阵内的几个块的值。

示例：

#include <Eigen/Dense>
#include <iostream>
 
using namespace std;
 
int main()
{
  Eigen::MatrixXf m(4,4);
  m <<  1, 2, 3, 4,
        5, 6, 7, 8,
        9,10,11,12,
       13,14,15,16;
  cout << "中间的块" << endl;
  cout << m.block<2,2>(1,1) << endl << endl;
  for (int i = 1; i <= 3; ++i)
  {
    cout << "大小为 " << i << "x" << i << " 的块" << endl;
    cout << m.block(0,0,i,i) << endl << endl;
  }
}

输出：

中间的块
 6  7
10 11

大小为 1x1 的块
1

大小为 2x2 的块
1 2
5 6

大小为 3x3 的块
 1  2  3
 5  6  7
 9 10 11

在上面的例子中，.block()函数用作右值，即它只是被读取。然而，块也可以用作左值，这意味着可以给块赋值。

这在下面的例子中进行了说明。这个例子还演示了数组中的块操作，其工作方式与矩阵中的块操作完全相同。

示例：

#include <Eigen/Dense>
#include <iostream>
 
int main()
{
  Eigen::Array22f m;
  m << 1,2,
       3,4;
  Eigen::Array44f a = Eigen::Array44f::Constant(0.6);
  std::cout << "这里是数组a:\n" << a << "\n\n";
  a.block<2,2>(1,1) = m;
  std::cout << "现在的a是将m复制到它的中央2x2块:\n" << a << "\n\n";
  a.block(0,0,2,3) = a.block(2,1,2,3);
  std::cout << "现在的a是将右下角的2x3块复制到左上角的2x3块:\n" << a << "\n\n";
}

输出：

这里是数组a:
0.6 0.6 0.6 0.6
0.6 0.6 0.6 0.6
0.6 0.6 0.6 0.6
0.6 0.6 0.6 0.6

现在的a是将m复制到它的中央2x2块:
0.6 0.6 0.6 0.6
0.6   1   2 0.6
0.6   3   4 0.6
0.6 0.6 0.6 0.6

现在的a是将右下角的2x3块复制到左上角的2x3块:
  3   4 0.6 0.6
0.6 0.6 0.6 0.6
0.6   3   4 0.6
0.6 0.6 0.6 0.6

虽然.block()方法可以用于任何块操作，但对于特殊情况有其他方法，提供了更专业的API和/或更好的性能。在性能方面，重要的是在编译时尽可能多地提供信息。例如，如果您的块是矩阵中的单独整列，则使用下面描述的专用.col()函数可以让Eigen知道，从而获得优化机会。

2、操作整列和整行

单个列和行是块的特殊情况。Eigen提供了方法轻松地访问它们：.col()和.row()。

块操作	方法
第i行	`matrix.row(i);`
第j列	`matrix.col(j);`

col()和row()的参数是要访问的列或行的索引。如同在Eigen中一样，索引从0开始。

示例：

#include <Eigen/Dense>
#include <iostream>
 
using namespace std;
 
int main()
{
  Eigen::MatrixXf m(3,3);
  m << 1,2,3,
       4,5,6,
       7,8,9;
  cout << "这里是矩阵m:" << endl << m << endl;
  cout << "第2行: " << m.row(1) << endl;
  m.col(2) += 3 * m.col(0);
  cout << "在第三列加上第一列的3倍之后，矩阵m变为:\n";
  cout << m << endl;
}

输出：

这里是矩阵m:
1 2 3
4 5 6
7 8 9
第2行: 4 5 6
在第三列加上第一列的3倍之后，矩阵m变为:
 1  2  6
 4  5 18
 7  8 30

该示例还演示了块表达式（此处为列）可以像其他任何表达式一样用于算术运算。

3、与角相关的操作

Eigen还提供了对矩阵或数组的一角或一侧的块进行操作的特殊方法。例如.topLeftCorner()可用于引用矩阵左上角的块。

不同的可能性总结在下表中：

块操作	构造动态尺寸块表达式的版本	构造固定尺寸块表达式的版本
左上角p乘q的块	`matrix.topLeftCorner(p,q);`	`matrix.topLeftCorner<p,q>();`
左下角p乘q的块	`matrix.bottomLeftCorner(p,q);`	`matrix.bottomLeftCorner<p,q>();`
右上角p乘q的块	`matrix.topRightCorner(p,q);`	`matrix.topRightCorner<p,q>();`
右下角p乘q的块	`matrix.bottomRightCorner(p,q);`	`matrix.bottomRightCorner<p,q>();`
包含前q行的块	`matrix.topRows(q);`	`matrix.topRows<q>();`
包含后q行的块	`matrix.bottomRows(q);`	`matrix.bottomRows<q>();`
包含前p列的块	`matrix.leftCols(p);`	`matrix.leftCols<p>();`
包含后q列的块	`matrix.rightCols(q);`	`matrix.rightCols<q>();`
包含从i开始的q列的块	`matrix.middleCols(i,q);`	`matrix.middleCols<q>(i);`
包含从i开始的q行的块	`matrix.middleRows(i,q);`	`matrix.middleRows<q>(i);`

下面是一个简单的例子，说明了上述操作的使用：

示例：

#include <Eigen/Dense>
#include <iostream>
 
using namespace std;
 
int main()
{
  Eigen::Matrix4f m;
  m << 1, 2, 3, 4,
       5, 6, 7, 8,
       9, 10,11,12,
       13,14,15,16;
  cout << "m.leftCols(2) =" << endl << m.leftCols(2) << endl << endl;
  cout << "m.bottomRows<2>() =" << endl << m.bottomRows<2>() << endl << endl;
  m.topLeftCorner(1,3) = m.bottomRightCorner(3,1).transpose();
  cout << "赋值后，m = " << endl << m << endl;
}

输出：

m.leftCols(2) =
 1  2
 5  6
 9 10
13 14

m.bottomRows<2>() =
 9 10 11 12
13 14 15 16

赋值后，m = 
 8 12 16  4
 5  6  7  8
 9 10 11 12
13 14 15 16

4、向量的块操作

Eigen还提供了一组专门为向量和一维数组设计的块操作：

块操作	构造动态尺寸块表达式的版本	构造固定尺寸块表达式的版本
包含前n个元素的块	`vector.head(n);`	`vector.head<n>();`
包含最后n个元素的块	`vector.tail(n);`	`vector.tail<n>();`
包含从位置i开始的n个元素的块	`vector.segment(i,n);`	`vector.segment<n>(i);`

下面是一个例子：

示例：

#include <Eigen/Dense>
#include <iostream>
 
using namespace std;
 
int main()
{
  Eigen::ArrayXf v(6);
  v << 1, 2, 3, 4, 5, 6;
  cout << "v.head(3) =" << endl << v.head(3) << endl << endl;
  cout << "v.tail<3>() = " << endl << v.tail<3>() << endl << endl;
  v.segment(1,4) *= 2;
  cout << "在'v.segment(1,4) *= 2'之后，v =" << endl << v << endl;
}

输出：

v.head(3) =
1
2
3

v.tail<3>() = 
4
5
6

在'v.segment(1,4) *= 2'之后，v =
 1
 4
 6
 8
10
 6

五、切片和索引

这一小节介绍了Eigen库中的切片和索引操作，这些操作可以用于子集行和列的选择。从Eigen 3.4版本开始引入了这些功能，它支持Block API提供的所有功能，并且更加灵活。特别地，它支持切片操作，这些操作可以在矩阵或向量中均匀间隔地选择行、列或元素。熟悉Python中pandas的同学可以类比pandas中的.iloc操作。

所有这些操作通过DenseBase::operator()(const RowIndices&, const ColIndices&)方法处理。每个参数可以是：

整数，用于索引单个行或列，包括符号索引。
符号Eigen::all，表示增序矩阵或向量的所有行或列。
ArithmeticSequence，由Eigen::seq、Eigen::seqN或Eigen::placeholders::lastN函数构造，用于表示等差数列。
任何整数的1D向量/数组，包括Eigen的向量/数组、表达式、std::vector、std::array以及普通的C数组：int[N]。

此外，它还接受任何实现以下两个成员函数的对象：

<integral type> operator[](<integral type>) const;
<integral type> size() const;

其中，<integral type>是与Eigen::Index兼容的任何整数类型（例如std::ptrdiff_t）。

基本切片

通过Eigen::seq或Eigen::seqN函数可以在矩阵或向量中均匀间隔地选择行、列或元素。它们的签名如下所示：

函数	描述	示例
`seq(firstIdx,lastIdx)`	表示从`firstIdx`到`lastIdx`的整数序列	`seq(2,5)` <=> {2,3,4,5}
`seq(firstIdx,lastIdx,incr)`	使用增量`incr`从一个索引到下一个索引	`seq(2,8,2)` <=> {2,4,6,8}
`seqN(firstIdx,size)`	表示从`firstIdx`开始的`size`个整数序列	`seqN(2,5)` <=> {2,3,4,5,6}
`seqN(firstIdx,size,incr)`	使用增量`incr`从一个索引到下一个索引	`seqN(2,3,3)` <=> {2,5,8}

其中，firstIdx和lastIdx参数还可以使用Eigen::last符号，表示基础矩阵/向量的最后一行、列或元素的索引。

以下是一些示例，展示了如何使用切片操作来选择矩阵或向量中的子集：

· 选择从第i行开始，有n列的左下角块：

A(seq(i,Eigen::last), seqN(0,n))

等同于：

A.bottomLeftCorner(A.rows()-i,n)

· 选择从第i行、第j列开始，大小为m行n列的块：

A(seqN(i,m), seqN(j,n))

相当于：

A.block(i,j,m,n)

· 选择从(i0, j0) 到 (i1, j1)的块：

A(seq(i0,i1), seq(j0,j1))

相当于：

A.block(i0,j0,i1-i0+1,j1-j0+1)

· 选择A的偶数列：

A(all, seq(0,Eigen::last,2))

· 选择A的前n个奇数行：

A(seqN(1,n,2), all)

Eigen的切片和索引功能通过DenseBase::operator()方法提供了强大的灵活性，可以用于任何行、列或元素的选择，支持符号索引、等差数列、最后索引符号以及各种整数向量和数组。这些功能使得对矩阵和向量进行复杂的选择和操作变得简单和高效。

六、高级初始化

本节讨论了几种初始化矩阵的高级方法，详细介绍了之前引入的逗号初始化器。同时解释了如何获取特殊的矩阵，如单位矩阵和零矩阵。

1、逗号初始化器

Eigen提供了逗号初始化器语法，允许用户轻松设置矩阵、向量或数组的所有系数。只需按照从左上角开始从左到右、从上到下的顺序列出系数即可。在使用逗号初始化器之前需要指定对象的大小。如果列出的系数太少或太多，Eigen会报错。

示例:

Matrix3f m;
m << 1, 2, 3,
     4, 5, 6,
     7, 8, 9;
std::cout << m;

输出：

1 2 3
4 5 6
7 8 9

此外，初始列表中的元素本身也可以是向量或矩阵，常见用法是将两个行向量连接在一起：

RowVectorXd vec1(3);
vec1 << 1, 2, 3;

RowVectorXd vec2(4);
vec2 << 1, 4, 9, 16;

RowVectorXd joined(7);
joined << vec1, vec2;
std::cout << joined;

输出：

1 2 3 1 4 9 16

同样的技术可以用于初始化具有块结构的矩阵：

MatrixXf matA(2, 2);
matA << 1, 2,
        3, 4;

MatrixXf matB(4, 4);
matB << matA, matA/10,
        matA/10, matA;

std::cout << matB << std::endl;

输出：

1   2  0.1 0.2
3   4  0.3 0.4
0.1 0.2  1   2
0.3 0.4  3   4

逗号初始化器还可以用于填充块表达式，如m.row(i)。

2、特殊矩阵和数组

Matrix和Array类具有静态方法，如Zero()，用于将所有系数初始化为零。有三个变体：

第一个变体不带参数，只能用于固定大小对象。
第二个变体需要一个参数，用于一维动态大小对象。
第三个变体需要两个参数，用于二维对象。

ArrayXXf a1 = ArrayXXf::Zero(3, 4);
std::cout << a1 << std::endl;

输出：

0 0 0 0
0 0 0 0
0 0 0 0

类似地，静态方法Constant(value)将所有系数设置为指定的值。如果需要指定对象的大小，额外的参数应该在值参数之前传入，例如MatrixXd::Constant(rows, cols, value)：

MatrixXd mat = MatrixXd::Constant(3, 3, 5.0);
std::cout << mat << std::endl;

输出：

5 5 5
5 5 5
5 5 5

总的来说，Eigen提供了强大的初始化和操作功能，可以通过逗号初始化器轻松设置矩阵、向量和数组的系数。此外，它还提供了方便的静态方法来生成特殊的矩阵，如零矩阵、常数矩阵以及随机矩阵。这些功能使得在线性代数和数值计算中快速、高效地使用Eigen变得可能。

七、Eigen的归约、访问器和广播机制

1、归约

在Eigen库中，归约是一种将矩阵或数组转换为单一标量值的函数。一个常用的归约函数是 .sum()，它返回矩阵或数组中所有系数的和。

示例：

#include <iostream>
#include <Eigen/Dense>

using namespace std;
int main()
{
  Eigen::Matrix2d mat;
  mat << 1, 2,
         3, 4;
  cout << "Here is mat.sum():       " << mat.sum()       << endl;
  cout << "Here is mat.prod():      " << mat.prod()      << endl;
  cout << "Here is mat.mean():      " << mat.mean()      << endl;
  cout << "Here is mat.minCoeff():  " << mat.minCoeff()  << endl;
  cout << "Here is mat.maxCoeff():  " << mat.maxCoeff()  << endl;
  cout << "Here is mat.trace():     " << mat.trace()     << endl;
}

输出：

Here is mat.sum():       10
Here is mat.prod():      24
Here is mat.mean():      2.5
Here is mat.minCoeff():  1
Here is mat.maxCoeff():  4
Here is mat.trace():     5

2、范数计算

向量的（欧几里得，也叫ℓ2）平方范数可以通过 squaredNorm() 获得。它等于向量自身的点积，或者等于系数的平方和。Eigen还提供了 norm() 方法，返回 squaredNorm() 的平方根。

这些操作也可以作用于矩阵，在这种情况下，一个 n-by-p 矩阵被视为大小为 (n*p) 的向量。例如，norm() 方法返回“Frobenius”或“Hilbert-Schmidt”范数。如果需要其他系数-wise的 ℓp 范数，可以使用 lpNorm<p>() 方法。模板参数 p 可以取特殊值 Infinity ，表示 ℓ∞ 范数，即系数的绝对值的最大值。

示例：

#include <Eigen/Dense>
#include <iostream>

int main()
{
  Eigen::VectorXf v(2);
  Eigen::MatrixXf m(2,2);
  
  v << -1, 2;
  
  m << 1, -2,
       -3, 4;
 
  std::cout << "v.squaredNorm() = " << v.squaredNorm() << std::endl;
  std::cout << "v.norm() = " << v.norm() << std::endl;
  std::cout << "v.lpNorm<1>() = " << v.lpNorm<1>() << std::endl;
  std::cout << "v.lpNorm<Infinity>() = " << v.lpNorm<Eigen::Infinity>() << std::endl;
 
  std::cout << std::endl;
  std::cout << "m.squaredNorm() = " << m.squaredNorm() << std::endl;
  std::cout << "m.norm() = " << m.norm() << std::endl;
  std::cout << "m.lpNorm<1>() = " << m.lpNorm<1>() << std::endl;
  std::cout << "m.lpNorm<Infinity>() = " << m.lpNorm<Eigen::Infinity>() << std::endl;
}

输出：

v.squaredNorm() = 5
v.norm() = 2.23607
v.lpNorm<1>() = 3
v.lpNorm<Infinity>() = 2

m.squaredNorm() = 30
m.norm() = 5.47723
m.lpNorm<1>() = 10
m.lpNorm<Infinity>() = 4

3、布尔归约

以下归约操作用于布尔值：

all() 如果所有系数都为 true，返回 true。
any() 如果至少有一个系数为 true，返回 true。
count() 返回系数为 true 的数量。

这些通常与系数-wise比较和相等操作符结合使用。例如，array > 0 是一个与 array 同大小的数组，其中位置为 true 的对应系数大于 0。因此，(array > 0).all() 测试数组的所有系数是否大于 0。

示例：

#include <Eigen/Dense>
#include <iostream>

int main()
{
  Eigen::ArrayXXf a(2,2);
  
  a << 1, 2,
       3, 4;
 
  std::cout << "(a > 0).all()   = " << (a > 0).all() << std::endl;
  std::cout << "(a > 0).any()   = " << (a > 0).any() << std::endl;
  std::cout << "(a > 0).count() = " << (a > 0).count() << std::endl;
  std::cout << std::endl;
  std::cout << "(a > 2).all()   = " << (a > 2).all() << std::endl;
  std::cout << "(a > 2).any()   = " << (a > 2).any() << std::endl;
  std::cout << "(a > 2).count() = " << (a > 2).count() << std::endl;
}

输出：

(a > 0).all()   = 1
(a > 0).any()   = 1
(a > 0).count() = 4

(a > 2).all()   = 0
(a > 2).any()   = 1
(a > 2).count() = 2

4、访问器

访问器在需要获取矩阵或数组中某个系数的位置时很有用。最简单的示例是 maxCoeff(&x,&y) 和 minCoeff(&x,&y)，它们可以用来找到矩阵或数组中最大或最小系数的位置。

示例：

#include <iostream>
#include <Eigen/Dense>

int main()
{
  Eigen::MatrixXf m(2,2);
  
  m << 1, 2,
       3, 4;
 
  // 获取最大值的位置
  Eigen::Index maxRow, maxCol;
  float max = m.maxCoeff(&maxRow, &maxCol);
 
  // 获取最小值的位置
  Eigen::Index minRow, minCol;
  float min = m.minCoeff(&minRow, &minCol);
 
  std::cout << "Max: " << max <<  ", at: " <<
     maxRow << "," << maxCol << std::endl;
  std::cout << "Min: " << min << ", at: " <<
     minRow << "," << minCol << std::endl;
}

输出：

Max: 4, at: 1,1
Min: 1, at: 0,0

5、部分归约

看到这里有的同学可能会想要，为什么没有部分按行或者按列计算。别着急，接下来就要介绍这种功能。部分归约可以通过 colwise() 或 rowwise() 实现。

示例：

#include <iostream>
#include <Eigen/Dense>

using namespace std;
int main()
{
  Eigen::MatrixXf mat(2,4);
  mat << 1, 2, 6, 9,
         3, 1, 7, 2;
  
  std::cout << "Column's maximum: " << std::endl
   << mat.colwise().maxCoeff() << std::endl;
}

输出：

Column's maximum: 
3 2 7 9

同样的操作也可以按行进行：

#include <iostream>
#include <Eigen/Dense>

using namespace std;
int main()
{
  Eigen::MatrixXf mat(2,4);
  mat << 1, 2, 6, 9,
         3, 1, 7, 2;
  
  std::cout << "Row's maximum: " << std::endl
   << mat.rowwise().maxCoeff() << std::endl;
}

输出：

Row's maximum: 
9
7

6、广播

广播的概念类似于部分归约，不同之处在于广播构建一个表达式，其中一个向量（列或行）通过在一个方向上复制它来解释为矩阵。

#include <iostream>
#include <Eigen/Dense>

using namespace std;
int main()
{
  Eigen::MatrixXf mat(2,4);
  Eigen::VectorXf v(2);
  
  mat << 1, 2, 6, 9,
         3, 1, 7, 2;
         
  v << 0,
       1;
  
  mat.colwise() += v;
  
  cout << mat << endl;
}

输出：

1 2 6 9
4 2 8 3

同样的，广播也可以按列进行：

#include <iostream>
#include <Eigen/Dense>

using namespace std;
int main()
{
  Eigen::MatrixXf mat(2,4);
  Eigen::RowVectorXf v(4);
  
  mat << 1, 2, 6, 9,
         3, 1, 7, 2;
         
  v << 0, 1, 2, 3;
  
  mat.rowwise() += v;
  
  cout << mat << endl;
}

输出：

1 3 8 12
3 2 9 5

八、重塑

从版本3.4开始，Eigen提供了方便的方法将矩阵重塑为不同大小的矩阵或向量。所有情况都通过DenseBase::reshaped(NRowsType, NColsType)和DenseBase::reshaped()函数处理。这些函数不执行就地重塑，而是返回输入表达式的视图。

1、二维视图

更常见的重塑转换通过reshaped(nrows, ncols)处理。下面是一个将4x4矩阵重塑为2x8矩阵的示例：

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::Matrix4i m = Eigen::Matrix4i::Random();
    std::cout << "Here is the matrix m:" << std::endl << m << std::endl;
    std::cout << "Here is m.reshaped(2, 8):" << std::endl << m.reshaped(2, 8) << std::endl;
}

输出：

Here is the matrix m:
 7  9 -5 -3
-2 -6  1  0
 6 -3  0  9
 6  6  3  9
Here is m.reshaped(2, 8):
 7  6  9 -3 -5  0 -3  9
-2  6 -6  6  1  3  0  9

默认情况下，无论输入表达式的存储顺序如何，输入系数总是按照列主要顺序解释。

2、一维线性视图

重塑的一个非常常见的用法是为给定的二维矩阵或表达式创建一维线性视图。在这种情况下，可以省略尺寸，如以下示例所示：

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::Matrix4i m = Eigen::Matrix4i::Random();
    std::cout << "Here is the matrix m:" << std::endl << m << std::endl;
    std::cout << "Here is m.reshaped().transpose():" << std::endl << m.reshaped().transpose() << std::endl;
    std::cout << "Here is m.reshaped<RowMajor>().transpose():" << std::endl << m.reshaped<Eigen::RowMajor>().transpose() << std::endl;
}

输出：

Here is the matrix m:
 7  9 -5 -3
-2 -6  1  0
 6 -3  0  9
 6  6  3  9
Here is m.reshaped().transpose():
 7 -2  6  6  9 -6 -3  6 -5  1  0  3 -3  0  9  9
Here is m.reshaped<RowMajor>().transpose():  
 7  9 -5 -3 -2 -6  1  0  6 -3  0  9  6  6  3  9

这个快捷方式总是返回一个列向量，默认情况下输入系数总是按照列主要顺序解释。

3、就地重塑

上述示例创建了重塑视图，但如何就地重塑给定矩阵呢？当然，这项任务仅适用于动态尺寸的矩阵和数组。在许多情况下，这可以通过PlainObjectBase::resize(Index, Index)实现：

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::MatrixXi m = Eigen::Matrix4i::Random();
    std::cout << "Here is the matrix m:" << std::endl << m << std::endl;
    std::cout << "Here is m.reshaped(2, 8):" << std::endl << m.reshaped(2, 8) << std::endl;
    m.resize(2,8);
    std::cout << "Here is the matrix m after m.resize(2,8):" << std::endl << m << std::endl;
}

输出：

Here is the matrix m:
 7  9 -5 -3
-2 -6  1  0
 6 -3  0  9
 6  6  3  9
Here is m.reshaped(2, 8):
 7  6  9 -3 -5  0 -3  9
-2  6 -6  6  1  3  0  9
Here is the matrix m after m.resize(2,8):
 7  6  9 -3 -5  0 -3  9
-2  6 -6  6  1  3  0  9

但是要注意，与reshaped不同，resize的结果取决于输入存储顺序。它的行为类似于reshaped<AutoOrder>：

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> m = Eigen::Matrix4i::Random();
    std::cout << "Here is the matrix m:" << std::endl << m << std::endl;
    std::cout << "Here is m.reshaped(2, 8):" << std::endl << m.reshaped(2, 8) << std::endl;
    std::cout << "Here is m.reshaped<AutoOrder>(2, 8):" << std::endl << m.reshaped<Eigen::AutoOrder>(2, 8) << std::endl;
    m.resize(2,8);
    std::cout << "Here is the matrix m after m.resize(2,8):" << std::endl << m << std::endl;
}

输出：

Here is the matrix m:
 7 -2  6  6
 9 -6 -3  6
-5  1  0  3
-3  0  9  9
Here is m.reshaped(2, 8):
 7 -5 -2  1  6  0  6  3
 9 -3 -6  0 -3  9  6  9
Here is m.reshaped<AutoOrder>(2, 8):
 7 -2  6  6  9 -6 -3  6
-5  1  0  3 -3  0  9  9
Here is the matrix m after m.resize(2,8):
 7 -2  6  6  9 -6 -3  6
-5  1  0  3 -3  0  9  9

最后，将重塑矩阵赋值给自身目前是不支持的，因为会导致别名问题，结果是未定义行为。以下是不允许的：

A = A.reshaped(2,8);

可以这样做：

A = A.reshaped(2,8).eval();

九、STL迭代器和算法

1、稠密矩阵和数组操作

从版本3.4开始，Eigen的稠密矩阵和数组提供了与STL兼容的迭代器。如下所示，这使得它们与range-for循环和STL算法自然兼容。

2、遍历一维数组和向量

任何一维表达式都提供了一对begin()/end()方法来遍历它们。

这直接支持C++11的range for循环：

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::VectorXi v = Eigen::VectorXi::Random(4);
    std::cout << "Here is the vector v:\n";
    for(auto x : v) std::cout << x << " ";
    std::cout << "\n";
}

输出：

Here is the vector v:
7 -2 6 6

一维表达式也可以轻松传递给STL算法：

#include <Eigen/Dense>
#include <iostream>
#include <algorithm>

int main() {
    Eigen::Array4i v = Eigen::Array4i::Random().abs();
    std::cout << "Here is the initial vector v:\n" << v.transpose() << "\n";
    std::sort(v.begin(), v.end());
    std::cout << "Here is the sorted vector v:\n" << v.transpose() << "\n";
}

输出：

Here is the initial vector v:
7 2 6 6
Here is the sorted vector v:
2 6 6 7

类似于std::vector，一维表达式也提供了一对cbegin()/cend()方法，以便在非const对象上方便地获取const迭代器。

3、遍历二维数组和矩阵的系数

STL迭代器本质上是为遍历一维结构设计的。这就是为什么begin()/end()方法对于二维表达式是禁用的。通过重塑为一维线性视图，仍然可以轻松遍历二维表达式的所有系数：

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::Matrix2i A = Eigen::Matrix2i::Random();
    std::cout << "Here are the coeffs of the 2x2 matrix A:\n";
    for(auto x : A.reshaped())
        std::cout << x << " ";
    std::cout << "\n";
}

输出：

Here are the coeffs of the 2x2 matrix A:
7 -2 6 6

4、遍历二维数组和矩阵的行或列

还可以获取二维表达式的行或列的迭代器。这些迭代器通过rowwise()和colwise()代理提供。下面是一个示例，排序矩阵的每一行：

#include <Eigen/Dense>
#include <iostream>
#include <algorithm>

int main() {
    Eigen::ArrayXXi A = Eigen::ArrayXXi::Random(4,4).abs();
    std::cout << "Here is the initial matrix A:\n" << A << "\n";
    for(auto row : A.rowwise())
        std::sort(row.begin(), row.end());
    std::cout << "Here is the sorted matrix A:\n" << A << "\n";
}

输出：

Here is the initial matrix A:
7 9 5 3
2 6 1 0
6 3 0 9
6 6 3 9
Here is the sorted matrix A:
3 5 7 9
0 1 2 6
0 3 6 9
3 6 6 9

十、使用原始缓冲区：Map类

本节解释了如何使用“原始” C/C++数组。这在多种情况下都非常有用，特别是在将向量和矩阵从其他库导入到Eigen时。

有时你可能会有一个预定义的数字数组，想在 Eigen 中将其用作向量或矩阵。虽然一种选择是复制这些数据，但通常你可能想将此内存重用为 Eigen 类型。幸运的是，使用 Map 类可以很容易地实现这一点。

1、Map 类型和声明 Map 变量

一个 Map 对象的类型由其对应的 Eigen 类型定义：

Map<Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime> >

注意，在默认情况下，Map 只需要一个模板参数。要构造一个 Map 变量，你需要两个其他信息：指向定义系数数组的内存区域的指针，以及所需的矩阵或向量的形状。例如，要定义一个编译时确定大小的浮点矩阵，可以这样做：

Map<MatrixXf> mf(pf, rows, columns);

其中 pf 是一个指向内存数组的 float *。一个固定大小的只读整数向量可以声明为：

Map<const Vector4i> mi(pi);

其中 pi 是一个 int *。在这种情况下，不必将大小传递给构造函数，因为它已由 Matrix/Array 类型指定。注意，Map 没有默认构造函数；必须传递一个指针来初始化对象。

Map 足够灵活，可以适应各种不同的数据表示。还有两个其他（可选）的模板参数：

Map<typename MatrixType, int MapOptions, typename StrideType>

MapOptions 指定指针是对齐的还是未对齐的。默认值是未对齐。StrideType 允许你使用 Stride 类为内存数组指定自定义布局。例如，可以指定数据数组是以行优先格式组织的：

#include <Eigen/Dense>
#include <iostream>

int main() {
    int array[8];
    for(int i = 0; i < 8; ++i) array[i] = i;
    std::cout << "Column-major:\n" << Eigen::Map<Eigen::Matrix<int,2,4> >(array) << std::endl;
    std::cout << "Row-major:\n" << Eigen::Map<Eigen::Matrix<int,2,4,Eigen::RowMajor> >(array) << std::endl;
    std::cout << "Row-major using stride:\n" <<
      Eigen::Map<Eigen::Matrix<int,2,4>, Eigen::Unaligned, Eigen::Stride<1,4> >(array) << std::endl;
}

输出：

Column-major:
0 2 4 6
1 3 5 7
Row-major:
0 1 2 3
4 5 6 7
Row-major using stride:
0 1 2 3
4 5 6 7

2、使用 Map 变量

你可以像使用其他 Eigen 类型一样使用 Map 对象：

#include <Eigen/Dense>
#include <iostream>

int main() {
    typedef Eigen::Matrix<float,1,Eigen::Dynamic> MatrixType;
    typedef Eigen::Map<MatrixType> MapType;
    typedef Eigen::Map<const MatrixType> MapTypeConst;   // 只读映射
    const int n_dims = 5;
    
    MatrixType m1(n_dims), m2(n_dims);
    m1.setRandom();
    m2.setRandom();
    float *p = &m2(0);  // 获取存储 m2 数据的地址
    MapType m2map(p, m2.size());   // m2map 与 m2 共享数据
    MapTypeConst m2mapconst(p, m2.size());  // m2 的只读访问器
    
    std::cout << "m1: " << m1 << std::endl;
    std::cout << "m2: " << m2 << std::endl;
    std::cout << "Squared euclidean distance: " << (m1 - m2).squaredNorm() << std::endl;
    std::cout << "Squared euclidean distance, using map: " <<
      (m1 - m2map).squaredNorm() << std::endl;
    m2map(3) = 7;   // 这将改变 m2，因为它们共享同一个数组
    std::cout << "Updated m2: " << m2 << std::endl;
    std::cout << "m2 coefficient 2, constant accessor: " << m2mapconst(2) << std::endl;
    /* m2mapconst(2) = 5; */   // 这会导致编译时错误
}

输出：

m1:   0.68 -0.211  0.566  0.597  0.823
m2: -0.605  -0.33  0.536 -0.444  0.108
Squared euclidean distance: 3.26
Squared euclidean distance, using map: 3.26
Updated m2: -0.605  -0.33  0.536      7  0.108
m2 coefficient 2, constant accessor: 0.536

所有 Eigen 函数都编写为接受 Map 对象，就像其他 Eigen 类型一样。然而，在编写自己的接受 Eigen 类型的函数时，这不会自动发生：Map 类型与其 Dense 不同。

3、更改映射数组

可以在声明 Map 对象后使用 C++ 的“placement new”语法更改 Map 对象的数组：

#include <Eigen/Dense>
#include <iostream>

int main() {
    int data[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
    Eigen::Map<Eigen::RowVectorXi> v(data, 4);
    std::cout << "The mapped vector v is: " << v << "\n";
    new (&v) Eigen::Map<Eigen::RowVectorXi>(data + 4, 5);
    std::cout << "Now v is: " << v << "\n";
}

输出：

The mapped vector v is: 1 2 3 4
Now v is: 5 6 7 8 9

尽管看起来如此，但这并不会调用内存分配器，因为语法指定了存储结果的位置。

这种语法使得在声明 Map 对象时无需先知道映射数组在内存中的位置：

Map<Matrix3f> A(NULL);  // 不要试图使用这个矩阵！
VectorXf b(n_matrices);
for (int i = 0; i < n_matrices; i++) {
  new (&A) Map<Matrix3f>(get_matrix_pointer(i));
  b(i) = A.trace();
}

十一、混叠

在Eigen中，混叠指的是在赋值语句中，同一个矩阵（或数组或向量）同时出现在赋值操作符的左侧和右侧。例如，语句 mat = 2 * mat; 或 mat = mat.transpose(); 都涉及混叠。第一个例子中的混叠是无害的，但第二个例子中的混叠会导致意外结果。本页面解释了什么是混叠，何时有害，以及如何处理它。

下面是一个简单的混叠示例：

Example	Output
MatrixXi mat(3,3); 
mat << 1, 2, 3,   4, 5, 6,   7, 8, 9;
cout << "这里是矩阵mat:\n" << mat << endl;
 
// 这个赋值语句展示了混叠问题
mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);
cout << "赋值后，mat = \n" << mat << endl;

这里是矩阵mat:
1 2 3
4 5 6
7 8 9
赋值后，mat = 
1 2 3
4 1 2
7 4 1

输出的结果并不是预期的。问题在于赋值语句：

mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);

这个赋值语句表现了混叠：系数 mat(1,1) 同时出现在赋值操作符左侧的 mat.bottomRightCorner(2,2) 块和右侧的 mat.topLeftCorner(2,2) 块中。赋值后，右下角的 (2,2) 位置应具有赋值前 mat(1,1) 的值，即5。然而，输出显示 mat(2,2) 实际上是1。问题在于Eigen对 mat.topLeftCorner(2,2) 使用了惰性求值。结果类似于：

mat(1,1) = mat(0,0);
mat(1,2) = mat(0,1);
mat(2,1) = mat(1,0);
mat(2,2) = mat(1,1);

因此，mat(2,2) 被赋值为新的 mat(1,1) 值，而不是旧值。下一节解释了如何通过调用 eval() 解决这个问题。

在尝试缩小矩阵时，混叠问题更自然地发生。例如，表达式 vec = vec.head(n) 和 mat = mat.block(i,j,r,c) 表现了混叠。

通常，混叠无法在编译时检测到：如果第一个例子中的 mat 稍大一些，那么块就不会重叠，也不会有别名问题。然而，Eigen确实在运行时检测到了一些别名实例。以下示例展示了在矩阵和向量算术中的混叠：

Example	Output
Matrix2i a; a << 1, 2, 3, 4;
cout << "这里是矩阵a:\n" << a << endl;
 
a = a.transpose(); // !!! 不要这样做 !!!
cout << "混叠效应的结果:\n" << a << endl;

这里是矩阵a:
1 2
3 4
混叠效应的结果:
1 2
2 4

同样，输出展示了混叠问题。然而，默认情况下，Eigen使用运行时断言检测此问题，并退出并显示类似以下消息：

void Eigen::DenseBase<Derived>::checkTransposeAliasing(const OtherDerived&) const 
[with OtherDerived = Eigen::Transpose<Eigen::Matrix<int, 2, 2, 0, 2, 2> >, Derived = Eigen::Matrix<int, 2, 2, 0, 2, 2>]: 
Assertion `(!internal::check_transpose_aliasing_selector<Scalar,internal::blas_traits<Derived>::IsTransposed,OtherDerived>::run(internal::extract_data(derived()), other)) 
&& "aliasing detected during transposition, use transposeInPlace() or evaluate the rhs into a temporary using .eval()"' failed.

用户可以通过定义 EIGEN_NO_DEBUG 宏来关闭Eigen的运行时断言（如用于检测混叠问题的断言），并且上述程序是用关闭此宏编译的，以便展示混叠问题。

1、解决混叠问题

如果你理解混叠问题的原因，那么解决它的办法就显而易见了：Eigen必须将右侧表达式完全评估到一个临时矩阵/数组中，然后再赋值给左侧。eval() 函数正是这样做的。

例如，下面是第一个例子的更正版本：

Example	Output
MatrixXi mat(3,3); 
mat << 1, 2, 3,   4, 5, 6,   7, 8, 9;
cout << "这里是矩阵mat:\n" << mat << endl;
 
// eval() 解决了混叠问题
mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2).eval();
cout << "赋值后，mat = \n" << mat << endl;

这里是矩阵mat:
1 2 3
4 5 6
7 8 9
赋值后，mat = 
1 2 3
4 1 2
7 4 5

现在，赋值后 mat(2,2) 的值为5，如期望的一样。

同样的解决方法也适用于第二个例子中的转置：只需将 a = a.transpose(); 替换为 a = a.transpose().eval();。然而，在这种常见情况下，有更好的解决方案。Eigen提供了专门的 transposeInPlace() 函数，用于将矩阵替换为其转置。如下所示：

Example	Output
MatrixXf a(2,3); a << 1, 2, 3, 4, 5, 6;
cout << "这里是初始矩阵a:\n" << a << endl;
 
 
a.transposeInPlace();
cout << "转置后:\n" << a << endl;

这里是初始矩阵a:
1 2 3
4 5 6
转置后:
1 4
2 5
3 6

如果有 xxxInPlace() 函数可用，那么最好使用它，因为它能更清晰地表示你正在做什么。这也可能允许Eigen更积极地优化。以下是一些提供的 xxxInPlace() 函数：

原始函数	就地函数
MatrixBase::adjoint()	MatrixBase::adjointInPlace()
DenseBase::reverse()	DenseBase::reverseInPlace()
LDLT::solve()	LDLT::solveInPlace()
LLT::solve()	LLT::solveInPlace()
TriangularView::solve()	TriangularView::solveInPlace()
DenseBase::transpose()	DenseBase::transposeInPlace()

在使用类似 vec = vec.head(n) 表达式缩小矩阵或向量的特殊情况下，你可以使用 conservativeResize()。

2、混叠与组件操作

如上所述，如果在赋值操作符的左右两侧都出现了相同的矩阵或数组，则可能会出现危险情况，并且通常需要显式地评估右侧表达式。然而，应用组件操作（如矩阵加法、标量乘法和数组乘法）是安全的。

以下示例只有组件操作。因此，即使在赋值的两侧出现相同的矩阵，也不需要 eval()。

Example	Output
MatrixXf mat(2,2); 
mat << 1, 2,  4, 7;
cout << "这里是矩阵mat:\n" << mat << endl << endl;
 
mat = 2 * mat;
cout << "执行 'mat = 2 * mat' 后，mat = \n" << mat << endl << endl;
 
 
mat = mat - MatrixXf::Identity(2,2);
cout << "减法后，它变成\n" << mat << endl << endl;
 
 
ArrayXXf arr = mat;
arr = arr.square();
cout << "平方后，它变成\n" << arr << endl << endl;
 
// 将所有操作合并为一个语句：
mat << 1, 2,  4, 7;
mat = (2 * mat - MatrixXf::Identity(2,2)).array().square();
cout << "所有操作一次完成，结果是\n" << mat << endl << endl;

这里是矩阵mat:
1 2
4 7

执行 'mat = 2 * mat' 后，mat = 
 2  4
 8 14

减法后，它变成
 1  4
 8 13

平方后，它变成
  1  16
 64 169

所有操作一次完成，结果是
  1  16
 64 169

一般来说，如果右侧表达式的 (i,j) 条目只依赖于左侧矩阵或数组的 (i,j) 条目，而不依赖于任何其他条目，则赋值是安全的。在这种情况下，不需要显式地评估右侧表达式。

3、混叠与矩阵乘法

矩阵乘法是Eigen中唯一默认假定存在混叠的操作，前提是目标矩阵没有被重新调整大小。因此，如果 matA 是一个方矩阵，那么语句 matA = matA * matA; 是安全的。Eigen中的所有其他操作都假定没有混叠问题，要么因为结果被赋值给不同的矩阵，要么因为它是组件操作。

Example	Output
MatrixXf matA(2,2); 
matA << 2, 0,  0, 2;
matA = matA * matA;
cout << matA;

4 0
0 4

然而，这有一定代价。当执行表达式 matA = matA * matA 时，Eigen在一个临时矩阵中评估乘积，计算后将其赋值给 matA。这没问题。但Eigen在乘积被赋值给不同矩阵时也会这样做（例如，matB = matA * matA）。在这种情况下，直接将乘积评估到 matB 中比先评估到临时矩阵然后复制到 matB 更有效。

用户可以使用 noalias() 函数指示没有混叠，如下所示：matB.noalias() = matA * matA。这允许Eigen将矩阵乘积 matA * matA 直接评估到 matB 中。

Example	Output
MatrixXf matA(2,2), matB(2,2); 
matA << 2, 0,  0, 2;
 
// 简单但效率不高
matB = matA * matA;
cout << matB << endl << endl;
 
// 更复杂但更有效
matB.noalias() = matA * matA;
cout << matB;

当然，当实际存在混叠时，不应使用 noalias()。否则可能会得到错误的结果：

Example	Output
MatrixXf matA(2,2); 
matA << 2, 0,  0, 2;
matA.noalias() = matA * matA;
cout << matA;

4 0
0 4

此外，从Eigen 3.3开始，如果目标矩阵被重新调整大小并且乘积未直接赋值给目标矩阵，则不假定存在混叠。因此，以下示例也是错误的：

Example	Output
MatrixXf A(2,2), B(3,2);
B << 2, 0,  0, 3, 1, 1;
A << 2, 0, 0, -2;
A = (B * A).c

十二、储存顺序

对于矩阵和二维数组，有两种不同的存储顺序：列优先（column-major）和行优先（row-major）。本页面解释了这些存储顺序以及如何指定使用哪种存储顺序。

1、列优先和行优先存储

矩阵的元素形成一个二维网格。然而，当矩阵存储在内存中时，元素必须以某种方式线性排列。有两种主要方式：按行和按列。

下面的示例通过Eigen代码演示了这一点。它使用 PlainObjectBase::data() 函数返回矩阵第一个元素的内存位置指针。

Example	Output
Matrix<int, 3, 4, ColMajor> Acolmajor;
Acolmajor << 8, 2, 2, 9,
             9, 1, 4, 4,
             3, 5, 4, 5;
cout << "矩阵 A:" << endl;
cout << Acolmajor << endl << endl; 
 
cout << "内存中的存储顺序（列优先）:" << endl;
for (int i = 0; i < Acolmajor.size(); i++)
  cout << *(Acolmajor.data() + i) << "  ";
cout << endl << endl;
 
Matrix<int, 3, 4, RowMajor> Arowmajor = Acolmajor;
cout << "内存中的存储顺序（行优先）:" << endl;
for (int i = 0; i < Arowmajor.size(); i++)
  cout << *(Arowmajor.data() + i) << "  ";
cout << endl;

矩阵 A:
8 2 2 9
9 1 4 4
3 5 4 5

内存中的存储顺序（列优先）:
8  9  3  2  1  5  2  4  4  9  4  5  

内存中的存储顺序（行优先）:
8  2  2  9  9  1  4  4  3  5  4  5

2、Eigen中的存储顺序

可以通过为 Matrix 或 Array 指定 Options 模板参数来设置矩阵或二维数组的存储顺序。正如 Matrix 类所解释的，Matrix 类模板有六个模板参数，其中三个是必选的（Scalar，RowsAtCompileTime 和 ColsAtCompileTime），三个是可选的（Options，MaxRowsAtCompileTime 和 MaxColsAtCompileTime）。如果 Options 参数设置为 RowMajor，则矩阵或数组按行优先存储；如果设置为 ColMajor，则按列优先存储。上述Eigen程序中的代码就是通过这种机制来指定存储顺序的。

如果未指定存储顺序，则Eigen默认按列优先存储条目。即使使用方便的类型定义（如 Matrix3f，ArrayXXd 等）也是如此。

使用一种存储顺序的矩阵和数组可以赋值给使用另一种存储顺序的矩阵和数组，如上述程序中 Arowmajor 使用 Acolmajor 进行初始化时那样。Eigen会自动重新排列条目。更一般地说，我们可以在表达式中混合使用行优先和列优先的矩阵。

3、选择哪种存储顺序？

那么，在你的程序中应该使用哪种存储顺序呢？对此问题没有简单的答案；这取决于你的应用程序。以下是一些需要记住的要点：

你的用户可能期望你使用特定的存储顺序。或者，你可能会使用除Eigen之外的其他库，而这些库可能需要特定的存储顺序。在这些情况下，在整个程序中使用这种存储顺序可能是最简单和最快的。
按行遍历矩阵的算法在矩阵按行优先存储时会更快，因为数据局部性更好。同样，按列遍历对于列优先矩阵更快。可能值得做一些实验，以找出对于你的特定应用程序哪种方式更快。
Eigen默认使用列优先存储。因此，Eigen库的大部分开发和测试都是在列优先矩阵上进行的。这意味着，尽管我们致力于透明地支持列优先和行优先存储顺序，但Eigen库可能在使用列优先矩阵时效果最好。

标签：Eigen,矩阵,C++,matrix,3.4,数组,include,mat
From： https://blog.csdn.net/m0_70254540/article/details/139813333

C++矩阵库：Eigen 3.4.90 中文使用文档 （一）