@taqikema 2020-02-05T16:32:51.000000Z 字数 8299 阅读 3501

Eigen学习笔记

Eigen 学习笔记

Eigen是使用 C++开发的矩阵计算开源库，因为其中使用了模板类实现，不能分离式编译，所以直接对用户提供了源码。因此，Eigen的使用也非常简单，直接下载源码，在工程中包括相应的头文件即可。
Eigen的官网上有详细的说明文档，具体的类型、用法和示例在上面都能找到，强烈建议想要学习 Eigen使用和查找某种语法时，直接去翻阅官方文档即可。不过，因为官方文档是英文的，有博主翻译了大部分的内容，并将其放在了 GitHub上。优点是中文写的，看起来很快，并且通过示例+注释方式讲解具体用法；缺点是有些原文中关键或需要格外注意的地方被作者删掉了。有需要的同学自取（中文翻译）。
想要快速上手的话，也可以参照下面两篇博客。
Eigen: C++开源矩阵计算工具——Eigen的简单用法
 C++矩阵处理工具——Eigen

具体的用法官方文档都有详细解释，所以后面部分主要是我在阅读官方文档过程中，记录了一些自己觉得重要或容易搞混的事项。

0.下标从 0开始！
1.resize操作，如果新矩阵的大小与旧矩阵一致，则不执行任何操作。否则会重新设置大小，并清空原有元素。想要不清空原有元素并 resize，使用 conservativeResize方法。考虑到 API的一致性，固定矩阵也支持了 resize操作。新大小与旧大小一致，无操作；否则，产生一个断言失败。
2.固定尺寸与动态尺寸。什么时候应该使用固定尺寸（例如Matrix4f），什么时候应该使用动态尺寸（例如MatrixXf）？简单的答案是：在可能的地方使用固定尺寸（分配在栈上，大小太大，可能导致内存溢出）来显示非常小的尺寸，在需要的地方使用动态尺寸（分配在堆上）来显示较大的尺寸。
3.Eigen定义了以下Matrix typedef：

 MatrixNt for Matrix<type, N, N>. For example, MatrixXi for Matrix<int, Dynamic, Dynamic>.
 VectorNt for Matrix<type, N, 1>. For example, Vector2f for Matrix<float, 2, 1>.
 RowVectorNt for Matrix<type, 1, N>. For example, RowVector3d for Matrix<double, 1, 3>.
 // N可以是任何一个2，3，4，或X（意思Dynamic）。
 // t可以是i（表示int），f（表示float），d（表示double），cf（表示complex <float>）或cd（表示complex <double>）的任何一种。

4.因为“别名问题”，不能使用 a = a.transpose()，需要使用的话，请使用 m.transposeInPlace()。
5.m=m*m是可行的，不会引发“别名问题”，这是因为 Eigen在实现矩阵相乘操作时使用了一个临时变量。如下所示，

tmp = m*m;
m = tmp;

6.向量叉积只适用于大小为3的向量，点积适用于任意向量。
7.minCoeff/maxCoeff方法包含参数，不仅能返回最值本身，还能返回下标。注意，传递给访问者的参数是指向要存储行和列位置的变量的指针，这些变量应为Index类型。
8. 数组支持array + scalar的形式，将标量添加到数组中的每个系数。这提供了不能直接用于Matrix对象的功能。
9. 矩阵将乘法解释为矩阵乘积，而数组将乘法解释为按系数乘积。因此，当且仅当两个数组具有相同的维数时，它们才能相乘.
10. 数组的运算单位为每一个元素！，a.min(b)，a/b数组对应元素比较，得到由较小值构成的数组。
11. 矩阵适用线性代数运算，数组适用元素运算。矩阵和数组可以通过 .array和 .matrix相互转化。Eigen禁止在表达式中混合矩阵和数组。例如，不能直接矩阵和数组相加。运算符+的操作数要么都是矩阵，要么都是数组。即使是矩阵，也可以使用.cwiseProduct()方法来计算元素乘积。
12. 块表达式既可以用作右值，也可以用作左值。
13. 逗号初始化器，必须先设置大小！另外，其初始化列表的元素本身还可以是向量或矩阵。
14. finished方法的作用是立即生成该对象。
15. 范数操作。

Code:
    VectorXf v(2);
    MatrixXf m(2, 2), n(2, 2);
    v << -1,
          2;
    m << 1, -2,
        -3, 4;
    // 向量范数
    cout << "v.squaredNorm() = " << v.squaredNorm() << endl;
    cout << "v.norm() = " << v.norm() << endl;
    cout << "v.lpNorm<1>() = " << v.lpNorm<1>() << endl;
    cout << "v.lpNorm<Infinity>() = " << v.lpNorm<Infinity>() << endl;
    cout << endl;
    // 矩阵范数
    cout << "m.squaredNorm() = " << m.squaredNorm() << endl;
    cout << "m.norm() = " << m.norm() << endl;
    cout << "m.lpNorm<1>() = " << m.lpNorm<1>() << endl;
    cout << "m.lpNorm<Infinity>() = " << m.lpNorm<Infinity>() << endl;
Output:
    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

16.column-wise返回行向量，row-wise返回列向量。
17.广播，相当于使向量（列或行）通过在一个方向上复制而被解释为矩阵。在使用Matrix进行操作时，广播操作只能应用于Vector类型的对象。
18.Map 类实现C++中的数组内存和Eigen对象的交互，指向定义系数数组的内存区域的指针，以及所需的矩阵或矢量形状。如 Map<Matrix<int, 2, 4>>(array)、Map<Matrix<float, 1, Dynamic>> m2map(p, m.size())和 Map<const Matrix<float, 1, Dynamic>> m2map(p, m.size())。
19.通过C++ 位置new语法，可以改变 Map对象的元素，具体用法如下所示：

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

20.Eigen没有提供方便的方法来切片或重塑矩阵。但是，可以使用Map类轻松模拟这些功能。下面的例子演示了如何将2x6矩阵重塑为6x2矩阵。

Code:
MatrixXf M1(2,6);    // Column-major storage
M1 << 1, 2, 3,  4,  5,  6,
      7, 8, 9, 10, 11, 12;
Map<MatrixXf> M2(M1.data(), 6,2);
cout << "M2:" << endl << M2 << endl;
Output:
M2:
 1  4
 7 10
 2  5
 8 11
 3  6
 9 12

21.根据存储顺序，可以使用 outerStride/innerStride来选取矩阵（向量）中每隔 P列中的几列（个）元素值。目前还没有全看懂 Stride的用法

Code:
MatrixXf M1 = MatrixXf::Random(3,8);
cout << "Column major input:" << endl << M1 << "\n";
Map<MatrixXf,0,OuterStride<> > M2(M1.data(), M1.rows(), (M1.cols()+2)/3, OuterStride<>(M1.outerStride()*3));
cout << "1 column over 3:" << endl << M2 << "\n";
typedef Matrix<float,Dynamic,Dynamic,RowMajor> RowMajorMatrixXf;
RowMajorMatrixXf M3(M1);
cout << "Row major input:" << endl << M3 << "\n";
Map<RowMajorMatrixXf,0,Stride<Dynamic,3> > M4(M3.data(), M3.rows(), (M3.cols()+2)/3, Stride<Dynamic,3>(M3.outerStride(),3));
cout << "1 column over 3:" << endl << M4 << "\n";
Output:
Column major input:
   0.68   0.597   -0.33   0.108   -0.27   0.832  -0.717  -0.514
 -0.211   0.823   0.536 -0.0452  0.0268   0.271   0.214  -0.726
  0.566  -0.605  -0.444   0.258   0.904   0.435  -0.967   0.608
1 column over 3:
   0.68   0.108  -0.717
 -0.211 -0.0452   0.214
  0.566   0.258  -0.967
Row major input:
   0.68   0.597   -0.33   0.108   -0.27   0.832  -0.717  -0.514
 -0.211   0.823   0.536 -0.0452  0.0268   0.271   0.214  -0.726
  0.566  -0.605  -0.444   0.258   0.904   0.435  -0.967   0.608
1 column over 3:
   0.68   0.108  -0.717
 -0.211 -0.0452   0.214
  0.566   0.258  -0.967

22.在Eigen中，混叠(aliasing)是指相同的矩阵（或数组或向量）出现在赋值运算符的左侧和右侧的赋值语句;例如，A = AB , a = a^Tb A=A*A。产生混叠的原因是 Eigen采用惰性求值。混叠可能是有害的，也可能是无害的,有害的混叠，可能导致不正确的结果，无害的混叠可以产生正确的结果。

无害混叠是我们不需要评估，它对系数计算无害。这包括标量乘法和矩阵或数组加法。
将两个矩阵相乘时，Eigen会假定发生混叠(注意，Eigen3.3以后的版本中需要区分目标矩阵大小是否改变了)。如果您知道没有混叠，则可以使用noalias()，这样执行效率会更高。
在所有其他情况下，Eigen假定不存在混叠问题，因此如果实际上发生混叠，则会给出错误的结果。为防止这种情况，您必须使用eval（）或xxxInPlace（）函数之一。
有害混叠，使用.eval()方法，可以解决混叠问题。具体的说，.eval()方法生成临时对象，然后再执行赋值操作，而如果Eigen中假定该操作是混淆的，那么它会自动使用.eval()方法，而不需要我们显示调用。
错误示例：

a = a.transpose();

正确示例：

a = a.transpose().eval();
或 a.transposeInPlace()

可以的话，最好使用 xxxInPlace()函数，Eigen提供了以下xxxInPlace()函数：

Original function                       In-place function
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()。
23.当矩阵以行优先顺序存储时，逐行遍历矩阵的算法会更快，因为数据位置更好。同样，对于主要列矩阵，逐列遍历更快。而因为 Eigen默认矩阵以列顺序存储，所以大部分算法都是对于按列存储的矩阵，效率更高。
24.当结构体中具有固定大小的Eigen对象成员时，不能直接 new该结构体类型的对象。原因是 Eigen对固定大小且大小是 16字节倍数的对象，会进行16字节对齐。解决办法是将 EIGEN\_MAKE\_ALIGNED\_OPERATOR_NEW宏放在类的 public部分。具体示例如下所示：
错误示例：

class Foo
{
  ...
  Eigen::Vector2d v;
  ...
};
...
Foo *foo = new Foo;

正确示例：

class Foo
{
  double x;
  Eigen::Vector2d v;
public:
  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
};

25.将Stl容器（例如std :: vector，std :: map等）与Eigen对象或包含Eigen对象的类一起使用时，也需要额外注意。解决办法为：1.使用aligned_allocator；2.vector容器还需要额外添加#include 。
错误示例：

std::map<int, Eigen::Vector4f>

正确示例：

std::map<int, Eigen::Vector4f, std::less<int>, 
         Eigen::aligned_allocator<std::pair<const int, Eigen::Vector4f> > >

26.通过值传递Eigen对象时，程序会出错。这是必须使用引用来传递 Eigen对象！
27.编译器会对堆栈对齐做出错误假设。Windows上使用GCC（例如MinGW或TDM-GCC）时可能会出现该问题，即编译器不能总是爆炸16字节对齐。这应该是 GCC的 bug，并且在 GCC4.5中已经修复。
28.不同分解方式的特点：
多种分解方式
29.计算特征值之前，先使用 info函数确定计算是否成功！
30.最小二乘求解的最准确方法是SVD分解。Eigen提供了两种实现。推荐的对象是BDCSVD类，它可以很好地解决较大的问题，并自动退回到JacobiSVD类以解决较小的问题。另外，在进行矩阵分解时，可以先构造对象，在计算分解后的矩阵，这时可以调用 .compute方法。
31.任何秩计算都取决于对任意阈值的选择，可以在调用 rank()方法之前在分解对象上调用setThreshold（）来进行设置。
32.计算线性方程组的最小二乘解，可供选择的三种方法是SVD分解，QR分解和正态方程。其中，SVD分解通常最准确但最慢，正则方程(normal equations)最快但最不准确，QR分解（有3种QR分解类：HouseholderQR（无pivoting，因此快速但不稳定），ColPivHouseholderQR（列枢轴，因此较慢但更准确）和FullPivHouseholderQR（全枢轴，因此最慢且最稳定）。）介于两者之间。
33.complete pivoting（即full pivoting，全主元），就是在矩阵分解或高斯消元的过程中，主元是在未进行变换的所有行和列之间进行选择。也就是需要同时进行行交换和列交换。partial pivoting（列主元）就是只在当前进行变换的列中选择主元，只需要进行行交换。
34.从Eigen 3.3开始，LU，Cholesky和QR分解可以就地进行操作，即直接在给定的输入矩阵内进行。当处理大量矩阵或可用内存非常有限（嵌入式系统）时，此功能特别有用。为此，必须使用Ref <>矩阵类型实例化各个分解类，并且必须使用输入矩阵作为参数来构造分解对象。示例如下：

MatrixXd A(2, 2);
A << 2, -1, 1, 3;
cout << "Here is the input matrix A before decomposition:\n"
             << A << endl;
// Output is:
// Here is the input matrix A before decomposition:
//  2 -1
//  1  3
 PartialPivLU<Ref<MatrixXd>> lu(A);
cout << "Here is the input matrix A after decomposition:\n"
             << A << endl;
// Output is:
// Here is the input matrix A after decomposition:
//   2  -1
// 0.5 3.5
//然后，该lu对象可以像往常一样使用，例如解决Ax = b问题。但是在验证残差时，必须声明一个新矩阵 A0：
MatrixXd A0(2, 2);
A0 << 2, -1, 1, 3;
VectorXd b(2);
b << 1, 2;
VectorXd x = lu.solve(b);
cout << "Residual: " << (A0 * x - b).norm() << endl;
// Output is:
// Residual: 0
//由于在A和lu之间共享内存，因此修改矩阵A将导致lu无效。可以通过修改的内容A并尝试再次解决初始问题来轻松验证这一点：
A << 3, 4, -2, 1;
x = lu.solve(b);
// Output is:
// Residual: 15.8114
//请注意，调用compute不会更改该lu对象引用的内存。因此，如果使用A1不同于的另一个矩阵调用计算方法A，则A1不会修改内容。这仍然A是将用于存储矩阵的L和U因子的内容A1。可以很容易地验证如下：
MatrixXd A1(2, 2);
A1 << 5, -2, 3, 4;
lu.compute(A1);
cout << "Here is the input matrix A1 after decomposition:\n"
             << A1 << endl;
// Output is:
// Here is the input matrix A1 after decomposition:
//  5 -2
//  3  4
//矩阵A1是不变的，因此可以求解A1 *x = b，直接检查残差而无需任何副本A1 :
x = lu.solve(b);
cout << "Residual: " << (A1 * x - b).norm() << endl;
// Output is:
/ /Residual: 2.48253e-16

35.对于不同分解方案的执行速度，有如下结论：
1.LLT始终是最快的解决方案。
2.对于很大程度上过度约束的问题，Cholesky/LU分解的成本主要由对称协方差矩阵的计算决定。
3.对于较大的问题，只有实现高速缓存友好的分块策略的分解才能很好地进行扩展。这些包括LLT，PartialPivLU，HouseholderQR和BDCSVD。这解释了为什么对于4k x 4k矩阵，HouseholderQR比LDLT更快。将来，LDLT和ColPivHouseholderQR也将实施分块策略。
4.CompleteOrthogonalDecomposition基于ColPivHouseholderQR，因此可以达到相同的性能水平。
36.Eigen的Geometry模块提供了两种不同的几何变换：
1. 抽象的变换，例如旋转（轴角或四元数表示），平移，缩放。这些转换未表示为矩阵，但是您仍然可以将它们与表达式中的矩阵和向量混合，并根据需要将它们转换为矩阵。
2. 射影或仿射变换矩阵：请参见Transform类。这些确实是矩阵。注意, 如果要使用OpenGL 4x4矩阵，则需要Affine3f和Affine3d。由于Eigen默认为列序存储，因此您可以直接使用Transform :: data()方法将转换矩阵传递给OpenGL。可以从抽象转换构造一个Transform，但是不能在定义的同时直接赋值，这是因为 C++会默认调用隐式转换构造函数，而 Eigen并没有定义这么一个函数。即，以下的写法是错误的。

Transform t = AngleAxis(angle,axis);    /* 错误 */
Transform t(AngleAxis(angle,axis));     /* 正确 */

Eigen学习笔记

内容目录