潘薇鸿

keep on fighting

0%

发表于 更新于 分类于

本章介绍如何利用不同的矩阵分解方法（LU，QR，SVD，特征分解等）求解线性系统

矩阵分解

矩阵分解（decomposition/factorization）是指将矩阵拆分为若干矩阵的乘积，包括三角分解（LU）、满秩分解、QR分解、Jordan分解和奇异值分解（SVD）等。Eigen提供了常用的三种矩阵分解：LU、QR、SVD

三角分解 LU

将原方阵（square matrix）分解成一个上三角矩阵和一个下三角形矩阵。主要用于简化一个大矩阵的行列式的计算、求逆矩阵、求解方程组。这种分解法所得到的上下三角形矩阵并非唯一。

LLT分解（Cholesky分解）

若A为正定矩阵，则A有如下唯一的分解形式：

其中 $L$ 为下三角矩阵，$L^$ 为其共轭转置矩阵，当A为实矩阵时，$L^=L^T$

LDLT分解

若A为一对称矩阵且其任意k阶主子阵均不为零，即半正定或半负定，则A有如下惟一的分解形式：

其中 $L$ 为下三角形单位矩阵（即主对角线元素皆为1），$D$ 为对角矩阵，$L^T$ 为 $L$ 的转置矩阵。LDLT分解法是Cholesky分解法的改进，Cholesky分解法虽然不需要选主元，但其运算过程中涉及到开方问题，而LDLT分解法则避免了这一问题，可用于求解线性方程组。

QR分解

对列向量线性无关的矩阵 A ：

其中，$Q$ 为 $m\times m$ 的酉矩阵（实数域内称为正交矩阵），$R$ 为上三角矩阵

奇异值分解 SVD

对任意矩阵 A：

其中 $U$ 和 $V$ 均为酉矩阵（实数域内称为正交矩阵，$V^*=V^T$）且不唯一， $\Sigma$ 为对角矩阵，对角线上的元素称为A的奇异值。用于解最小二乘和数据压缩。

用矩阵分解求解线性方程组

Eigen为矩阵提供了如图所示的分解方法，调用第二列的方法后会返回一个第一列对应的类对象，再通过调用该对象的 solve() 方法实现对方程的求解，具体见后文代码示例。

矩阵分解函数

LU 分解

  • partialPivLu()

    基于部分消元的LU分解（必须是可逆方阵）。$A=PLU$ ，其中 $L$ 为单位下三角矩阵，$U$ 为上三角矩阵，$P$ 为置换矩阵，即只进行行变换

  • fullPivLu()

    基于全消元的LU分解（对任意矩阵）。$A=P^{-1}LUQ^{-1}$ ，$L,U,P$ 同上，$Q$ 为置换矩阵，即同时进行行变换和列变换。速度比部分消元慢

  • llt()

    标准Cholesky分解（必须是对称正定矩阵）。$A=LL^*$

  • ldlt()

    有主元的鲁棒Cholesky分解（必须是半正定或半负定矩阵）。$A=P^TLDL^*P$

QR 分解

  • householderQr()

    使用 Househoder 变换实现 $A=QR$ 的分解，其中 $Q$ 为酉矩阵，$R$ 为上三角矩阵。无主元，快，但不稳定

  • colPivHouseholderQr()

    基于行变换（列主元消元法）的 Householder 变换，$AP=QR$ 。较慢，但更精确

  • fullPivHouseholderQr()

    同时进行行变换和列变换的 Householder 变换，$PAP’=QR$ 。最慢，但最稳定

  • completeOrthogonalDecomposition()

    完全正交分解。可以看作是QR分解的推广

SVD 分解

  • bdcSvd()

    最快的SVD算法（推荐使用）

  • jacobiSvd()

    对小矩阵速度快，大矩阵慢

表格给出了上述所有分解方法的效率对比

代码示例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
#include <Eigen/Eigen>
#include <ctime>
#include <iostream>

using namespace std;
using namespace Eigen;

#define MATRIX_SIZE 10

int main()
{
     /* 假设矩阵A为方阵，待求解方程为 Ax=b */
     // 首先定义一个随机矩阵A，随机列向量b和待求解向量x
     MatrixXd A;
     Matrix<double, MATRIX_SIZE, 1> b;
     A = MatrixXd::Random(MATRIX_SIZE, MATRIX_SIZE);
     b = MatrixXd::Random(MATRIX_SIZE, 1);
     // 定义时间变量，方便比较效率
     clock_t start_t;
     double time;

     /* 方法一：直接右乘A的逆矩阵 x = A^{-1}*b */
     Matrix<double, MATRIX_SIZE, 1> x_inv;
     start_t = clock();
     x_inv = A.inverse()*b;
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_inv time = " << time << " ms" << endl;
     cout << "x_inv = " << x_inv.transpose() << endl;

     /* 方法二：普通LU分解 
      * Eigen 提供了四种LU分解的方法
      * - partialPivLu(): A = PLU -> 要求A必须可逆
      * - fullPivLu(): A = P^{-1}LUQ^{-1}
      * - llt(): A = LL^T -> 要求A必须对称正定
      * - ldlt(): A = LDL^T -> 要求A必须半正定或半负定
      */
     // 2.1 partialPivLu()
     Matrix<double, MATRIX_SIZE, 1> x_plu;
     start_t = clock();
     x_plu = A.partialPivLu().solve(b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_plu time = " << time << " ms" << endl;
     cout << "x_plu = " << x_plu.transpose() << endl; 
     // 2.2 fullPivLu()
     Matrix<double, MATRIX_SIZE, 1> x_flu;
     start_t = clock();
     x_flu = A.fullPivLu().solve(b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_flu time = " << time << " ms" << endl;
     cout << "x_flu = " << x_flu.transpose() << endl; 

     /* 方法三：LLT分解 */
     Matrix<double, MATRIX_SIZE, 1> x_llt;
     start_t = clock();
     // 因为llt分解要求 A 必须为正定矩阵，所以构造新的方程：(A^T*A)*x = A^T*b 来求解x
     x_llt = (A.transpose()*A).llt().solve(A.transpose()*b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_llt time = " << time << " ms" << endl;
     cout << "x_llt = " << x_llt.transpose() << endl; 

     /* 方法四：LDLT分解 */
     Matrix<double, MATRIX_SIZE, 1> x_ldlt;
     start_t = clock();
     x_ldlt = (A.transpose()*A).llt().solve(A.transpose()*b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_ldlt time = " << time << " ms" << endl;
     cout << "x_ldlt = " << x_ldlt.transpose() << endl; 

     /* 方法五：QR分解
      * Eigen 提供了四种QR分解的方法
      * - householderQr()
      * - colPivHouseholderQr()
      * - fullPivHouseholderQr()
      * - completeOrthogonalDecomposition()
      * 常用2和3
      */
     // 5.1 colPivHhouseholderQr()
     Matrix<double, MATRIX_SIZE, 1> x_cqr;
     start_t = clock();
     x_cqr = A.colPivHouseholderQr().solve(b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_cqr time = " << time << " ms" << endl;
     cout << "x_cqr = " << x_cqr.transpose() << endl; 
     // 5.2 fullPivHhouseholderQr()
     Matrix<double, MATRIX_SIZE, 1> x_fqr;
     start_t = clock();
     x_fqr = A.fullPivHouseholderQr().solve(b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_fqr time = " << time << " ms" << endl;
     cout << "x_fqr = " << x_plu.transpose() << endl; 

     /* 方法六：SVD分解
      * Eigen 提供了两种SVD分解方法
      * - bdcSvd()
      * - jacobiSvd()
      */
     // 6.1 bdcSvd()
     Matrix<double, MATRIX_SIZE, 1> x_bsvd;
     start_t = clock();
     x_bsvd = A.bdcSvd(ComputeThinU|ComputeThinV).solve(b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_bsvd time = " << time << " ms" << endl;
     cout << "x_bsvd = " << x_bsvd.transpose() << endl; 
     // 6.2 jacobiSvd()
     Matrix<double, MATRIX_SIZE, 1> x_jsvd;
     start_t = clock();
     x_jsvd = A.jacobiSvd(ComputeThinU|ComputeThinV).solve(b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_jsvd time = " << time << " ms" << endl;
     cout << "x_jsvd = " << x_jsvd.transpose() << endl; 
}

计算结果：

计算结果

计算特征值与特征向量

image-20210527002238291

如图，Eigen为矩阵提供了四种特征分解的方法，第二列表示适用矩阵的条件，第三四列分别表示运算速度和精度

以 SelfAdjointEigenSolver 为例，代码如下，其中 info() 函数用来检查特征值/特征向量是否收敛：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include <iostream>
#include <Eigen/Dense>
 
using namespace std;
using namespace Eigen;
 
int main()
{
   Matrix2f A;
   A << 1, 2, 2, 3;
   cout << "Here is the matrix A:\n" << A << endl;
   SelfAdjointEigenSolver<Matrix2f> eigensolver(A);
   if (eigensolver.info() != Success) abort();
   cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << endl;
   cout << "Here's a matrix whose columns are eigenvectors of A \n"
        << "corresponding to these eigenvalues:\n"
        << eigensolver.eigenvectors() << endl;
}

运算结果如下：

1
2
3
4
5
6
7
8
9
10
Here is the matrix A:
1 2
2 3
The eigenvalues of A are:
-0.236
  4.24
Here's a matrix whose columns are eigenvectors of A 
corresponding to these eigenvalues:
-0.851 -0.526
 0.526 -0.851

计算逆矩阵与行列式

尽管逆和行列式是基本的数学概念，但它们在数值线性代数中不像在纯数学中那么常用。求逆通常由 solve() 代替，而行列式一般并不是检测矩阵是否可逆的好方法。不过对于比较小的矩阵，逆和行列式还是比较有用的。

尽管Eigen提供了上述矩阵分解的方法，我们仍然可以直接调用 inverse() 方法和 determinant() 方法。若矩阵尺寸较小（不超过4x4），那么Eigen可以避免使用LU分解，而是使用数学公式，这样更高效。

求解最小二乘

求解最小二乘最精确的方法是使用SVD分解。Eigen提供了两种方法（见上文），推荐使用 bdcSvd() ，对大规模问题友好，同时对小规模问题能自动降级到 jacobiSvd() 方法；

另一种方法是使用 Cholesky 分解或 QR 分解，速度可能会快点，但准确率会下降；

第三种方法是将方程转化为 $A^TAx=A^Tb$ 并使用 ldlt() 方法求解。不过，如果 A 是一个病态矩阵的话，这种方法就不会太好，因为 $A^TA$ 的条件数是 $A$ 的平方，这意味着这种方法会比其他方法损失两倍数值。

将求解和构造分开

上面的例子将矩阵分解对象的构造和计算写在同一句中，当然也有方法将它们分开写，每种矩阵分解类都有一个默认的构造函数和一个 compute(matrix) 方法，该方法可以对不同的 matrix 重复调用：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include <iostream>
#include <Eigen/Dense>
 
using namespace std;
using namespace Eigen;
 
int main()
{
   Matrix2f A, b;
   LLT<Matrix2f> llt;
   A << 2, -1, -1, 3;
   b << 1, 2, 3, 1;
   cout << "Here is the matrix A:\n" << A << endl;
   cout << "Here is the right hand side b:\n" << b << endl;
   cout << "Computing LLT decomposition..." << endl;
   llt.compute(A);
   cout << "The solution is:\n" << llt.solve(b) << endl;
   A(1,1)++;
   cout << "The matrix A is now:\n" << A << endl;
   cout << "Computing LLT decomposition..." << endl;
   llt.compute(A);
   cout << "The solution is now:\n" << llt.solve(b) << endl;
}

输出：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Here is the matrix A:
 2 -1
-1  3
Here is the right hand side b:
1 2
3 1
Computing LLT decomposition...
The solution is:
1.2 1.4
1.4 0.8
The matrix A is now:
 2 -1
-1  4
Computing LLT decomposition...
The solution is now:
    1  1.29
    1 0.571

也可以通过向构造函数中传入矩阵大小来预分配分解矩阵所需要的内存空间大小：

1
2
3
HouseholderQR<MatrixXf> qr(50,50);
MatrixXf A = MatrixXf::Random(50,50);
qr.compute(A); // no dynamic memory allocation

秩显式分解

某些分解方法可以同时计算出矩阵的秩，这些通常也是在面对非满秩矩阵（方阵时称为奇异矩阵）时表现最好的方法。包括 FullPivLU 、ColPivHouseholderQR 、FullPivHouseholderQR 、BDCSVD 、JacobiSVD 、SelfAdjointEigenSolver 、ComplexEigenSolver 和 EigenSolver

这些分解类至少提供了一个 rank() 方法用来求矩阵的秩。此外，还提供了一些方便的方法，例如 isInvertible() ；有些还提供了计算矩阵的内核(零空间)和图像(列空间)方法，例如 FullPivLU ：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#include <iostream>
#include <Eigen/Dense>
 
using namespace std;
using namespace Eigen;
 
int main()
{
   Matrix3f A;
   A << 1, 2, 5,
        2, 1, 4,
        3, 0, 3;
   cout << "Here is the matrix A:\n" << A << endl;
   FullPivLU<Matrix3f> lu_decomp(A);
   cout << "The rank of A is " << lu_decomp.rank() << endl;
   cout << "Here is a matrix whose columns form a basis of the null-space of A:\n"
        << lu_decomp.kernel() << endl;
   cout << "Here is a matrix whose columns form a basis of the column-space of A:\n"
        << lu_decomp.image(A) << endl; // yes, have to pass the original A
}

输出：

1
2
3
4
5
6
7
8
9
10
11
12
13
Here is the matrix A:
1 2 5
2 1 4
3 0 3
The rank of A is 2
Here is a matrix whose columns form a basis of the null-space of A:
 0.5
   1
-0.5
Here is a matrix whose columns form a basis of the column-space of A:
5 1
4 2
3 3

当然，任何秩计算都取决于对随机阈值的选择，因为实际上没有哪个浮点矩阵是完全秩亏的。Eigen会选择一个合理的默认阈值。我们也可以自己设定想要的正确阈值：在调用 rank() 或任何其他需要使用这个阈值的方法之前，调用 setThreshold() 来设置阈值。分解计算本身，即 compute() 方法，与阈值选取无关，所以在改变阈值之后，不需要重新计算矩阵分解。

例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include <iostream>
#include <Eigen/Dense>
 
using namespace std;
using namespace Eigen;
 
int main()
{
   Matrix2d A;
   A << 2, 1,
        2, 0.9999999999;
   FullPivLU<Matrix2d> lu(A);
   cout << "By default, the rank of A is found to be " << lu.rank() << endl;
   lu.setThreshold(1e-5);
   cout << "With threshold 1e-5, the rank of A is found to be " << lu.rank() << endl;
}

输出：

1
2
By default, the rank of A is found to be 2
With threshold 1e-5, the rank of A is found to be 1

发表于 分类于

混叠问题指的是在进行赋值操作的时候，赋值表达式的两边存在重叠的矩阵区域，例如 mat = 2*mat; 或 mat = mat.transpose(); 。前一个表达式中的混叠问题不会对结果产生影响，但后一个表达式中的混叠问题会导致不可预料的结果。本节将解释什么是混叠问题，何时有害，如何处理。

举例

1
2
3
4
5
6
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 << "After the assignment:\n" << mat << endl;

输出：

1
2
3
4
5
6
7
8
mat:
1 2 3
4 5 6
7 8 9
After the assignment: 
1 2 3
4 1 2
7 4 1

可见，输出并不像我们所期望的一样将左上角的 $2\times2$ 子矩阵赋值给右下角的 $2\times2$ 子矩阵，这种情况正是由混叠问题引起的，赋值表达式的两端存在重叠的矩阵区域，即 $mat(1,1)$ 。之所以会造成这个问题，是因为Eigen中的赋值操作采用了延迟计算（也称惰性计算，Lazy Evaluation），即赋值运算并不会马上执行，而是在用到该变量时才进行计算，上述赋值运算等价于如下过程：

1
2
3
4
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) 已经不再是旧的值，而已经被赋上了新的值。

尝试进行矩阵压缩的时候很容易出现混叠问题，例如 vec = vec.head(n); 或 mat = mat.block(i,j,row,col); 等。

通常情况下，混叠问题在编译过程中是不会被识别出来的，但运行时会以异常信息 (run-time assertion) 做出提示，可以使用 EIGEN_NO_DEBUG 退出debug模式关闭异常提示，但最好不要。

如何解决混叠问题

显然，最简单的解决方法就是把右侧的计算结果赋给一个临时变量，计算结束后再将这个临时变量赋值给左侧的变量，Eigen中使用 eval() 函数实现。用 eval() 纠正上例中出现的问题：

1
2
3
4
5
6
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).eval();
cout << "After the assignment:\n" << mat << endl;

输出：

1
2
3
4
5
6
7
8
mat:
1 2 3
4 5 6
7 8 9
After the assignment: 
1 2 3
4 1 2
7 4 5

又例如对矩阵作转置时，直接使用 a = a.transpose(); 会导致混叠问题：

1
2
3
4
5
Matrix2i a; a << 1, 2, 3, 4;
cout << "Matrix a:\n" << a << endl;
 
a = a.transpose(); // !!! do NOT do this !!!
cout << "Transposed a:\n" << a << endl;

输出：

1
2
3
4
5
6
Matrix a:
1 2
3 4
Transposed a:
1 2
2 4

与上面相同的解决方法，即使用 a = a.transpose().eval() 代替。此外，Eigen提供了一种专门针对这种情况的特殊函数 transposeInPlace() 以实现将当前矩阵替换为其转置矩阵：

1
2
3
4
5
6
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;

输出：

1
2
3
4
5
6
7
Here is the initial matrix a:
1 2 3
4 5 6
and after being transposed:
1 4
2 5
3 6

类似的原地处理函数还有：

  • 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() 实现原地处理。

混叠问题与元素级操作

在元素级操作中，尽管赋值表达式两边也有可能出现同一个矩阵（例如矩阵加法，数组乘法，标量乘法等），但不会对结果产生影响，即元素级操作中的混叠是安全的，因此也就不需要 eval() ：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
MatrixXf mat(2,2); 
mat << 1, 2,  4, 7;
cout << "Here is the matrix mat:\n" << mat << endl << endl;
 
mat = 2 * mat;
cout << "After 'mat = 2 * mat', mat = \n" << mat << endl << endl;
 
 
mat = mat - MatrixXf::Identity(2,2);
cout << "After the subtraction, it becomes\n" << mat << endl << endl;
 
 
ArrayXXf arr = mat;
arr = arr.square();
cout << "After squaring, it becomes\n" << arr << endl << endl;
 
// Combining all operations in one statement:
mat << 1, 2,  4, 7;
mat = (2 * mat - MatrixXf::Identity(2,2)).array().square();
cout << "Doing everything at once yields\n" << mat << endl << endl;

输出：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
Here is the matrix mat:
1 2
4 7

After 'mat = 2 * mat', mat = 
 2  4
 8 14

After the subtraction, it becomes
 1  4
 8 13

After squaring, it becomes
  1  16
 64 169

Doing everything at once yields
  1  16
 64 169

总而言之，当表达式右侧的 (i,j) 元素仅仅依赖于表达式左侧对应的 (i,j) 元素时，赋值操作是安全的，也就不再需要多余的 eval() 操作。

混叠问题与矩阵乘法

在目标矩阵未调整大小的情况下，矩阵乘法是Eigen中唯一默认进行混叠处理的操作。换言之，若矩阵 matA 为平方型矩阵（行数 = 列数），则赋值语句 matA = matA*matA 是安全的。

除此之外的所有操作，Eigen都默认不存在由混叠导致的错误，要么是因为结果被赋值给了一个不同的矩阵，要么是因为当前进行的是元素级操作。

不过Eigen对矩阵乘法进行的这种处理是需要付出一定计算资源的代价的。在进行上述乘法的过程中，Eigen首先将乘积计算入一个临时矩阵，结束计算后再将这个临时矩阵赋值给 matA ，在当前情况下这是合理的。但相同的操作也会对 matB = matA * matA 进行，这就会导致一定程度上的计算资源浪费。在这种情况下，Eigen提供了 noalias() 函数进行直接赋值而不用经过中间的临时矩阵：matB.noalias() = matA * matA 。这种操作一定要慎用！！！

总结

当表达式两端同时出现相同的矩阵/数组时会出现混叠

  • 元素级操作时无害，包括标量/数组乘法和矩阵/数组加法
  • 两个行列相等的矩阵相乘，Eigen默认进行混叠处理；若明确不会出现混叠，使用 noalias()
  • 其他情况下Eigen默认不进行混叠处理，因此容易在某些情况下产生错误答案。使用 eval() 或 xxxInPlace() 函数来避免此类情况

发表于 更新于 分类于

Eigen 尚未提供显式方法实现变形与切片，但可以用Map类模拟实现

变形

变形指的是在保持原有矩阵系数不变的情况下，调整矩阵的行列大小，Map类提供了一种在原有存储的基础上创建不同矩阵视图的方法，需要注意的是，原矩阵中的数据存储顺序会影响数据在线性视图中的顺序：

1
2
3
4
5
6
7
8
9
10
11
MatrixXf Ma(3,3);
M1 << 1,2,3,
	  4,5,6,
	  7,8,9;

Map<RowVectorXf> v1(M1.data(), M1.size()); // M1.data()将按列优先取出数据
cout << "v1:" << endl << v1 << endl;

Matrix<float, Dynamic, Dynamic, RowMajor> M2(M1);
Map<RowVectorXf> v2(M2.data(), M2.size()); // M2.data()将按行优先取出数据
cout << "v2:" << endl << v2 << endl;

输出：

1
2
3
4
v1:
1 4 7 2 5 8 3 6 9
v2:
1 2 3 4 5 6 7 8 9
  • 例：2 $\times$ 6 矩阵转置为6 $\times$ 2：
1
2
3
4
5
Matrix<float, Dynamic, Dynamic, RowMajor> M1(2,6);
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;

输出：

1
2
3
4
5
6
7
M2:
1  7
2  8
3  9
4 10
5 11
6 12

注意到实现矩阵转置的时候，必须采用行优先保存，否则变形出来的矩阵就不是转置矩阵了！

切片

切片指的是取出矩阵中均匀间隔的行、列或元素

  • 例：间隔取出向量中的某些元素
1
2
3
4
RowVectorXf v = RowVectorXf::LinSpaced(20,0,19);
cout << "Input:" << endl << v << endl;
Map<RowVectorXf,0,InnerStride<2> > v2(v.data(), v.size()/2);
cout << "Even:" << v2 << endl;

输出：

1
2
3
Input:
 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19
Even: 0  2  4  6  8 10 12 14 16 18
  • 例：根据具体的存储顺序，使用足够的步幅 stride 依次从三列中取出一列
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
MatrixXf M1 = MatrixXf::Random(3, 8);
cout << "Column Major Input:\n"
     << M1 << endl;

Map<MatrixXf, 0, OuterStride<>> M2(
    M1.data(), M1.rows(), (M1.cols() + 2) / 3,
    OuterStride<>(M1.outerStride() * 3));
// 第二个参数 MapOptions=0 即使用默认值Unaligned
// 第三个参数 OuterStride<>说明初始化时只需声明外步幅的大小
cout << "对列优先矩阵取出每3列中的第1列:\n"
     << M2 << endl;

typedef Matrix<float, Dynamic, Dynamic, RowMajor> RowMajorMatrixXf;
RowMajorMatrixXf M3(M1);
cout << "Row Major Input:\n"
     << M3 << endl;

Map<RowMajorMatrixXf, 0, Stride<Dynamic, 3>> M4(
    M3.data(), M3.rows(), (M3.cols() + 2) / 3,
    Stride<Dynamic, 3>(M3.outerStride(), 3)); // 意味着行间步幅1，列间步幅3
// 第三个参数 Stride<OuterStride<>, InnerStride<>> 指定外步幅和内步幅
cout << "对行优先矩阵取出每3列中的第1列:\n"
     << M4 << endl;

输出：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Column Major Input:
  0.680375    0.59688  -0.329554    0.10794  -0.270431    0.83239  -0.716795  -0.514226
 -0.211234   0.823295   0.536459 -0.0452059  0.0268018   0.271423   0.213938  -0.725537
  0.566198  -0.604897  -0.444451   0.257742   0.904459   0.434594  -0.967399   0.608353
对列优先矩阵取出每3列中的第1列:
  0.680375    0.10794  -0.716795
 -0.211234 -0.0452059   0.213938
  0.566198   0.257742  -0.967399
Row Major Input:
  0.680375    0.59688  -0.329554    0.10794  -0.270431    0.83239  -0.716795  -0.514226
 -0.211234   0.823295   0.536459 -0.0452059  0.0268018   0.271423   0.213938  -0.725537
  0.566198  -0.604897  -0.444451   0.257742   0.904459   0.434594  -0.967399   0.608353
对行优先矩阵取出每3列中的第1列:
  0.680375    0.10794  -0.716795
 -0.211234 -0.0452059   0.213938
  0.566198   0.257742  -0.967399

发表于 分类于

本章介绍如何使用 Map 类通过不复制数据的方法将原始 C/C++ 数组转换成矩阵或向量

Map类模板参数与构造函数

Map类模板参数如下：

1
Map<typename MatrixType, int MapOptions = Unaligned, typename StrideType>

  • MatrixType : required ，声明是哪种类型的矩阵或向量

  • MapOptions : optional ，声明指针是对齐的还是未对齐的

  • StrideType : optional ，使用 Stride 类为内存中的数组自定义布局

    Stride 类有两个模板参数：

    1
    Stride<OuterStride, InnerStride>
    • OuterStride : 指的是 列优先矩阵的两个连续列 或 行优先矩阵的两个连续行 之间的指针增量
    • InnerStride : 指的是 列优先矩阵的某一列内两个连续行 或 行优先矩阵的某一列内两个连续列 之间的指针增量
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    int array[8];
    for(int i = 0; i < 8; ++i) array[i] = i;
    // 默认列优先
    cout << "Column-major:\n" 
         << Map<Matrix<int,2,4> >(array) << endl;
    // 通过定义矩阵类型实现行优先
    cout << "Row-major:\n" 
         << Map<Matrix<int,2,4,RowMajor> >(array) << endl;
    // 通过Stride参数实现行优先
    cout << "Row-major using stride:\n" 
         << Map<Matrix<int,2,4>, Unaligned, Stride<1,4> >(array) << endl;

    输出：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    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

    Stride 类还能更灵活：

    1
    2
    3
    4
    5
    6
    int array[24];
    for(int i = 0; i < 24; ++i) array[i] = i;
    cout << Map<MatrixXi, 0, Stride<Dynamic,2> >
            (array, 3, 3, Stride<Dynamic,2>(8, 2))
         << endl;
    // 列间隔8，行间隔2

    输出：

    1
    2
    3
    0  8 16
    2 10 18
    4 12 20

在构造一个Map类型对象的时候，还需要另外两个信息：指向数组内存空间的指针；目标矩阵/向量的大小。当矩阵/变量类型为固定大小的类型时，第二个信息可以省略。

  • 构造一个Map类型的对象用来转换动态大小的 float 矩阵：其中 pf 是一个 float* 类型的指针，指向数组内存空间，rows 和 cols 则定义了目标矩阵的维度

    1
    Map<MatrixXf> mf(pf, rows, cols);

  • 构造一个Map类型的对象用来转换固定大小的 int 只读向量：其中 pi 是一个 int* 类型的指针，因为 Vector4i 已经说明了向量维度，所以不再需要向构造函数传入向量的维度信息

    1
    Map<const Vector4i> mi(pi);

需要注意的是，Map类并没有默认构造函数，在构造对象的时候，必须告诉构造函数所要转化的数组的地址指针。

Map类型对象的使用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
typedef Matrix<float, 1, Dynamic> MatrixType;
typedef Map<MatrixType> MapType;
typedef Map<const MatrixType> MapTypeConst; // a read-only map
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()); // 通过矩阵m2mapconst只能实现对m2的访问，而不能修改

cout << "m1: " << m1 << endl;
cout << "m2: " << m2 << endl;
cout << "Squared euclidean distance: " << (m1-m2).squaredNorm() << endl;
cout << "Squared euclidean distance, using map: " 
     << (m1-m2map).squaredNorm() << endl;

m2map(3) = 7;   // 这一步会同时改变m2和m2mapconst的值，因为这三者共享内存空间
cout << "Updated m2: " << m2 << endl;
cout << "m2 coefficient 2, constant accessor: " << m2mapconst(3) << endl;

输出：

1
2
3
4
5
6
m1:  0.680375 -0.211234  0.566198   0.59688  0.823295
m2: -0.604897 -0.329554  0.536459 -0.444451   0.10794
Squared euclidean distance: 3.26291
Squared euclidean distance, using map: 3.26291
Updated m2: -0.604897 -0.329554  0.536459         7   0.10794
m2 coefficient 2, constant accessor: 7

虽然Eigen的所有函数都会像接受其他Eigen类型一样接受Map对象，但在实现自己的函数过程中，从Map对象到其对应的Dense类对象的转换并不会自动完成。具体见 以Eigen类型为参数编写函数

修改映射数组

可以使用C++的 placement new 语法实现在生命Map对象后修改其数组。简而言之，placement new 语法实现的是从给定指针所指向的地址开始创建一个新的对象，一般来说该指针指向的是一片提前申请好的内存，所以也不会进行内存分配。例如：

1
2
3
4
5
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";

输出：

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

用这种方法可以在不知道映射数组在内存中位置的情况下构造一个Map对象：

1
2
3
4
5
6
7
Map<Matrix3f> A(NULL);  // don't try to use this matrix yet!
VectorXf b(n_matrices);
for (int i = 0; i < n_matrices; i++)
{
    new (&A) Map<Matrix3f>(get_matrix_pointer(i));
    b(i) = A.trace();
}

发表于 更新于 分类于

毕设做铁路轨道检测，发现Lane Detection的大部分开源代码都是跑的公用数据集CULane（或Tusimple），而我们自己采集的数据用labelme标注完以后不是CULane的标准格式。为了能直接用公开代码训练我们自己的数据集，需要先把labelme标注的json数据转化成CULane标准格式。

项目地址：[labelme_to_culane]

CULane数据集

根据官网的介绍，CULane数据集文件包括：

1
2
3
4
5
6
7
8
- driver_23_30frame
- driver_161_90frame
- driver_182_30frame
- driver_37_30frame
- driver_100_30frame
- driver_193_90frame
- laneseg_label_w16
- list

  • 三个训练集&验证集ground truth图像与标注：

    1
    2
    3
    - driver_23_30frame
    - driver_161_90frame
    - driver_182_30frame

    每个文件夹下都包含若干视频中抽帧形成的子文件夹，每个子文件夹又包含若干图像数据文件（eg. xxxx.jpg）和一个对应的标注文件（xxxx.lines.txt），以182_30文件为例：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    - driver_182_30frame
      - 05312327_0001.MP4
        |- 00000.jpg
        |- 00000.lines.txt
        |- 00030.jpg
        |- 00030.lines.txt
        |- 00060.jpg
        |- 00060.lines.txt
        ...
      ...

    标注文件中每一行都给出了车道线上关键点的对应（x，y）坐标，以 05130.lines.txt 为例：

    1
    2
    3
    4
    104.092 590 123.101 580 144.826 570 165.645 560 187.369 550 208.189 540 229.913 530 250.733 520 272.457 510 293.277 500 315.001 490 335.82 480 357.545 470 379.269 460 399.746 450 421.471 440 442.29 430 464.015 420 485.373 410 506.192 400 527.551 390 549.275 380 570.095 370 591.452 360 613.176 350 633.996 340 655.72 330 677.001 320 698.726 310 719.545 300 740.828 290 762.552 280 783.371 270 
    661.843 590 665.482 580 669.121 570 673.124 560 677.128 550 681.131 540 684.77 530 688.773 520 692.775 510 696.778 500 700.417 490 704.419 480 708.422 470 712.06 460 716.093 450 720.096 440 724.099 430 727.918 420 731.922 410 735.562 400 739.746 390 743.385 380 747.387 370 751.206 360 755.208 350 759.211 340 763.215 330 766.854 320 770.858 310 774.862 300 778.5 290 782.503 280 786.505 270 
    1457.13 590 1438.25 580 1418.46 570 1397.78 560 1377.09 550 1356.41 540 1335.96 530 1315.27 520 1294.59 510 1273.91 500 1253.22 490 1232.77 480 1212.09 470 1192.3 460 1171.62 450 1150.94 440 1130.47 430 1109.79 420 1089.11 410 1068.42 400 1047.74 390 1027.84 380 1007.16 370 986.473 360 965.79 350 945.106 340 924.422 330 903.617 320 883.833 310 863.149 300 842.466 290 821.782 280 801.099 270 
    1676.79 450 1632.17 440 1587.27 430 1542.37 420 1497.07 410 1452.17 400 1406.89 390 1361.98 380 1317.08 370 1272.76 360 1227.86 350 1182.96 340 1137.76 330 1092.86 320 1047.96 310 1003.06 300 957.841 290 912.94 280 868.04 270

    其中，关键点选取的y坐标，是固定的，称其为行锚点 ，默认图像大小为288 * 800：

    1
    culane_row_anchor = [121, 131, 141, 150, 160, 170, 180, 189, 199, 209, 219, 228, 238, 248, 258, 267, 277, 287]

  • 测试集ground truth图像及标注：

    1
    2
    3
    - driver_37_30frame
    - driver_100_30frame
    - driver_193_90frame

    格式与训练集/验证集相同

  • 训练集/验证集车道线分割label图像：

    1
    - laneseg_label_w16

    该文件夹下为训练集/验证集中所有ground truth图像所对应的逐像素label图，即每个像素点的值表示它所在的车道线id，0表示背景。

  • 训练集/验证集/测试集列表

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    - 
      - list
        |- test_split
           |- test0_normal.txt
           |- test1_crowd.txt
           |- test2_hlight.txt
           |- test3_shadow.txt
           |- test4_noline.txt
           |- test5_arrow.txt
           |- test6_curve.txt
           |- test7_cross.txt
           |- test8_night.txt
        |- test.txt
        |- train.txt
        |- train_gt.txt
        |- val.txt
        |- val_gt.txt

    • trian.txt val.txt test.txt 都是对应训练集/验证集/测试集图像的相对路径，以 test.txt 为例：

      1
      2
      3
      4
      /driver_100_30frame/05251517_0433.MP4/00000.jpg
      /driver_100_30frame/05251517_0433.MP4/00030.jpg
      /driver_100_30frame/05251517_0433.MP4/00060.jpg
      ...

    • test_split 文件夹下的文件是将测试集 test.txt 分为不同的环境后重组，本质也是图像的相对路径

    • train_gt.txt 和 val_gt.txt 中每行的格式为：

      gt图像相对路径 逐像素label图像相对路径 四个0/1数字表示\<左2><左1><右1><右2>四个车道是否存在

      1
      2
      3
      /driver_23_30frame/05171102_0766.MP4/00020.jpg /laneseg_label_w16/driver_23_30frame/05171102_0766.MP4/00020.png 1 1 1 0
      /driver_23_30frame/05171102_0766.MP4/00050.jpg /laneseg_label_w16/driver_23_30frame/05171102_0766.MP4/00050.png 1 1 1 0
      /driver_23_30frame/05171102_0766.MP4/00080.jpg /laneseg_label_w16/driver_23_30frame/05171102_0766.MP4/00080.png 1 1 1 0

从json标注文件到label.png

labelme提供了将json文件转为label.png的方法：

1
labelme_json_to_dataset [-h] [-o output_file] json_file.json

此方法实现对一个json文件的转换，将生成以下文件：

1
2
3
4
5
6
- output_file
  |- json_file_json
     |- img.png
     |- label.png
     |- label_names.txt
     |- lable_viz.png

实现json文件批量转化：

1
2
3
4
5
6
import os
import os.path as osp

for name in os.listdir(json_files_path):
    json_data = osp.join(json_files_path, name)
    os.system("labelme_json_to_dataset json_data -o " + output_path)

如果在标注的时候不是用的纯数字，labelme在进行转换的时候会首先对标注进行字符串排序，排序后的结果在 label_names.txt 中；label.png 虽然看上去是彩色图，但用 PIL.Image 读入后会发现，这其实是个逐像素分类的图，像素值由 label_names.txt 中的标注顺序生成（下标/像素从0开始，其中0表示的是背景，即 label_names.txt 中的第一个默认分类 _background_ ）。存在的像素值对应 *_gt.txt 中的值为1，否则为0

1
2
3
4
5
6
import numpy as np
from PIL import Image

label = Image.open(label_path)
label_data = ny.asarray(label).copy()
label.save(target_path)

这样就得到了我们需要的 xxxx_label.png 文件

从label.png到 xxxx.lines.txt 和 *_gt.txt 文件

  • 根据 label.png 和 culane_row_anchor 获取轨道线上关键点的坐标生成 xxxx.lines.txt 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
# image = Image.open(image_path)
label = Image.open(label_path)
w, h = label.size

# 将行采样点按输入图像高度放缩
if h != 288:
    scale_f = lambda x : int((x * 1.0/288) * h)
    sample_tmp = list(map(scale_f,culane_row_anchor)) 
    # 根据提供的函数对指定序列做映射

lines = "" 
for i,r in enumerate(sample_tmp):
    label_r = np.asarray(label)[int(round(r))] 
    # 取出label图像中行坐标为int(round(r))的一行        
    for lane_idx in range(1, 5):
        line = ""
        pos = np.where(label_r == lane_idx)[0]
        if len(pos) == 0:            
            continue        
        pos = np.mean(pos)        
        line = line + str(pos) + " " + str(r) + " "
        # print(line)
        # cv2.circle(image, (int(round(pos)), int(round(r))), 1, (0,0,255),2)       
    lines = lines + line + "\n"
with open(lines_file_path, 'w') as f:
    f.write(lines)
  • 根据 label.png 生成 *_gt.txt 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
lines = ""
for label_path in os.listdir(label_file_path):
    label = Image.open(label_path)
    label_data = np.asarray(label).copy()
    np.unique(label_data)
    line = gt_image_path + ' ' + label_path
    for i in range(1, 5):
        if i in label_data:
            line = line + ' 1'
        else:
            line = line + ' 0'
    lines = lines + line + '\n'
with open(gt_file_path, 'w') as f:
    f.write(lines)

发表于 更新于 分类于

本章介绍Eigen中的归约( Reductions )，迭代器( Visitors ) 和广播( Broadcasting )机制，以及它们是如何应用与矩阵和数组的。

归约

归约指的是一类以矩阵或数组作为输入，返回一个标量值的函数。

常用归约函数

.sum() .prod() .mean() .minCoeff() .maxCoeff() .trace()

范数计算

  • .norm() : L-2 范数（所有系数平方和开根号）
  • .squaredNorm() : L-2 范数的平方（所有系数平方和）
  • .lpNorm<p>() : P范数（所有系数绝对值的p次幂之和的p次根），当 p=Infinity 时表示所有系数绝对值的最大值

逻辑归约函数

  • .all() : Matrix 或 Array 的所有系数均为 true 时返回 true
  • .any() : Matrix 或 Array 存在某个系数为 true 时返回 true
  • .count() : 返回 Matrix 或 Array 中为 true 的系数总个数
1
2
3
4
5
6
7
8
9
10
11
12
13
int main()
{
	ArrayXXf a(2,2); 
  	a << 1,2,
         3,4; 
  	cout << "(a > 0).all()   = " << (a > 0).all() << endl;
  	cout << "(a > 0).any()   = " << (a > 0).any() << endl;
  	cout << "(a > 0).count() = " << (a > 0).count() << endl;
  	cout << endl;
  	cout << "(a > 2).all()   = " << (a > 2).all() << endl;
  	cout << "(a > 2).any()   = " << (a > 2).any() << endl;
  	cout << "(a > 2).count() = " << (a > 2).count() << endl;
}

输出：

1
2
3
4
5
6
7
(a > 0).all()   = 1
(a > 0).any()   = 1
(a > 0).count() = 4

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

部分归约

部分归约是可以对矩阵或数组按列或行进行操作的归约，对每个列或行应用归约运算，然后返回具有相应值的列或行向量。部分归约由函数 .colwise() 和 rowwise() 实现：

1
2
3
4
5
6
7
8
9
10
int main()
{
	Eigen::MatrixXf mat(2,4);
  	mat << 1, 2, 6, 9,
           3, 1, 7, 2;
  	cout << "Column's maximum: " << endl
   		 << mat.colwise().maxCoeff() << endl;
    cout << "Row's maximum: " << endl
   		 << mat.rowwise().maxCoeff() << endl;
}

输出：

1
2
3
4
5
Column's maximum: 
3 2 7 9
Row's maximum: 
9
7

显然：列操作返回行向量，行操作返回列向量

将部分归约与其他操作结合

例：寻找系数和最大的列

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
int main()
{
  	MatrixXf mat(2,4);
  	mat << 1, 2, 6, 9,
           3, 1, 7, 2;
  
  	MatrixXf::Index   maxIndex;
  	float maxNorm = mat.colwise().sum().maxCoeff(&maxIndex);
  
  	cout << "Maximum sum at position " << maxIndex << endl;
 
  	cout << "The corresponding vector is: " << endl;
  	cout << mat.col( maxIndex ) << endl;
  	cout << "And its sum is is: " << maxNorm << endl;
}

输出：

1
2
3
4
5
Maximum sum at position 2
The corresponding vector is: 
6
7
And its sum is is: 13

迭代器

最常用的就是利用visitor机制获取矩阵/数组中最大值/最小值的位置，交给一个visitor的参数是指向存储行列位置变量的指针，变量类型为 Index ：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
int main()
{
  	Eigen::MatrixXf m(2,2);
  
  	m << 1, 2,
         3, 4;
 
  	//get location of maximum
  	MatrixXf::Index maxRow, maxCol;
  	float max = m.maxCoeff(&maxRow, &maxCol);
 
 	//get location of minimum
  	MatrixXf::Index minRow, minCol;
  	float min = m.minCoeff(&minRow, &minCol);
 
  	cout << "Max: " << max <<  ", at: " 
         << maxRow << "," << maxCol << endl;
  	cout << "Min: " << min << ", at: " 
         << minRow << "," << minCol << endl;
}

输出：

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

广播

广播机制类似于部分规约，区别在于广播构造了一个表达式，其中向量被解释成了一个矩阵，这个矩阵由向量在行/列方向上复制自身形成。

加减（+、-、+=、-=）

只能用于Vector类型的对象！！！

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
int main()
{
  	Eigen::MatrixXf mat(2,4);
  	Eigen::VectorXf v(2);
  	mat << 1, 2, 6, 9,
           3, 1, 7, 2;   
  	v << 0,
         1;
       
  	//add v to each column of m
  	mat.colwise() += v;
  
  	cout << "Broadcasting result: " << endl;
  	cout << mat << endl;
    
    Eigen::VectorXf w(4);
    w << 0,1,2,3;
       
  	//add w to each row of m
  	mat.rowwise() += v.transpose();
	  
  	cout << "Broadcasting result: " << endl;
  	cout << mat << endl;
}

输出：

1
2
3
4
5
6
Broadcasting result: 
1 2 6 9
4 2 8 3
Broadcasting result: 
1 3  8 12
4 3 10  6

在执行 mat.colwise() += v 时，相当于进行了两步操作，首先将向量 v 在行方向上复制了4次形成了一个新的4x2的矩阵，再将这个矩阵和矩阵 mat 进行加法操作：

只有向量才能进行按行列加减，不能是矩阵，否则会报编译错误；

Array类同理，只能对 ArrayXf 执行这种操作

乘除（、/、\=、/=）

按行/按列执行系数级乘法除法运算。

只能用于 Array 类型的对象！！！

如果想让矩阵第 i 列乘上向量的第 i 个系数，应该写成 mat = mat * v.asDiagonal()

将广播与其他操作结合

例：在矩阵 m 中找到与向量 v 距离最近的列向量

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#include <iostream>
#include <Eigen/Dense>
 
using namespace std;
using namespace Eigen;
 
int main()
{
  	Eigen::MatrixXf m(2,4);
  	Eigen::VectorXf v(2);
  
 	 m << 1, 23, 6, 9,
          3, 11, 7, 2;
       
  	v << 2,
         3;
 
  	MatrixXf::Index index;
  	// find nearest neighbour
  	(m.colwise() - v).colwise().squaredNorm().minCoeff(&index);
 
  	cout << "Nearest neighbour is column " << index << ":" << endl;
  	cout << m.col(index) << endl;
}

输出：

1
2
3
Nearest neighbour is column 0:
1
3

发表于 更新于 分类于

本章介绍矩阵初始化的一些高级方法，包括逗号初始化和一些特殊矩阵（单位矩阵、零矩阵）初始化

逗号初始化

  • 系数初始化顺序：从左上角开始，按从左到右，从上到下的顺序排列

  • 必须预先规定好矩阵/向量/数组的大小，否则会报错

  • 初始化列表的元素除了数值外，也可以是向量或矩阵，常用来连接几个向量或矩阵，同样必须预先规定好大小：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int main()
{
	RowVectorXd vec1(3);
	vec1 << 1, 2, 3;
	cout << "vec1 = " << vec1 << endl;
 
	RowVectorXd vec2(4);
	vec2 << 1, 4, 9, 16;
	cout << "vec2 = " << vec2 << endl;
 
	RowVectorXd joined(7);
	joined << vec1, vec2;
	cout << "joined = " << joined << endl;
}

输出：

1
2
3
vec1 = 1 2 3
vec2 =  1  4  9 16
joined =  1  2  3  1  4  9 16
  • 可以用同样的方法按块结构初始化矩阵：
1
2
3
4
5
6
7
8
int main()
{
	MatrixXf matA(2, 2);
	matA << 1, 2, 3, 4;
	MatrixXf matB(4, 4);
	matB << matA, matA/10, matA/10, matA;
	cout << matB << endl;
}

输出：

1
2
3
4
  1   2 0.1 0.2
  3   4 0.3 0.4
0.1 0.2   1   2
0.3 0.4   3   4
  • 逗号初始化也可以用来为块表达式赋值：
1
2
3
4
5
6
7
8
int main()
{
	Matrix3f m;
	m.row(0) << 1, 2, 3;
	m.block(1,0,2,2) << 4, 5, 7, 8;
	m.col(2).tail(2) << 6, 9;                   
	cout << m;
}

输出：

1
2
3
1 2 3
4 5 6
7 8 9

特殊矩阵与数组

参考手册 Quick Reference Guide - Predefined Matrices

Zero()

初始化所有系数为0，共有三种形式：

  • Zero() : 不接收任何变量，仅用于固定大小的对象；
  • Zero(size) : 接收一个变量，仅用于一维动态大小的对象；
  • Zero(rows,cols) : 接收两个变量，仅用于二维动态大小的对象。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
int main()
{
	cout << "A fixed-size array:\n";
	Array33f a1 = Array33f::Zero();
	cout << a1 << "\n\n";

	cout << "A one-dimensional dynamic-size array:\n";
	ArrayXf a2 = ArrayXf::Zero(3);
	cout << a2 << "\n\n";

	cout << "A two-dimensional dynamic-size array:\n";
	ArrayXXf a3 = ArrayXXf::Zero(3, 4);
	cout << a3 << "\n";
}

输出：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
A fixed-size array:
0 0 0
0 0 0
0 0 0

A one-dimensional dynamic-size array:
0
0
0

A two-dimensional dynamic-size array:
0 0 0 0
0 0 0 0
0 0 0 0

Ones()

初始化所有系数为1，三种形式与 Zero() 相同

Random()

初始化所有系数为随机数，三种形式与 Zero() 相同

Constant(value)

初始化所有系数为value，同样三种形式：

  • Constant(value)
  • Constant(size,value)
  • Constant(rows,cols,value)

Identity()

返回单位矩阵，仅用于 Matrix 类，包括 Identity() 和 Identity(rows,cols)

LinSpaced(size, low, high)

仅用于向量或一维数组，返回一个指定大小的向量，其系数在参数 low 和 high 之间等距分布

1
2
3
4
5
6
7
8
9
10
int main()
{
	ArrayXXf table(10, 4);
	table.col(0) = ArrayXf::LinSpaced(10, 0, 90);
	table.col(1) = M_PI / 180 * table.col(0);
	table.col(2) = table.col(1).sin();
	table.col(3) = table.col(1).cos();
	cout << "  Degrees   Radians      Sine    Cosine\n";
	cout << table << std::endl;
}

输出：

1
2
3
4
5
6
7
8
9
10
11
Degrees   Radians      Sine    Cosine
      0         0         0         1
     10     0.175     0.174     0.985
     20     0.349     0.342      0.94
     30     0.524       0.5     0.866
     40     0.698     0.643     0.766
     50     0.873     0.766     0.643
     60      1.05     0.866       0.5
     70      1.22      0.94     0.342
     80       1.4     0.985     0.174
     90      1.57         1 -4.37e-08

从这些例子可以看出，上述方法都会返回一个矩阵或数组再赋值给变量（或表达式），方便起见，Eigen还定义了一些函数来方便地执行这些操作，例如 setZero() MatrixBase::setIdentity() DenseBase::setLinSpaced() 等。下面这个例程展示了如何使用静态方法、逗号初始化和 setXxx() 函数来初始化矩阵：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
const int size = 6;
MatrixXd mat1(size, size);
mat1.topLeftCorner(size/2, size/2)     = MatrixXd::Zero(size/2, size/2);
mat1.topRightCorner(size/2, size/2)    = MatrixXd::Identity(size/2, size/2);
mat1.bottomLeftCorner(size/2, size/2)  = MatrixXd::Identity(size/2, size/2);
mat1.bottomRightCorner(size/2, size/2) = MatrixXd::Zero(size/2, size/2);
std::cout << mat1 << std::endl << std::endl;
 
MatrixXd mat2(size, size);
mat2.topLeftCorner(size/2, size/2).setZero();
mat2.topRightCorner(size/2, size/2).setIdentity();
mat2.bottomLeftCorner(size/2, size/2).setIdentity();
mat2.bottomRightCorner(size/2, size/2).setZero();
std::cout << mat2 << std::endl << std::endl;
 
MatrixXd mat3(size, size);
mat3 << MatrixXd::Zero(size/2, size/2), MatrixXd::Identity(size/2, size/2),
        MatrixXd::Identity(size/2, size/2), MatrixXd::Zero(size/2, size/2);
std::cout << mat3 << std::endl;

输出：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0

0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0

0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0

除此之外，Eigen还定义了一些用于生成基向量的方法：

  • Vector3f::UnitX() : (1, 0, 0)
  • Vector3f::UnitY() : (0, 1, 0)
  • Vector3f::UnitZ() : (0, 0, 1)
  • VectorXf::Unit(size,i) :
    • 例如 ：VectorXf::Unit(4,1) == Vector4f(0,1,0,0) == Vector4f::UnitY()

用作临时对象

无论是逗号初始化程序还是上面介绍的几种类内静态方法，都可以用作表达式中的临时对象。

逗号初始化程序用作临时对象时，需要调用 finished() 方法获取实际的矩阵对象。

1
2
3
4
5
6
7
8
9
10
11
12
13
int main()
{
	MatrixXd m = MatrixXd::Random(3,3);
    m = (m + MatrixXd::Constant(3,3,1.2)) * 50;
    cout << "m =" << endl << m << endl;
    VectorXd v(3);
    v << 1, 2, 3;
    cout << "m * v =" << endl << m * v << endl;
    MatrixXf mat = MatrixXf::Random(2, 3);
	cout << "mat =\n" << mat << endl;
	mat = (MatrixXf(2,2) << 0, 1, 1, 0).finished() * mat;
	cout << "mat =\n" << mat << endl;
}

输出：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
m =
  10 55.9 14.7
23.2 63.3 77.9
85.6 31.9 77.9
m * v =
166
383
383
mat =
     -1   0.511  0.0655
 -0.737 -0.0827  -0.562
mat =
 -0.737 -0.0827  -0.562
     -1   0.511  0.0655

发表于 更新于 分类于

本章介绍 Matrix 和 Array 中的块操作。

块操作 .block()

取一个从 (i,j) 开始的，大小为 (p,q) 的块，有两种语法：

  • matrix.block(i,j,p,q); => 构造一个动态大小的块
  • matrix.block<p,q>(i,j); => 构造一个固定大小的块

两种语法的执行结果相同，都能用在固定大小/动态大小的矩阵/数组上。唯一不同的是当块大小已知且比较小的时候，返回固定大小的块能让代码执行得更快一点。

尽管 .block() 方法可以用于任何块操作，但对于特殊情况，使用更具针对性的方法能提供更好的性能，下面将介绍三种特殊方法。

行列操作

  • 取第 i 行：matrix.row(i)
  • 取第 j 列：matrix.col(j)

与角相关的操作

同样也有返回动态大小/固定大小两种语法

  • 左上角 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>();
  • 取后 p 列：matrix.rightCols(q); matrix.rightCols<q>(); 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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;
}

输出：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
m.leftCols(2) =
 1  2
 5  6
 9 10
13 14

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

After assignment, m = 
 8 12 16  4
 5  6  7  8
 9 10 11 12
13 14 15 16

对向量和一维数组的块操作

同样也有返回动态大小/固定大小两种语法

  • 包含前 n 个元素：vector.head(n); vector.head<n>();
  • 包含后 n 个元素：vector.tail(n); vector.tail<n>();
  • 从 i 开始，包含 n 个元素：vector.segment(i,n); vector.segment<n>(i); 
1
2
3
4
5
6
7
8
9
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;
}

输出：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
	
v.head(3) =
1
2
3

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

after 'v.segment(1,4) *= 2', v =
 1
 4
 6
 8
10
 6

发表于 更新于 分类于

本章介绍Eigen中的Array类。

Matrix类主要面向线性代数，Array类则提供更普遍意义上的数组操作，允许进行没有具体线代意义的系数级运算，例如给每个系数加上一个常数，或将两个数组按系数相乘等。

类型

Array的六个模板参数和Matrix一模一样，并且也定义了一些常用类型

  • ArrayNt ：表示一个一维数组，N 为数组大小，t 为数值类型
  • ArrayNNt ：表示一个二维数组
1
2
3
4
typedef Array<float, Dynamic, 1> ArrayXf;
typedef Array<float, 3, 1> Array3f;
typedef Array<double, Dynamic, Dynamic> ArrayXXd;
typedef Array<double, 3, 3> Array33d;

访问系数

与Matrix相同，Array也是通过重载括号运算符 () 实现对数组系数的读写操作；重载 << 操作符逗号初始化数组或打印数组。

加减法

与Matrix相同，当两个数组大小相同时操作有效，且加减法按系数进行；

与Matrix不同，Array也支持 array+scalar 的操作，即向Array中的每个系数加上同一个常数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include <iostream>
#include <eigen3/Eigen/Eigen>
 
using namespace Eigen;
using namespace std;
 
int main()
{
  ArrayXXf a(3,3);
  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
  cout << "a + b = " << endl << a + b << endl << endl;
 
  // Subtracting a scalar from an array
  cout << "a - 2 = " << endl << a - 2 << endl;
}

输出：

1
2
3
4
5
6
7
8
9
a + b = 
 2  4  6
 5  7  9
 8 10 12

a - 2 = 
-1  0  1
 2  3  4
 5  6  7

乘法

  • array * scalar 或 scalar * array ：同Matrix
  • array * array ：与Matrix不同。Matrix进行的是矩阵乘法，Array则是系数对应相乘。因此当且仅当两个数组维度相同时，才能相乘。
1
2
3
4
5
6
7
8
9
10
int main()
{
  ArrayXXf a(2,2);
  ArrayXXf b(2,2);
  a << 1,2,
       3,4;
  b << 5,6,
       7,8;
  cout << "a * b = " << endl << a * b << endl;
}

输出：

1
2
3
a * b = 
 5 12
21 32

其他系数级操作

参考Eigen手册 Quick Reference Guide : Coefficient-wise & Array operators

Array 与 Matrix 类型转换

  • Array -> Matrix : A.matrix()
  • Matrix -> Array : M.array()

Eigen禁止在一个表达式中混用Array和Matrix，唯一例外的是 赋值运算 ，它允许Array和Matrix的互相赋值。

Eigen也为Matrix提供了一些 cwise*() 的方法使其允许直接进行系数级操作，类似 .array().*()

发表于 更新于 分类于

这章介绍Eigen中矩阵、向量和标量之间的运算。对Matrix类，运算符重载仅支持线性代数运算，例如矩阵x矩阵，而像向量+标量的操作是不允许的。Array类则允许执行所有的数组操作，见 The Array Class and Coefficient-wise Operations 。

加减法

Eigen不会自动执行类型转换，所以左右维度和数值类型必须保证相同。可用的操作符包括：

  • 一元操作符 - : -a
  • 二元操作符 + - : a+b, a-b
  • 复合操作符 += -= : a+=b, a-=b

标量乘除

矩阵所有系数乘/除以该标量。可用操作符包括：

  • 二元操作符 * : matrix scalar, scalar matrix
  • 二元操作符 / : matrix / scalar
  • 复合操作符 *= /= : matrix *= scalar, matrix /= scalar

转置与共轭

  • transpose() : $\boldsymbol{a}^T$ （转置）
  • conjugate() : $\overline{\boldsymbol{a}}$ （共轭）
  • adjoint() : $\boldsymbol{a}^*$ （共轭转置）

显然，对实矩阵而言，conjugate() 并不执行什么操作，所以 adjoint() 完全等同于 transpose() 。

需要非常注意的一点是，基本运算符 transpose() 和 adjoint() 仅仅返回一个代理对象，而并未进行实质的运算。例如，执行 b = a.transpose() 时，Eigen在将结果写入b的同时执行转置；从而引出另一个问题，如果执行 a = a.transpose() ，Eigen同样会在转置完成前开始将结果写入矩阵a，因此这句话不会像预期一样将a变为它的转置，反而会在执行过程中得到 aliasing issue 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
#include <iostream>
#include <eigen3/Eigen/Eigen>

using namespace std;

int main()
{
    Eigen::Matrix3i a;
    Eigen::Matrix3i b;
    a << 1, 2, 3,
        4, 5, 6,
        7, 8, 9;
    cout << "origin matrix a:\n"
         << a << endl;

    b = a.transpose();
    cout << "b as the transpose of matrix a:\n"
         << b << endl;

    a = a.transpose(); // !!! DO NOT DO THIS !!!
    cout << "Here is the result of the aliasing effect:\n"
         << a << endl;

    return 0;
}

输出：

1
2
3
4
5
6
7
8
9
10
origin matrix a:
1 2 3
4 5 6
7 8 9
b as the transpose of matrix a:
1 4 7
2 5 8
3 6 9
main: /usr/include/eigen3/Eigen/src/Core/Transpose.h:378: static void Eigen::internal::checkTransposeAliasing_impl<Derived, OtherDerived, MightHaveTransposeAliasing>::run(const Derived&, const OtherDerived&) [with Derived = Eigen::Matrix<int, 3, 3>; OtherDerived = Eigen::Transpose<Eigen::Matrix<int, 3, 3> >; bool MightHaveTransposeAliasing = true]: Assertion `(!check_transpose_aliasing_run_time_selector <typename Derived::Scalar,blas_traits<Derived>::IsTransposed,OtherDerived> ::run(extract_data(dst), other)) && "aliasing detected during transposition, use transposeInPlace() " "or evaluate the rhs into a temporary using .eval()"' failed.
Aborted (core dumped)

所谓混叠问题，即 aliasing issue ，指的是在赋值表达式的左右两边存在矩阵的重叠区域，这种情况下，有可能得到非预期的结果，具体见 Aliasing 。

如果想执行就地转置（in-place transposition），调用 transposeInPlace() 方法：

1
2
3
a.transposeInPlace();
cout << "After being transposed:\n"
     << a << endl;

输出：

1
2
3
4
After being transposed:
1 4 7
2 5 8
3 6 9

同理，也有为复矩阵准备的就地共轭转置 adjointInPlace() 方法

矩阵-矩阵乘法乘法

向量是一种特殊的矩阵，因此矩阵-向量、向量-向量乘法也包括其中。可用操作符有：

  • 二元操作符 * : a*b
  • 复合操作符 *= : a*=b

执行操作 m = m * m 是安全的，因为Eigen将矩阵乘法当作一种特殊形式进行处理，会在编译时引入一个临时变量，所以不必担心会出现 aliasing issue：

1
2
tmp = m * m;
m = tmp;

当然，如果你很清楚自己的矩阵乘法不会有引起 aliasing issue 的可能，也可以调用 noalias() 方法来避免引入临时变量：

1
c.noalias() += a*b;

向量点乘和叉乘

dot() 和 corss() 方法分别实现向量点乘和叉乘，当然，点乘也可以通过 u.adjoint()*v 实现，只不过返回的结果是个1x1的矩阵：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#include <iostream>
#include <eigen3/Eigen/Eigen>

using namespace std;
using namespace Eigen;
int main()
{
    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 scalar
    cout << "Dot product via a matrix product: " << dp << endl;
    cout << "Cross product:\n"
         << v.cross(w) << endl;
    return 0;
}

输出：

1
2
3
4
5
6
Dot product: 8
Dot product via a matrix product: 8
Cross product:
 1
-2
 1

需要非常注意的是，叉乘只能用于大小为3的向量 ；点乘则可用于任意大小的向量

基本算数归约运算

Eigen提供了一些归约运算（reduction operation），将给定的矩阵或向量归约为单个值，例如 sum() , prod() , maxCoeff() , minCoeff() , trace() 分别表示所有矩阵系数的和、积、最大值、最小值，以及矩阵的迹，矩阵的迹也可以通过 m.diagonal().sum() 求得。

通过向 maxCoeff() 或 minCoeff() 中传入 std::ptrdiff_t 类型的参数指针，可以得到对应最大值/最小值的坐标：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <iostream>
#include <eigen3/Eigen/Eigen>

using namespace std;
using namespace Eigen;
int main()
{
    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;
    return 0;
}

输出：

1
2
3
4
5
6
7
8
Here is the matrix m:
 0.680375   0.59688 -0.329554
-0.211234  0.823295  0.536459
 0.566198 -0.604897 -0.444451
Its minimum coefficient (-0.604897) is at position (2,1)

Here is the vector v:  115899597  -48539462  276748203 -290373134
Its maximum coefficient (276748203) is at position 2

操作有效性

Eigen会在执行操作的时候检查有效性。

重要的编译错误信息会用大写突出出来。

但更多时候（例如检查动态大小的矩阵的操作），编译过程不会检查出错误。Eigen将在 debug mode 下使用 runtime assertions 让程序在运行过程中碰到非法操作的时候退出，同时报出错误信息。如果关掉assertions，程序很有可能会崩溃。