Eigen 库学习指南1-稠密矩阵和数组的运算(一)

2023年01月18日 · 8985字 · 18分钟

翻译: Qi Deng

简介 🔗

Eigen 的介绍 🔗

Eigen 是一个高级 C++ 库，包含用于线性代数、矩阵和向量运算、几何变换、数值求解器和相关算法的模板头文件。 Eigen 是自 3.1.1 版起根据 Mozilla Public License 2.0 许可的开源软件。早期版本是根据 GNU 宽松通用公共许可证获得许可的。 [2] 1.0 版于 2006 年 12 月发布

使用了 expression templates metaprogramming 来实现，具体来说就是: 这意味着它在编译时构建表达式树并生成自定义代码来评估这些树。使用表达式模板和浮点运算的成本模型，库执行自己的循环展开 loop unrolling 和矢量化 vectorization。

Eigen itself can provide BLAS and a subset of LAPACK interfaces

Eigen 的特性 🔗

Eigen is versatile(通用途的)
- 它支持所有矩阵大小，从小的固定大小矩阵到任意大的稠密矩阵，甚至是稀疏矩阵。
- 它支持所有标准数字类型，包括 std::complex、整数，并且可以轻松扩展为自定义数字类型。
- 它支持各种矩阵分解和几何特征。
- 其不受支持的模块生态系统提供了许多专门的功能，例如非线性优化、矩阵函数、多项式求解器、FFT 等等。
Eigen is fast
- 表达式模板允许智能地删除临时对象并在适当时启用延迟评估。
- 对 SSE 2/3/4、AVX、AVX2、FMA、AVX512、ARM NEON（32 位和 64 位）、PowerPC AltiVec/VSX（32 位和 64 位）、ZVector (s390x/ zEC13) SIMD 指令集，并且从 3.4 MIPS MSA 开始可以优雅地回退到非矢量化代码。
- 固定大小的矩阵已完全优化：避免了动态内存分配，并在有意义时展开循环。
- 对于大型矩阵，特别注意了缓存友好性。
Eigen is reliable
- 算法经过精心挑选以确保可靠性。可靠性权衡被清楚地记录下来，并且可以使用极其安全的分解。
- Eigen 通过其自己的测试套件（超过 500 个可执行文件）、标准 BLAS 测试套件和部分 LAPACK 测试套件进行了全面测试。
Eigen is elegant
- 由于表达式模板，API 非常简洁且富有表现力，同时让 C++ 程序员感觉很自然。
- 在 Eigen 之上实现算法感觉就像只是复制伪代码。
Eigen has good compiler support
- Eigen 拥有良好的编译器支持，因为我们针对许多编译器运行我们的测试套件以保证可靠性并解决任何编译器错误。 Eigen 版本 3.4 是标准的 C++03，并保持合理的编译时间。 3.4 之后的版本将是 C++14

The matrix class 🔗

默认都实在 C++14 之后

前三个模板参数 🔗

在 Eigen 中，Matrix Class 是基本类， Vector 只不过是 column 或者 row 数目为 1 的特殊矩阵。

Matrix 类的参数共有六个，一般后三个保持默认就可以，前三个需要自身定义:

Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>

Scalar 表示系数的类型，比如 float, double 或者 int64_t 当然也可以是定点数
RowsAtCompileTime 和 ColsAtCompileTime 是在编译时已知的矩阵的行数和列数。如果不知道，请使用下面的动态矩阵类。

Eigen 中提供了很多方便的 typedef 来涵盖常见的种类，比如:

typedef Matrix<float, 4, 4> Matrix4f;
typedef Matrix<int, 4, 4> Matrix4i;

Vectors 向量 🔗

在 Eigen 中，向量只是矩阵的一种特殊情况，只有 1 行或 1 列。他们有 1 列的情况是最常见的；这样的向量称为列向量(column-vectors)，通常缩写为向量。在另一种情况下，它们只有一行，它们被称为行向量(row-vectors)。

Eigen 中提供了很多方便的 typedef 来涵盖常见的向量种类

typedef Matrix<float, 3, 1> Vector3f;
`typedef Matrix<int, 1, 2> RowVector2i;

特殊值 Dynamic 🔗

Eigen 并不局限于编译时维度已知的矩阵。 RowsAtCompileTime 和 ColsAtCompileTime 模板参数可以采用特殊值 Dynamic，表示编译时大小未知，因此必须作为运行时变量处理。编译时时知道大小称之为 fixed size, 编译时不知道大小的称之为 dynamic size 。

Eigen 中提供呢许多方便的 typedef 来涵盖常见的动态矩阵或向量类型:

typedef Matrix<double, Dynamic, Dynamic> MatrixXd;
typedef Matrix<int, Dynamic, 1> VectorXi;
固定行，动态列: Matrix<float, 3, Dynamic>

构造函数 && 给矩阵赋值 🔗

默认构造函数
- Matrix3f a: 3x3 的系数，系数分配了空间，但没有初始化
- MatrixXf b: 0x0 的系数，系数还没有分配
带大小的构造函数，矩阵顺序 行数，列数，向量传递大小。分配系数数组空间，但是不会初始化。
- MatrixXf a(10,15);
- VectorXf b(30);
- 在固定大小的数组上使用这个方法也是可以的, Matrix3f a(3,3);
可以通过传递任意数量的系数来初始化任意大小的固定大小的列或行向量
- Vector2i a(1, 2); :包含元素 {1, 2} 的列向量
- Matrix<int, 5, 1> b {1, 2, 3, 4, 5}; : 包含元素 {1, 2, 3, 4, 5} 的行向量
- Matrix<int, 1, 5> c = {1, 2, 3, 4, 5};: 包含元素 {1, 2, 3, 4, 5} 的列向量
具有固定大小或运行时大小的矩阵和向量的一般情况下，系数必须按行分组并作为初始化列表(initializer list)的初始化列表传递

MatrixXi a { // 构造一个 2x2 矩阵
	{1, 2}, // 第一行
	{3, 4} // 第二行
};
Matrix<double, 2, 3> b {
      {2, 3, 4},
      {5, 6, 7},
};

对于列或行向量，隐式转置是允许的。这意味着可以从单行初始化列向量：

VectorXd {{1.5, 2.5, 3.5}}; // 具有 3 个系数的列向量
RowVectorXd b {{1.0, 2.0, 3.0, 4.0}}; // 一个有 4 个系数的行向量

可以使用所谓的逗号初始化语法方便地设置矩阵和矢量系数。

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

访问值 🔗

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;
Eigen::VectorXd v(2);
v(0) = 4;
v(1) = v(0) - 1;
std::cout << "Here is the vector v:\n" << v << std::endl;

请注意，语法 m(index) 不限于向量，它也可用于一般矩阵，这意味着在系数数组中进行基于索引的访问。然而，这取决于矩阵的存储顺序。所有特征矩阵默认为列优先存储顺序，但这可以更改为行优先，请参阅存储顺序。operator[] 也被重载用于向量中基于索引的访问，但请记住 C++ 不允许 operator[] 接受多个参数。我们将 operator[] 限制为向量，因为 C++ 语言中的尴尬会使 matrix[i,j] 编译为与 matrix[j] 相同的东西！

问一个问题，为什么这里不重载 [] 运算符号呢？

Resize 重新挑战大小 🔗

可以通过 rows()、cols() 和 size() 检索矩阵的当前大小。这些方法分别返回行数、列数和系数数。调整动态大小矩阵的大小是通过 resize() 方法完成的。

Eigen::MatrixXd m(2,5);
m.resize(4,3);

Eigen::VectorXd v(2);
v.resize(5);

如果实际矩阵大小不变，则 resize() 方法是空操作；否则它是破坏性的：系数的值可能会改变。如果您想要不改变系数的 resize() 的保守变体，请使用 conservativeResize()，请参阅此页面了解更多详细信息。

尝试让矩阵的大小发生改变，看看会发生什么？

为了 API 统一，所有这些方法在固定大小的矩阵上仍然可用。当然，您实际上无法调整固定大小的矩阵。尝试将固定大小更改为实际不同的值将触发断言失败；但是下面的代码是合法的：

#include <iostream>
#include <Eigen/Dense>
 
int main()
{
  Eigen::Matrix4d m;
  m.resize(4,4); // no operation
  std::cout << "The matrix m is of size "
            << m.rows() << "x" << m.cols() << std::endl;
}

赋值 🔗

赋值是使用 operator= 将矩阵复制到另一个矩阵的操作。 Eigen 会自动调整左侧矩阵的大小，使其与右侧矩阵的大小相匹配。例如：

MatrixXf a(2,2);
std::cout << "a is of size " << a.rows() << "x" << a.cols() << std::endl;
MatrixXf b(3,3);
a = b;
std::cout << "a is now of size " << a.rows() << "x" << a.cols() << std::endl;

当然，如果左侧是固定大小，则不允许调整它的大小。

固定大小与动态大小 🔗

什么时候应该使用固定大小（例如 Matrix4f），什么时候应该使用动态大小（例如 MatrixXf）？简单的答案是：尽可能对非常小的尺寸使用固定尺寸，对较大尺寸或必须使用动态尺寸。对于小尺寸，尤其是小于（大约）16 的尺寸，使用固定尺寸对性能非常有利，因为它允许 Eigen 避免动态内存分配和展开循环。在内部，固定大小的特征矩阵只是一个普通数组，即做

Matrix4f mymatrix;

真的等于只是做

float mymatrix[16];

所以这真的有零运行成本。相比之下，动态大小矩阵的数组总是分配在堆上，所以这样做：

MatrixXf mymatrix(rows,columns);

相当于做

float *mymatrix = new float[行*列];

除此之外，MatrixXf 对象将其行数和列数存储为成员变量。

当然，使用固定大小的限制是，只有在编译时知道大小时，这才有可能。此外，对于足够大的尺寸，比如大于（大约）32 的尺寸，使用固定尺寸的性能优势变得可以忽略不计。更糟糕的是，尝试在函数内部使用固定大小创建一个非常大的矩阵可能会导致堆栈溢出，因为 Eigen 会尝试将数组自动分配为局部变量，而这通常是在堆栈上完成的。最后，根据情况，当使用动态大小时，Eigen 也可以更积极地尝试矢量化（使用 SIMD 指令），请参阅矢量化。

可选模板参数 🔗

我们在本页开头提到 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 很有用，即使在编译时不知道矩阵的确切大小，但在编译时已知固定的上限。您可能想要这样做的最大原因是避免动态内存分配。例如，以下矩阵类型使用 12 个浮点数的普通数组，没有动态内存分配：

一些方便的 convenience typedefs 🔗

Eigen 定义了以下 Matrix 类型定义：

MatrixNt 表示 Matrix<type, N, N>。例如，MatrixXi 代表 Matrix<int, Dynamic, Dynamic>。
MatrixXNt 用于 Matrix<type, Dynamic, N>。例如，MatrixX3i 用于 Matrix<int, Dynamic, 3>。
MatrixNXt 用于 Matrix<type, N, Dynamic>。例如，Matrix4Xd 表示 Matrix<d, 4, Dynamic>。
Matrix<type, N, 1> 的 VectorNt。例如，Matrix<float, 2, 1> 的 Vector2f。
Matrix<type, 1, N> 的 RowVectorNt。例如，Matrix<double, 1, 3> 的 RowVector3d。

其中:

N 可以是 2、3、4 或 X（表示动态）中的任何一个。
t 可以是 i（表示 int）、f（表示 float）、d（表示 double）、cf（表示 complex<float>）或 cd（表示 complex<double>）中的任何一个。 typedef 仅为这五种类型定义的事实并不意味着它们是唯一受支持的标量类型。例如，支持所有标准整数类型，请参阅标量类型。

矩阵和向量的运算 🔗

Eigen 通过重载常见的 C++ 算术运算符（如 +、-、*）或通过特殊方法（如 dot()、cross() 等）提供矩阵/向量算术运算。对于 Matrix 类（矩阵和向量），运算符仅重载以支持线性代数运算。例如matrix1 * matrix2表示矩阵-矩阵乘积，向量+标量就是不允许的。如果您想执行各种数组运算，而不是线性代数，请参阅下一节。

加减 🔗

左右两边的矩阵需要有相同的行数和列数:

二元运算符 + 如 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 << "Doing a += b;" << std::endl;
  a += b;
  std::cout << "Now 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;
}

标量乘法和除法 🔗

二元运算符 * 如 matrix*scalar
二元运算符 * 如 scalar*scalar
二元运算符 / 如在 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 << "Doing v *= 2;" << std::endl;
  v *= 2;
  std::cout << "Now v =\n" << v << std::endl;
}

这是我们在此页面上解释的高级主题，但现在仅提及它很有用。在 Eigen 中，诸如 operator+ 之类的算术运算符本身不执行任何计算，它们只是返回一个描述要执行的计算的“表达式对象”。实际计算发生在稍后，当整个表达式被计算时，通常在 operator= 中。虽然这听起来很难，但任何现代优化编译器都能够优化抽象，结果是完美优化的代码。例如，当您执行以下操作时：

VectorXf a(50), b(50), c(50), d(50);
...
a = 3*b + 4*c + 5*d；

Eigen 将其编译为一个 for 循环，因此数组只被遍历一次。简化（例如忽略 SIMD 优化），此循环如下所示：

for(int i = 0; i < 50; ++i)
a[i] = 3*b[i] + 4*c[i] + 5*d[i];

Transposition and conjugation 转置和共轭 🔗

transpose转置 $a^T$
conjugate共轭 $\bar{a}$
adjoin伴随（即共轭转置）$a^{\star}$

MatrixXcf a = MatrixXcf::Random(2,2);
cout << "Here is the matrix a\n" << a << endl;
cout << "Here is the matrix a^T\n" << a.transpose() << endl;
cout << "Here is the conjugate of a\n" << a.conjugate() << endl;
cout << "Here is the matrix a^*\n" << a.adjoint() << endl;

对于实数矩阵，conjugate() 是空操作，因此 adjoint() 等价于 transpose()。

对于基本的算术运算符，transpose() 和 adjoint() 只是返回一个代理对象，而不进行实际的转置。如果执行 b = a.transpose()，则在将结果写入 b 的同时计算转置。然而，这里有一个复杂的问题。如果您执行 a = a.transpose()，则 Eigen 在转置计算完成之前开始将结果写入 a。因此，指令 a = a.transpose() 不会像预期的那样用它的转置替换 a。

Matrix2i a; a << 1, 2, 3, 4;
cout << "Here is the matrix a:\n" << a << endl;
 
a = a.transpose(); // !!! do NOT do this !!!
cout << "and the result of the aliasing effect:\n" << a << endl;

输出是什么?

这就是所谓的别名问题。在“调试模式”下，即当断言未被禁用时，会自动检测到此类常见陷阱。

对于就地转置，例如在 a = a.transpose() 中，只需使用 transposeInPlace() 函数：

MatrixXf a(2,3); a << 1, 2, 3, 4, 5, 6;
cout << "Here is the initial matrix a:\n" << a << endl;

a.transposeInPlace();
cout << "and after being transposed:\n" << a << endl;

还有用于复杂矩阵的 adjointInPlace() 函数。

矩阵-矩阵和矩阵-向量乘法 🔗

矩阵-矩阵乘法再次使用 operator* 完成。由于向量是矩阵的特例，它们在那里也被隐式处理，所以矩阵-向量乘积实际上只是矩阵-矩阵乘积的特例，向量-向量外积也是如此。因此，所有这些情况仅由两个操作员处理：

二元运算符 *\ 如 a*b
复合运算符 *= 与 a*=b 相同（这在右边相乘：a*=b 等同于 a = a*b）

#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 << "Here is mat*mat:\n" << mat*mat << std::endl;
  std::cout << "Here is mat*u:\n" << mat*u << std::endl;
  std::cout << "Here is u^T*mat:\n" << u.transpose()*mat << std::endl;
  std::cout << "Here is u^T*v:\n" << u.transpose()*v << std::endl;
  std::cout << "Here is u*v^T:\n" << u*v.transpose() << std::endl;
  std::cout << "Let's multiply mat by itself" << std::endl;
  mat = mat*mat;
  std::cout << "Now mat is mat:\n" << mat << std::endl;
}

注意：如果您阅读了上面关于表达式模板的段落并且担心执行 m=m*m 可能会导致别名问题，现在请放心：Eigen treats matrix multiplication as a special case and takes care to introducing a temporary here，所以它将 m=m*m 编译为：

tmp = m*m;
m = tmp;

如果您知道您的矩阵乘积可以安全地计算到目标矩阵中而没有别名问题，那么您可以使用 noalias() 函数来避免临时的变量，例如：

c.noalias() += a * b;

有关此主题的更多详细信息，请参阅别名页面。注意：对于担心性能的 BLAS 用户，表达式如 c.noalias() -= 2 * a.adjoint() * b;已完全优化并触发单个类似 gemm 的函数调用。

点积和叉积 🔗

对于点积和叉积，您需要 dot() 和 cross() 方法。当然，点积也可以像 u.adjoint()*v 一样得到一个1x1的矩阵。

#include <iostream>
#include <Eigen/Dense>
 
int main()
{
  Eigen::Vector3d v(1,2,3);
  Eigen::Vector3d w(0,1,2);
 
  std::cout << "Dot product: " << v.dot(w) << std::endl;
  double dp = v.adjoint()*w; // automatic conversion of the inner product to a scalar
  std::cout << "Dot product via a matrix product: " << dp << std::endl;
  std::cout << "Cross product:\n" << v.cross(w) << std::endl;
}

请记住，叉积仅适用于大小为 3 的向量。点积适用于任何大小的向量。使用复数时，Eigen 的点积在第一个变量中是共轭线性的，在第二个变量中是线性的。

基本算术归约运算 🔗

Eigen 还提供了一些缩减操作来将给定的矩阵或向量缩减为单个值，例如和（由 sum() 计算）、乘积 (prod()) 或最大值 (maxCoeff()) 和最小值 (minCoeff() ) 的所有系数。

#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;
}

函数 trace() 返回的矩阵迹是对角线系数的总和，也可以使用 a.diagonal().sum() 高效计算，我们将在后面看到。

还存在 minCoeff 和 maxCoeff 函数的变体，它们通过参数返回相应系数的坐标：

  Matrix3f m = Matrix3f::Random();
  std::ptrdiff_t i, j;
  float minOfM = m.minCoeff(&i,&j);
  cout << "Here is the matrix m:\n" << m << endl;
  cout << "Its minimum coefficient (" << minOfM 
       << ") is at position (" << i << "," << j << ")\n\n";
 
  RowVector4i v = RowVector4i::Random();
  int maxOfV = v.maxCoeff(&i);
  cout << "Here is the vector v: " << v << endl;
  cout << "Its maximum coefficient (" << maxOfV 
       << ") is at position " << i << endl;

操作的有效性 🔗

Eigen 检查您执行的操作的有效性。如果可能，它会在编译时检查它们，从而产生编译错误。这些错误消息可能又长又难看，但 Eigen 将重要消息写在 UPPERCASE_LETTERS_SO_IT_STANDS_OUT 中。例如

Matrix3f m;
Vector4f v;
v = m*v;      // Compile-time error: YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES

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

MatrixXf m(3,3);
VectorXf v(4);
v = m * v; // Run-time assertion failure here: "invalid matrix product"

可以参考这个页面。

Array 类和系数运算 🔗

Array 类提供通用数组，而不是用于线性代数的 Matrix 类。此外，Array 类提供了一种执行按系数运算的简单方法，这可能没有线性代数意义，例如将常数添加到数组中的每个系数或按系数将两个数组相乘。

Array types 🔗

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

Array<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>

最后三个模板参数是可选的。由于这与 Matrix 完全相同，我们这里不再解释，仅参考 Matrix 类。

Eigen 还为一些常见情况提供类型定义，其方式类似于矩阵类型定义，但有一些细微差别，因为“数组”一词用于一维和二维数组。我们采用 ArrayNt 形式的类型定义代表一维数组的约定，其中 N 和 t 是大小和标量类型，如本页解释的矩阵类型定义。对于二维数组，我们使用 ArrayNNt 形式的 typedef。下表显示了一些示例：

Type	Typedef
`Array<float,Dynamic,1>`	`ArrayXf`
`Array<float,3,1>`	`Array3f`
`Array<double,Dynamic,Dynamic>`	`ArrayXXd`
`Array<double,3,3>`	`Array33d`

访问 array 的值 🔗

括号运算符被重载以提供对数组系数的读写访问，就像矩阵一样。此外，<< 运算符可用于初始化数组（通过逗号初始化程序）或打印它们。

#include <Eigen/Dense>
#include <iostream>
 
int main()
{
  Eigen::ArrayXXf  m(2,2);
  
  // assign some values coefficient by coefficient
  m(0,0) = 1.0; m(0,1) = 2.0;
  m(1,0) = 3.0; m(1,1) = m(0,1) + m(1,0);
  
  // print values to standard output
  std::cout << m << std::endl << std::endl;
 
  // using the comma-initializer is also allowed
  m << 1.0,2.0,
       3.0,4.0;
     
  // print values to standard output
  std::cout << m << std::endl;
}

array 的加减 🔗

两个数组的加减法与矩阵相同。如果两个数组的大小相同，并且加法或减法是按系数进行的，则该操作有效。数组还支持 array + scalar 形式的表达式，它向数组中的每个系数添加一个标量。这提供了 Matrix 对象不能直接使用的功能。

#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;
       
  // Adding two arrays
  std::cout << "a + b = " << std::endl << a + b << std::endl << std::endl;
 
  // Subtracting a scalar from an array
  std::cout << "a - 2 = " << std::endl << a - 2 << std::endl;
}

array 的乘法 🔗

首先，当然你可以将一个数组乘以一个标量，这与矩阵的工作方式相同。数组与矩阵根本不同的地方在于将两者相乘。矩阵将乘法解释为矩阵乘积，而数组将乘法解释为系数乘积。因此，两个数组可以相乘当且仅当它们具有相同的维度。

#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;
       
  // Adding two arrays
  std::cout << "a + b = " << std::endl << a + b << std::endl << std::endl;
 
  // Subtracting a scalar from an array
  std::cout << "a - 2 = " << std::endl << a - 2 << std::endl;
}

其他系数运算 🔗

除了上述的加法、减法和乘法运算符之外，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;
}

在数组和矩阵表达式之间转换 🔗

什么时候应该使用 Matrix 类的对象，什么时候应该使用 Array 类的对象？您不能对数组应用 Matrix 运算，也不能对矩阵应用 Array 运算。因此，如果您需要进行线性代数运算，例如矩阵乘法，那么您应该使用矩阵；如果你需要做系数运算，那么你应该使用数组。然而，有时并没有那么简单，而是需要同时使用 Matrix 和 Array 操作。在这种情况下，您需要将矩阵转换为数组或相反。这允许访问所有操作，而不管选择将对象声明为数组还是矩阵。

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

Eigen 禁止在表达式中混合矩阵和数组。例如，您不能直接添加矩阵和数组； + 运算符的操作数要么都是矩阵，要么都是数组。但是，使用 .array() 和 .matrix() 很容易从一个转换为另一个。此规则的例外是赋值运算符：允许将矩阵表达式分配给数组变量，或将数组表达式分配给矩阵变量。

以下示例显示如何通过使用 .array() 方法对 Matrix 对象使用数组操作。例如，语句 result = m.array() * n.array() 接受两个矩阵 m 和 n，将它们都转换为数组，使用它们进行系数相乘并将结果赋给矩阵变量 result（这是合法的，因为 Eigen 允许将数组表达式分配给矩阵变量）。

事实上，这种用例非常普遍，以至于 Eigen 为矩阵提供了一个 const .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";
}

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

这是一个更高级的示例。表达式 (m.array() + 4).matrix() * m 将矩阵 m 中的每个系数加 4，然后计算结果与 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";
}

块操作 🔗

本页解释了块操作的要点。块是矩阵或数组的矩形部分。块表达式既可以用作右值也可以用作左值。与 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 << "Block in the middle" << endl;
  cout << m.block<2,2>(1,1) << endl << endl;
  for (int i = 1; i <= 3; ++i)
  {
    cout << "Block of size " << i << "x" << i << endl;
    cout << m.block(0,0,i,i) << endl << endl;
  }
}

在上面的示例中，.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 << "Here is the array a:\n" << a << "\n\n";
  a.block<2,2>(1,1) = m;
  std::cout << "Here is now a with m copied into its central 2x2 block:\n" << a << "\n\n";
  a.block(0,0,2,3) = a.block(2,1,2,3);
  std::cout << "Here is now a with bottom-right 2x3 block copied into top-left 2x3 block:\n" << a << "\n\n";
}

虽然 .block() 方法可用于任何块操作，但还有其他方法用于特殊情况，提供更专业的 API 和/或更好的性能。关于性能这个话题，重要的是你在编译时给 Eigen 尽可能多的信息。例如，如果您的块是矩阵中的一整列，使用下面描述的专门的 .col() 函数可以让 Eigen 知道这一点，这可以为其提供优化机会。

列和行 🔗

单独的列和行是块的特例。 Eigen 提供了轻松解决它们的方法：.col() 和 .row()。

块操作	方法
$i^{th}$ 行	`matrix.row(i)`
$j^{th}$ 列	`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 << "Here is the matrix m:" << endl << m << endl;
  cout << "2nd Row: " << m.row(1) << endl;
  m.col(2) += 3 * m.col(0);
  cout << "After adding 3 times the first column into the third column, the matrix m is:\n";
  cout << m << endl;
}

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

角相关操作 🔗

Eigen 还为与矩阵或数组的一个角或边之一齐平的块提供特殊方法。例如，.topLeftCorner() 可用于引用矩阵左上角的块。下表列出了所有可能的情况:

Block operation	Version constructing a dynamic-size block expression	Version constructing a fixed-size block expression
Top-left p by q block *	matrix.topLeftCorner(p,q);	matrix.topLeftCorner<p,q>();
Bottom-left p by q block *	matrix.bottomLeftCorner(p,q);	matrix.bottomLeftCorner<p,q>();
Top-right p by q block *	matrix.topRightCorner(p,q);	matrix.topRightCorner<p,q>();
Bottom-right p by q block *	matrix.bottomRightCorner(p,q);	matrix.bottomRightCorner<p,q>();
Block containing the first q rows *	matrix.topRows(q);	matrix.topRows();
Block containing the last q rows *	matrix.bottomRows(q);	matrix.bottomRows();
Block containing the first p columns *	matrix.leftCols(p);	matrix.leftCols ();
Block containing the last q columns *	matrix.rightCols(q);	matrix.rightCols();
Block containing the q columns starting from i *	matrix.middleCols(i,q);	matrix.middleCols(i);
Block containing the q rows starting from i *	matrix.middleRows(i,q);	matrix.middleRows(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 << "After assignment, m = " << endl << m << endl;
}

向量的块操作 🔗

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

Block operation	Version constructing a dynamic-size block expression	Version constructing a fixed-size block expression
Block containing the first n elements *	vector.head(n);	vector.head();
Block containing the last n elements *	vector.tail(n);	vector.tail();
Block containing n elements, starting at position i *	vector.segment(i,n);	vector.segment(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 << "after 'v.segment(1,4) *= 2', v =" << endl << v << endl;
}