C++之 Eigen-3.4.0 全方位教程：Chapter02-矩阵篇

矩阵类（The Matrix class）

在Eigen 中，所有矩阵和向量都是Matrix模板类的对象。向量只是矩阵的一种特殊情况，具有 1 行或 1 列。

Matrix的前三个模板参数

该矩阵类需要六个模板参数，但是现在了解前三个参数已经足够了。

剩下的三个参数有默认值，现在我们将保持不变，我们将在下面讨论。

Matrix的三个强制性模板参数是：

Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>

Scalar是标量类型，即系数的类型。也就是说，如果您想要一个浮点矩阵，请在此处选择float。有关所有支持的标量类型的列表以及如何扩展对新类型的支持，请参阅标量类型。
RowsAtCompileTime和ColsAtCompileTime是编译时已知的矩阵的行数和列数（如果在编译时不知道该数字，请参见下文了解该怎么做）。

我们提供了很多方便的 typedef 来涵盖通常的情况。例如，Matrix4f是一个 4x4 的浮点矩阵。这是Eigen定义的方式：

typedef Matrix<float, 4, 4> Matrix4f;

Eigen里面用到了很多的typedef简化名称长度，例如：

typedef Matrix<float, 3, 1> Vector3f;
typedef Matrix<int, 1, 2> RowVector2i;
typedef Matrix<double, Dynamic, Dynamic> MatrixXd;
typedef Matrix<int, Dynamic, 1> VectorXi;typedef Array<float,Dynamic,Dynamic>       ArrayXXf
typedef Array<double,Dynamic,1>            ArrayXd
typedef Array<int,1,Dynamic>               RowArrayXi
typedef Array<float,3,3>                   Array33f
typedef Array<float,4,1>                   Array4f
······

Tips: 虽然Eigen官方已经提供了许多命名模板，但是这是对于X <= 4 的情况比较试用，一般处理大数据时，为了尽可能的多实用固定矩阵，应该熟练掌握利用原始模板创建对象

Vector

如上所述，在Eigen 中，向量只是矩阵的一种特殊情况，具有 1 行或 1 列。他们有 1 列的情况是最常见的；这样的向量称为列向量，通常简称为向量。在它们有 1 行的另一种情况下，它们被称为行向量。

例如，方便的 typedef Vector3f是 3 个浮点数的（列）向量。它由Eigen定义如下：

typedef Matrix<float, 3, 1> Vector3f;

我们还为行向量提供方便的 typedef，例如：

typedef Matrix<int, 1, 2> RowVector2i;

The special value Dynamic

当然，Eigen不限于其维度在编译时已知的矩阵。在RowsAtCompileTime和ColsAtCompileTime模板参数可以采取特殊值Dynamic这表明大小在编译时是未知的，所以必须作为运行时变量来处理。在Eigen术语中，这样的大小被称为动态大小；而在编译时已知的大小称为固定大小。

例如，方便的 typedef MatrixXd 表示具有动态大小的双精度矩阵，定义如下：

typedef Matrix<double, Dynamic, Dynamic> MatrixXd;

同样，我们定义了一个动态大小的 typedef VectorXi如下：

typedef Matrix<int, Dynamic, 1> VectorXi;

还可以固定某行或某列的大小：

Matrix<float, 3, Dynamic>;

构造函数

默认构造函数始终可用，从不执行任何动态内存分配，也从不初始化矩阵系数。你可以做：

Matrix3f a; //a 是一个 3×3 矩阵，带有未初始化系数的普通 float[9] 数组，

MatrixXf b; //b 是一个动态大小的矩阵，其大小当前为 0×0，并且其系数数组尚未分配。

也可以使用获取尺寸的构造函数。

对于矩阵，总是首先传递行数。对于向量，只需传递向量大小。他们分配给定大小的系数数组，但不初始化系数本身：

MatrixXf a(10,15); //a 是一个 10x15 动态大小的矩阵，具有已分配但当前未初始化的系数。

VectorXf b(30); //b 是大小为 30 的动态大小向量，具有已分配但当前未初始化的系数。

为了提供跨固定大小和动态大小矩阵的统一 API，在固定大小矩阵上使用这些构造函数是合法的，即使在这种情况下传递大小是无用的。所以这是合法的：

Matrix3f a(3,3); //虽然没有意义，但是是合法的

矩阵和向量也可以从系数列表中初始化。在 C++11 之前，此功能仅限于固定大小的小列或最大大小为 4 的向量：

Vector4d c(5.0, 6.0, 7.0, 8.0); //C++11之前只能用于固定尺寸的矩阵和向量，且尺寸不大于4

如果启用 C++11，则可以通过传递任意数量的系数来初始化任意大小的固定大小的列或行向量：

Vector2i a(1, 2);                      // A column vector containing the elements {1, 2}
Matrix<int, 5, 1> b {1, 2, 3, 4, 5};   // A row-vector containing the elements {1, 2, 3, 4, 5}
Matrix<int, 1, 5> c = {1, 2, 3, 4, 5}; // A column vector containing the elements {1, 2, 3, 4, 5}

无论是固定还是动态矩阵，系数必须按行分组：

MatrixXi a {      // construct a 2x2 matrix{1, 2},     // first row{3, 4}      // second row
};
Matrix<double, 2, 3> b {{2, 3, 4},{5, 6, 7},
};

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

切记，只有向量才可以这样

VectorXd a {{1.5, 2.5, 3.5}};             // A column-vector with 3 coefficients
RowVectorXd b {{1.0, 2.0, 3.0, 4.0}};     // A row-vector with 4 coefficients

系数存取器

Eigen重载了圆括号()操作符，因此采用 (i,j) 取索引访问系数，对于向量只需要传递一个索引(i)即可，需要注意的是下标索引均是先行后列。

#include <iostream>
#include <Eigen/Dense>using namespace Eigen;int main()
{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;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] 相同的东西！

逗号初始化

矩阵和向量系数可以使用所谓的逗号初始化语法方便地设置，后续会介绍更高级的初始化方式。

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

调整大小

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

注意：

只有动态矩阵才能调整大小（为了 API 的一致性，所有这些方法在固定大小的矩阵上仍然可用。当然，您实际上无法调整固定大小的矩阵的大小。尝试将固定大小更改为实际不同的值将触发断言失败。）
resize() 如果没有改变矩阵大小，则不会对其进行操作；否则就会重置矩阵；如果您想要不改变系数的 resize() 的保守变体，请使用conservativeResize()

例1：

void Eigen_MatrixClass_Resizeing_001()
{MatrixXd m = Matrix2d::Ones();cout << m << endl;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;VectorXd v(2);v.resize(5);std::cout << "The vector v is of size " << v.size() << std::endl;std::cout << "As a matrix, v is of size "<< v.rows() << "x" << v.cols() << std::endl;
}

输出

1 1
1 1
The matrix m is of size 4x3
It has 12 coefficients
The vector v is of size 5
As a matrix, v is of size 5x1

例2：

void Eigen_MatrixClass_Resizeing_002()
{MatrixXd m = Matrix2d::Ones();cout << "调整前:" << endl << m << endl;m.resize(2, 2);cout << "不改变尺寸:" << endl << m << endl;m.conservativeResize(4, 3);cout << "保留调整后:" << endl << m << endl;RowVectorXd v(2);v << 1, 2;v.conservativeResize(5);cout << "保留调整后:" << endl << v << endl;
}

输出

调整前:
1 1
1 1
不改变尺寸:
1 1
1 1
保留调整后:1            1 -6.27744e+661            1 -6.27744e+66
-6.27744e+66 -6.27744e+66 -6.27744e+66
-6.27744e+66 -6.27744e+66 -6.27744e+66
保留调整后:1            2 -6.27744e+66 -6.27744e+66 -6.27744e+66

分配和调整大小

Eigen 采用 operator= 进行赋值，自动调整左侧矩阵的大小，使其与右侧矩阵的大小匹配。

结合前面的描述，对于 动态矩阵 = 随意矩阵 这样赋值是可以随意进行的，但是对于固定矩阵，只有 固定矩阵 = 同尺寸的矩阵 才可行，因为在调整大小时，resize()并不会对其进行操作，或者理解为 固定矩阵不能改变大小。

固定尺寸与动态尺寸

什么时候应该使用固定尺寸（例如Matrix4f），什么时候应该使用动态尺寸（例如MatrixXf）？

简单的答案是：尽可能对非常小的尺寸使用固定尺寸，对较大尺寸或必须使用的尺寸使用动态尺寸。对于小尺寸，尤其是小于（大约）16 的尺寸，使用固定尺寸对性能非常有益，因为它允许Eigen避免动态内存分配并展开循环。在内部，固定大小的特征矩阵只是一个普通数组。

固定尺寸

Matrix4f mymatrix; 等价于 float mymatrix[16];

零成本的运行，所以非常快
动态尺寸

MatrixXf mymatrix(rows,columns); 等价于 float *mymatrix = new float[rows*columns];

相比之下，动态大小矩阵的数组总是分配在堆上。

除此之外，MatrixXf 对象将其行数和列数存储为成员变量。
较大尺寸时，请用动态矩阵

对于足够大的尺寸，比如大于（大约）32 的尺寸，使用固定尺寸的性能优势变得可以忽略不计。

更糟糕的是，尝试在函数内部使用固定大小创建一个非常大的矩阵可能会导致栈溢出，因为Eigen会尝试将数组自动分配为局部变量，而这通常是在栈上完成的。

最后，根据情况，当使用动态大小时，Eigen也可以更积极地尝试矢量化（使用 SIMD 指令），请参阅Vectorization。

可选模板参数

我们在本章开头提到Matrix类需要六个模板参数，但到目前为止我们只讨论了前三个。其余三个参数是可选的。以下是模板参数的完整列表：

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

Options 是一个位域。在这里，我们只讨论一点：RowMajoror 1 . 它指定这种类型的矩阵使用行优先存储顺序；默认情况下，存储顺序为列优先，ColMajor or 0。请参阅有关StorageOrders的页面。例如，这种类型表示行主 3x3 矩阵：矩阵<float, 3, 3, RowMajor>
MaxRowsAtCompileTime and MaxColsAtCompileTime 预分配动态矩阵的行列上限，避免动态内存分配。

void Eigen_MatrixClass_TemplateParameters_001()
{Matrix<int, Dynamic, Dynamic, 0, 4, 4> m;cout << m.size() << endl;m.resize(2, 2);cout << m.size() << endl;m.resize(3, 3);cout << m.size() << endl;//m.resize(5, 3);  //error:(i,j) 中的 i and j 均小于 原有尺寸才行
}//输出
0
4
9

方便的类型定义

Eigen定义了以下矩阵类型定义：

MatrixNt 用于 Matrix<type, N, N>。例如，MatrixXi 表示 Matrix<int, Dynamic, 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（意思是Dynamic）中的任何一个。(大于4的时候需要自己改写)

t 可以是i（表示 int）、f（表示 float）、d（表示双精度）、cf（表示 complex<float>）或cd（表示 complex<double>）中的任何一种。

typedef 仅针对这五种类型定义的事实并不意味着它们是唯一受支持的标量类型。例如，支持所有标准整数类型，请参阅 Scalar types。

矩阵和向量算术（Matrix and vector arithmetic）

本章旨在提供有关如何使用Eigen在矩阵、向量和标量之间执行算术的概述和一些详细信息。

介绍

Eigen通过重载常见的 C++ 算术运算符（例如 +、-、*）或通过特殊方法（例如 dot()、cross() 等）提供矩阵/向量算术运算。对于Matrix类（矩阵和向量），运算符仅重载以支持线性代数运算。

例如，matrix1 * matrix2表示矩阵-矩阵乘积，这对于vector和scalar是不允许的。

如果要执行非线性代数的各种数组运算，得依靠数组来进行运算，后续会讲到。

加减

操作符左右对象必须具有相同的行列以及数据类型。

值得一提的是 一元运算符 - 就是 -a

void Eigen_MatrixArithmetic_AdditionAndSubtraction_001()
{Matrix2d a;a << 1, 2,3, 4;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;RowVector3d v(1, 2, 3);RowVector3d w(1, 0, 0);std::cout << "-v + w - v =\n" << -v + w - v << std::endl;
}

输出

a + b =
3 5
4 8
a - b =
-1 -12  0
Doing a += b;
Now a =
3 5
4 8
-v + w - v =
-1 -4 -6

标量乘除

so easy

void Eigen_MatrixArithmetic_ScalarMAndD_001()
{Matrix2d a;a << 1, 2,3, 4;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;
}

输出

a * 2.5 =
2.5   5
7.5  10
0.1 * v =
0.1
0.2
0.3
Doing v *= 2;
Now v =
2
4
6

关于表达式模板的说明

这虽然是一个高级主题，但是有必要在此提及。

在Eigen 中，诸如算术运算符operator+本身不执行任何计算，它们只返回一个描述要执行的计算的“表达式对象”。实际计算发生在稍后，当整个表达式被评估时，通常在operator=. 虽然这听起来很沉重，但任何现代优化编译器都能够优化掉这种抽象，结果是完美优化的代码。

简单点说就是，尽可能的将所有运算写在一起，这样才可以更好的优化代码。

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

全部写在一起的好处就是可以只遍历一次向量，下面是编译后的结果：

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

转置和共轭

转置 transpose()：aTa^TaT；共轭 conjugate()：aˉ\bar aaˉ ；伴随 adjoint ()：a∗a^*a∗；

void Eigen_MatrixArithmetic_Transposition_001()
{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;//实数矩阵 的共轭 会报错吗？ //不会MatrixXf b = MatrixXf::Random(2, 2);cout << "Here is the conjugate of b\n" << b.conjugate() << endl;
}

输出

Here is the matrix a(0.127171,-0.997497) (-0.0402539,0.170019)(0.617481,-0.613392)  (0.791925,-0.299417)
Here is the matrix a^T(0.127171,-0.997497)  (0.617481,-0.613392)
(-0.0402539,0.170019)  (0.791925,-0.299417)
Here is the conjugate of a(0.127171,0.997497) (-0.0402539,-0.170019)(0.617481,0.613392)    (0.791925,0.299417)
Here is the matrix a^*(0.127171,0.997497)    (0.617481,0.613392)
(-0.0402539,-0.170019)    (0.791925,0.299417)
Here is the conjugate of b0.64568 -0.6517840.49321  0.717887

对于实矩阵，conjugate()是无运算，因此adjoint()等价于transpose()。
aliasing issue:

例如b = a.transpose() ;，此时矩阵a在开始转置的同时就将结果写入了矩阵b，那么问题就来了，a = a.transpose() ; 就会导致，a在转置的过程中在不断的写入a，会导致程序崩坏，这种问题被称之为aliasing issue。
```
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
```
如果采用断点找到断言并注释掉进行debug，结果发现
```
Here is the matrix a:
1 2
3 4
and the result of the aliasing effect:
1 2
2 4
```
断言在Transpose.h文件中，记住用完得还原，因为有更好的方法。

转置自身

如果需要转置自身，使用 transposeInPlace() 即可；还有复数矩阵的伴随矩阵 transposeInPlace() 。

void Eigen_MatrixArithmetic_Transposition_002()
{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.transposeInPlace();cout << "and the result of the transposeInPlace():\n" << a << endl;
}

输出

Here is the matrix a:
1 2
3 4
and the result of the transposeInPlace():
1 3
2 4

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

对于操作符*不用担心 m = m * m 发生 aliasing issue，因为Eigen在计算矩阵乘法时都会引入一个临时变量temp：

m = m * m;//会被编译为
temp = m * m;
m = temp;

如果您知道您的矩阵乘积是安全的、不会造成 aliasing issue 的情况下，可以使用 noalias() 函数来避免创建临时变量；

c.noalias() += a * b;

优化问题

c.noalias() -= 2 * a.adjoint() * b;

这样的表达式已经完全进行了优化，放心使用

void Eigen_MatrixArithmetic_Multiplication_001()
{Matrix2d mat;mat << 1, 2,3, 4;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;
}

输出

Here is mat* mat:7 10
15 22
Here is mat * u:
1
1
Here is u^T * mat:
2 2
Here is u^T * v:
-2
Here is u * v^T:
-2 -02  0
Let's multiply mat by itself
Now mat is mat:7 10
15 22

点积和叉积

点积： dot() ；叉积： cross()

dot() 等价于 u.adjoint()*v

void Eigen_MatrixArithmetic_DotAndCross_001()
{Vector3d v(1, 2, 3);Vector3d w(0, 1, 2);cout << "Dot product: " << v.dot(w) << endl;double dp = v.adjoint() * w; // automatic conversion of the inner product to a scalarcout << "Dot product via a matrix product: " << dp << endl;cout << "Cross product:\n" << v.cross(w) <<    endl;
}

输出

Dot product: 8
Dot product via a matrix product: 8
Cross product:1
-21

注意

叉积仅适用于大小为3的向量！点积适用于任何大小的向量。

使用复数时，Eigen的点积在第一个变量中是共轭线性的，在第二个变量中是线性的。

基本算术归约运算

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

其中 trace() 等价于 a.diagonal().sum()

void Eigen_MatrixArithmetic_ReductionOperations_001()
{Eigen::Matrix2d mat;mat << 1, 2,3, 4;cout << "Here is mat.sum():       " << mat.sum() << endl;       // 10cout << "Here is mat.prod():      " << mat.prod() << endl;     // 24   cout << "Here is mat.mean():      " << mat.mean() << endl;      // 2.5  cout << "Here is mat.minCoeff():  " << mat.minCoeff() << endl;  // 1cout << "Here is mat.maxCoeff():  " << mat.maxCoeff() << endl;  // 4cout << "Here is mat.trace():     " << mat.trace() << endl;     // 5
}

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

void Eigen_MatrixArithmetic_ReductionOperations_002()
{Matrix3f m = Matrix3f::Random();std::ptrdiff_t i, j;//int i, j;float minOfM = m.minCoeff(&i, &j);cout << "这是矩阵 m:\n" << m << endl;cout << "它的最小系数 (" << minOfM<< ") 位于 (" << i << "," << j << ")\n\n";RowVector4i v = RowVector4i::Random();int maxOfV = v.maxCoeff(&i);cout << "这是向量 v：" << v << endl;cout << "它的最大系数 (" << maxOfV<< ") 在位置 " << i << endl;
}

输出

这是矩阵 m:0.842586    0.36906   0.853816
0.00411999  0.0541093  -0.511826-0.180944  0.0358592  -0.813776
它的最小系数 (-0.813776) 位于 (2,2)这是向量 v： -5082 -13404 -11318 -11799
它的最大系数 (-5082) 在位置 0

操作的有效性

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