少帥的博客

Eigen 快速参考指南

发表于 更新于 分类于
本文字数： 14k
本文为 Eigen 官方文档 Quick reference guide 的搬运与简单翻译.

原文链接: https://eigen.tuxfamily.org/dox/group__QuickRefPage.html

模块与头文件

Eigen 库分为一个核心 (Core) 模块和几个附加模块. 每个模块都有一个相应的头文件, 必须包含头文件才能使用对应的模块. 提供 Dense 和 Eigen 头文件是为了方便地一次访问多个模块.

Module Header file Contents
Core #include Matrix and Array classes, basic linear algebra (including triangular and selfadjoint products), array manipulation
Geometry #include Transform, Translation, Scaling, Rotation2D and 3D rotations (Quaternion, AngleAxis)
LU #include Inverse, determinant, LU decompositions with solver (FullPivLU, PartialPivLU)
Cholesky #include LLT and LDLT Cholesky factorization with solver
Householder #include Householder transformations; this module is used by several linear algebra modules
SVD #include SVD decompositions with least-squares solver (JacobiSVD, BDCSVD)
QR #include QR decomposition with solver (HouseholderQR, ColPivHouseholderQR, FullPivHouseholderQR)
Eigenvalues #include Eigenvalue, eigenvector decompositions (EigenSolver, SelfAdjointEigenSolver, ComplexEigenSolver)
Sparse #include Sparse matrix storage and related basic linear algebra (SparseMatrix, SparseVector) (see Quick reference guide for sparse matrices for details on sparse modules)
#include Includes Core, Geometry, LU, Cholesky, SVD, QR, and Eigenvalues header files
#include Includes Dense and Sparse header files (the whole Eigen library)

数组, 矩阵和向量类型

Eigen 提供了两类 Dense 对象: 由模板类 Matrix 表示的数学矩阵和向量; 以及由模板类 Array 表示的一般一维和二维数组.

1
2
typedef Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime, Options> MyMatrixType;
typedef Array<Scalar, RowsAtCompileTime, ColsAtCompileTime, Options> MyArrayType;
  • Scalar 为标量类型 (如 float, double, bool, int, 等).
  • RowsAtCompileTime 和 ColsAtCompileTime 为编译时确定的矩阵行数和列数, 或者为 Dynamic.
  • Options 可以是 ColMajor 或 RowMajor, 默认为 ColMajor. (更多选项参见类 Matrix)

所有的组合都是允许的: 你可以有一个固定行数和动态列数的矩阵, 等等. 以下内容均有效:

1
2
3
4
Matrix<double, 6, Dynamic>                  // 动态列数 (堆分配)
Matrix<double, Dynamic, 2>                  // 动态行数 (堆分配)
Matrix<double, Dynamic, Dynamic, RowMajor>  // 完全动态, 行优先 (堆分配)
Matrix<double, 13, 3>                       // 完全固定 (通常在栈上分配)

在大多数情况下, 你可以简单地对 matrices 和 arrays 使用一些方便的预定义类型. 例如:

Matrices

1
2
3
4
5
Matrix<float,Dynamic,Dynamic>   <=>   MatrixXf
Matrix<double,Dynamic,1>        <=>   VectorXd
Matrix<int,1,Dynamic>           <=>   RowVectorXi
Matrix<float,3,3>               <=>   Matrix3f
Matrix<float,4,1>               <=>   Vector4f

Arrays

1
2
3
4
5
Array<float,Dynamic,Dynamic>    <=>   ArrayXXf
Array<double,Dynamic,1>         <=>   ArrayXd
Array<int,1,Dynamic>            <=>   RowArrayXi
Array<float,3,3>                <=>   Array33f
Array<float,4,1>                <=>   Array4f

矩阵和数组之间的转换:

1
2
3
4
5
6
7
8
Array44f a1, a2;
Matrix4f m1, m2;
m1 = a1 * a2;                     // 元素乘积, 从数组到矩阵的隐式转换.
a1 = m1 * m2;                     // 矩阵乘积, 从矩阵到数组的隐式转换.
a2 = a1 + m1.array();             // 数组和矩阵的混合是被禁止的,
m2 = a1.matrix() + m1;            // 需要显式的转换.
ArrayWrapper<Matrix4f> m1a(m1);   // m1a 是 m1.array() 的一个别名, 它们共享相同的元素
MatrixWrapper<Array44f> a1m(a1);

在本文档的剩余部分, 我们将使用以下符号来强调特定于给定类型对象的特征:

  • * 仅适用与线性代数矩阵与向量
  • ** 仅适用于数组对象

基本矩阵操作

构造

1
2
3
4
5
6
7
8
9
10
11
12
/* 一维对象 */
Vector4d  v4;
Vector2f  v1(x, y);
Array3i   v2(x, y, z);
Vector4d  v3(x, y, z, w);
VectorXf  v5; // 空对象
ArrayXf   v6(size);
/* 二维对象 */
Matrix4f  m1; 
MatrixXf  m5; // 空对象
MatrixXf  m6(nb_rows, nb_columns);
// 注: 默认情况下, 不初始化元素

逗号初始化

1
2
3
4
5
6
7
/* 一维对象 */
Vector3f  v1;     v1 << x, y, z;
ArrayXf   v2(4);  v2 << 1, 2, 3, 4;
/* 二维对象 */
Matrix3f  m1;   m1 << 1, 2, 3,
                      4, 5, 6,
                      7, 8, 9;

逗号初始化(bis)

1
2
3
4
5
6
7
int rows=5, cols=5;
MatrixXf m(rows,cols);
m << (Matrix3f() << 1, 2, 3, 4, 5, 6, 7, 8, 9).finished(),
     MatrixXf::Zero(3,cols-3),
     MatrixXf::Zero(rows-3,3),
     MatrixXf::Identity(rows-3,cols-3);
cout << m

输出:
1 2 3 0 0
4 5 6 0 0
7 8 9 0 0
0 0 0 1 0
0 0 0 0 1

运行时信息

1
2
3
4
5
6
7
8
9
10
/* 一维对象 */
vector.size();
vector.innerStride();
vector.data();
/* 二维对象 */
matrix.rows();          matrix.cols();
matrix.innerSize();     matrix.outerSize();
matrix.innerStride();   matrix.outerStride();
matrix.data();
// Inner/Outer* 依赖于存储顺序

编译时信息

1
2
3
ObjectType::Scalar              ObjectType::RowsAtCompileTime
ObjectType::RealScalar          ObjectType::ColsAtCompileTime
ObjectType::Index               ObjectType::SizeAtCompileTime

改变规格

1
2
3
4
5
6
7
8
9
10
11
12
/* 一维对象 */
vector.resize(size);
vector.resizeLike(other_vector);
vector.conservativeResize(size);
/* 二维对象 */
matrix.resize(nb_rows, nb_cols);
matrix.resize(Eigen::NoChange, nb_cols);
matrix.resize(nb_rows, Eigen::NoChange);
matrix.resizeLike(other_matrix);
matrix.conservativeResize(nb_rows, nb_cols);
// 如果新规格匹配, 则不执行操作, 否则数据将丢失
// 通过数据保存调整规格

带范围检查的元素访问

1
2
3
4
5
6
7
8
/* 一维对象 */
vector(i)     vector.x()
vector[i]     vector.y()
              vector.z()
              vector.w()
/* 二维对象 */
matrix(i,j)
// 通过定义 NDEBUG 或 EIGEN_NO_DEBUG, 禁用范围检查

没有范围检查的元素访问

1
2
3
4
5
6
/* 一维对象 */
vector.coeff(i)
vector.coeffRef(i)
/* 二维对象 */
matrix.coeff(i,j)
matrix.coeffRef(i,j)

赋值/拷贝

1
2
3
object = expression;
object_of_float = expression_of_double.cast<float>();
// 目标会自动调整规格(如果可能的话)

预定义矩阵

固定规格的矩阵或向量

1
2
3
4
5
6
7
8
9
10
11
12
13
14
typedef {Matrix3f|Array33f} FixedXD;
FixedXD x;
 
x = FixedXD::Zero();
x = FixedXD::Ones();
x = FixedXD::Constant(value);
x = FixedXD::Random();
x = FixedXD::LinSpaced(size, low, high);
 
x.setZero();
x.setOnes();
x.setConstant(value);
x.setRandom();
x.setLinSpaced(size, low, high);

单位和基向量 *

1
2
3
4
5
6
7
8
x = FixedXD::Identity();
x.setIdentity();
 
Vector3f::UnitX() // 1 0 0
Vector3f::UnitY() // 0 1 0
Vector3f::UnitZ() // 0 0 1
Vector4f::Unit(i)
x.setUnit(i);

动态规格的矩阵

1
2
3
4
5
6
7
8
9
10
11
12
13
14
typedef {MatrixXf|ArrayXXf} Dynamic2D;
Dynamic2D x;
 
x = Dynamic2D::Zero(rows, cols);
x = Dynamic2D::Ones(rows, cols);
x = Dynamic2D::Constant(rows, cols, value);
x = Dynamic2D::Random(rows, cols);
N/A
 
x.setZero(rows, cols);
x.setOnes(rows, cols);
x.setConstant(rows, cols, value);
x.setRandom(rows, cols);
N/A

单位和基向量 *

1
2
3
4
x = Dynamic2D::Identity(rows, cols);
x.setIdentity(rows, cols);

N/A

动态规格的向量

1
2
3
4
5
6
7
8
9
10
11
12
13
14
typedef {VectorXf|ArrayXf} Dynamic1D;
Dynamic1D x;
 
x = Dynamic1D::Zero(size);
x = Dynamic1D::Ones(size);
x = Dynamic1D::Constant(size, value);
x = Dynamic1D::Random(size);
x = Dynamic1D::LinSpaced(size, low, high);
 
x.setZero(size);
x.setOnes(size);
x.setConstant(size, value);
x.setRandom(size);
x.setLinSpaced(size, low, high);

单位和基向量 *

1
2
3
4
5
6
N/A

VectorXf::Unit(size,i)
x.setUnit(size,i);
VectorXf::Unit(4,1) == Vector4f(0,1,0,0)
                    == Vector4f::UnitY()

请注意, 允许对动态规格的向量或矩阵调用任意 set* 函数, 而无需传递新的规格. 例如:

1
2
MatrixXi M(3,3);
M.setIdentity();

映射外部数组

连续内存

1
2
3
4
5
float data[] = {1,2,3,4};
Map<Vector3f> v1(data);       // uses v1 as a Vector3f object
Map<ArrayXf>  v2(data,3);     // uses v2 as a ArrayXf object
Map<Array22f> m1(data);       // uses m1 as a Array22f object
Map<MatrixXf> m2(data,2,2);   // uses m2 as a MatrixXf object

Strides 的典型用法

1
2
3
4
5
float data[] = {1,2,3,4,5,6,7,8,9};
Map<VectorXf,0,InnerStride<2> >  v1(data,3);                      // = [1,3,5]
Map<VectorXf,0,InnerStride<> >   v2(data,3,InnerStride<>(3));     // = [1,4,7]
Map<MatrixXf,0,OuterStride<3> >  m2(data,2,3);                    // both lines     |1,4,7|
Map<MatrixXf,0,OuterStride<> >   m1(data,2,3,OuterStride<>(3));   // are equal to:  |2,5,8|

算术运算

加减

1
2
mat3 = mat1 + mat2;           mat3 += mat1;
mat3 = mat1 - mat2;           mat3 -= mat1;

数乘

1
2
mat3 = mat1 * s1;             mat3 *= s1;           mat3 = s1 * mat1;
mat3 = mat1 / s1;             mat3 /= s1;

矩阵/向量乘法*

1
2
3
col2 = mat1 * col1;
row2 = row1 * mat1;           row1 *= mat1;
mat3 = mat1 * mat2;           mat3 *= mat1;

转置 伴随* 

1
2
mat1 = mat2.transpose();      mat1.transposeInPlace();
mat1 = mat2.adjoint();        mat1.adjointInPlace();

点乘 内积*

1
2
3
scalar = vec1.dot(vec2);
scalar = col1.adjoint() * col2;
scalar = (col1.adjoint() * col2).value();

外积* 

1
mat = col1 * col2.transpose();

范数 归一化* 

1
2
scalar = vec1.norm();         scalar = vec1.squaredNorm()
vec2 = vec1.normalized();     vec1.normalize(); // inplace

叉乘*

1
2
#include <Eigen/Geometry>
vec3 = vec1.cross(vec2);

元素优先的数组运算

除了上述运算符外, Eigen 还支持许多元素运算符和函数. 它们中的大多数在数组世界**中毫无疑问是有意义的. 以下操作符可以用于数组, 也可以通过 .array() 用于向量和矩阵:

算术运算符

1
2
array1 * array2     array1 / array2     array1 *= array2    array1 /= array2
array1 + scalar     array1 - scalar     array1 += scalar    array1 -= scalar

比较

1
2
3
4
array1 < array2     array1 > array2     array1 < scalar     array1 > scalar
array1 <= array2    array1 >= array2    array1 <= scalar    array1 >= scalar
array1 == array2    array1 != array2    array1 == scalar    array1 != scalar
array1.min(array2)  array1.max(array2)  array1.min(scalar)  array1.max(scalar)

三角函数, 幂函数和其它函数以及类 STL 的变体

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
array1.abs2()
array1.abs()                  abs(array1)
array1.sqrt()                 sqrt(array1)
array1.log()                  log(array1)
array1.log10()                log10(array1)
array1.exp()                  exp(array1)
array1.pow(array2)            pow(array1,array2)
array1.pow(scalar)            pow(array1,scalar)
                              pow(scalar,array2)
array1.square()
array1.cube()
array1.inverse()
 
array1.sin()                  sin(array1)
array1.cos()                  cos(array1)
array1.tan()                  tan(array1)
array1.asin()                 asin(array1)
array1.acos()                 acos(array1)
array1.atan()                 atan(array1)
array1.sinh()                 sinh(array1)
array1.cosh()                 cosh(array1)
array1.tanh()                 tanh(array1)
array1.arg()                  arg(array1)
 
array1.floor()                floor(array1)
array1.ceil()                 ceil(array1)
array1.round()                round(aray1)
 
array1.isFinite()             isfinite(array1)
array1.isInf()                isinf(array1)
array1.isNaN()                isnan(array1)

以下元素操作符可用于所有类型的表达式(矩阵, 向量和数组), 以及实数或复数标量类型:

Eigen API ------------------- 类STL API**

1
2
3
mat1.real()          real(array1)    // 读写, 对实表达式无操作
mat1.imag()          imag(array1)    // 对实表达式只读, 对复表达式读写
mat1.conjugate()     conj(array1)    // 对实表达式无操作

通过以下 cwise* 方法, 可以很容易地获得一些矩阵和向量的元素优先运算符:

Matrix API *

1
2
3
4
5
6
7
8
9
10
mat1.cwiseMin(mat2)         mat1.cwiseMin(scalar)
mat1.cwiseMax(mat2)         mat1.cwiseMax(scalar)
mat1.cwiseAbs2()
mat1.cwiseAbs()
mat1.cwiseSqrt()
mat1.cwiseInverse()
mat1.cwiseProduct(mat2)
mat1.cwiseQuotient(mat2)
mat1.cwiseEqual(mat2)       mat1.cwiseEqual(scalar)
mat1.cwiseNotEqual(mat2)

通过 Array 转换

1
2
3
4
5
6
7
8
9
10
mat1.array().min(mat2.array())    mat1.array().min(scalar)
mat1.array().max(mat2.array())    mat1.array().max(scalar)
mat1.array().abs2()
mat1.array().abs()
mat1.array().sqrt()
mat1.array().inverse()
mat1.array() * mat2.array()
mat1.array() / mat2.array()
mat1.array() == mat2.array()      mat1.array() == scalar
mat1.array() != mat2.array()

这两个 API 之间的主要区别在于, 基于 cwise* 方法的 API 返回矩阵世界中的表达式, 而第二个 API (基于.array()) 返回数组表达式. 回顾一下, .array() 没有成本, 它只改变可用的 API 和数据的解释.

使用 DenseBase::unaryExpr 和 std::ptr_fun (c++03, 在较新的c++版本中已弃用或删除), std::ref (c++11) 或 lambdas (c++11) 来应用任何用户定义的函数 foo 也非常简单:

1
2
3
mat1.unaryExpr(std::ptr_fun(foo));
mat1.unaryExpr(std::ref(foo));
mat1.unaryExpr([](double x) { return foo(x); });

请注意, 不可能将原始函数指针传递给 unaryExpr, 因此请像上面所示的那样进行转换.

简化

Eigen 提供了一些简化方法, 例如: minCoeff(), maxCoeff(), sum(), prod(), trace() , norm() , squaredNorm() *, all() , 和 any() . 所有的简化算符都可以按照 矩阵优先, 列优先 或 行优先 方式进行. 用法示例:

1
2
3
4
5
6
//       5 3 1
// mat = 2 7 8
//       9 4 6
mat.minCoeff();              // 1
mat.colwise().minCoeff();    // 2 3 1
mat.rowwise().minCoeff();    // [1 2 4]'

minCoeff 和 maxCoeff 的特殊版本:

1
2
3
int i, j;
s = vector.minCoeff(&i);        // s == vector[i]
s = matrix.maxCoeff(&i, &j);    // s == matrix(i,j)

all() 和 any() 的典型用例:

1
2
if((array1 > 0).all()) ...      // 如果 array1 的所有参数都大于 0 ...
if((array1 < array2).any()) ... // 如果存在一对 i,j 使得 array1(i,j) < array2(i,j) ...

子矩阵

Eigen 3.4 为子矩阵提供了一个改进的 API, 包括对数组的切片和索引:

对矩阵(或数组)的列或行进行读写访问:

1
2
mat1.row(i) = mat2.col(j);
mat1.col(j1).swap(mat1.col(j2));

对子向量的读写访问:

默认版本 --------------------- 当规格在编译时已知的优化版本

1
2
3
vec1.head(n)         vec1.head<n>()        // 前 n 个元素
vec1.tail(n)         vec1.tail<n>()        // 后 n 个元素
vec1.segment(pos,n)  vec1.segment<n>(pos)  // 范围 [pos:pos+n-1] 的n个元素

对子矩阵的读写访问:

默认版本 --------------------- 当规格在编译时已知的优化版本

1
2
3
4
5
6
7
8
9
10
11
12
// 从位置 (i,j) 开始的 rows x cols 子矩阵
mat1.block(i,j,rows,cols)          mat1.block<rows,cols>(i,j)
// 四个角落之一的 rows x cols 子矩阵
mat1.topLeftCorner(rows,cols)      mat1.topLeftCorner<rows,cols>()
mat1.topRightCorner(rows,cols)     mat1.topRightCorner<rows,cols>()
mat1.bottomLeftCorner(rows,cols)   mat1.bottomLeftCorner<rows,cols>()
mat1.bottomRightCorner(rows,cols)  mat1.bottomRightCorner<rows,cols>()
// 当分块占据两个角落时的特殊版本 block()
mat1.topRows(rows)                 mat1.topRows<rows>()
mat1.bottomRows(rows)              mat1.bottomRows<rows>()
mat1.leftCols(cols)                mat1.leftCols<cols>()
mat1.rightCols(cols)               mat1.rightCols<cols>()

其它运算

Eigen 3.4 支持一个新的 Reshape API.

反转

向量, 矩阵的行与列都可以进行反转 (参见 DenseBase::reverse(), DenseBase::reverseInPlace(), VectorwiseOp::reverse()).

1
2
vec.reverse()           mat.colwise().reverse()   mat.rowwise().reverse()
vec.reverseInPlace()

复制

向量, 矩阵的行与列可以在任何方向上进行复制 (参见 DenseBase::replicate(), VectorwiseOp::replicate()).

1
2
3
4
vec.replicate(times)                                          vec.replicate<Times>
mat.replicate(vertical_times, horizontal_times)               mat.replicate<VerticalTimes, HorizontalTimes>()
mat.colwise().replicate(vertical_times, horizontal_times)     mat.colwise().replicate<VerticalTimes, HorizontalTimes>()
mat.rowwise().replicate(vertical_times, horizontal_times)     mat.rowwise().replicate<VerticalTimes, HorizontalTimes>()

对角, 三角, 自伴矩阵

(矩阵世界*)

对角矩阵

通过向量构造对角矩阵

1
mat1 = vec1.asDiagonal();

声明一个对角矩阵

1
2
DiagonalMatrix<Scalar,SizeAtCompileTime> diag1(size);
diag1.diagonal() = vector;

作为向量访问矩阵的对角线和上/下对角 (读/写)

1
2
3
4
5
vec1 = mat1.diagonal();        mat1.diagonal() = vec1;      // main diagonal
vec1 = mat1.diagonal(+n);      mat1.diagonal(+n) = vec1;    // n-th super diagonal
vec1 = mat1.diagonal(-n);      mat1.diagonal(-n) = vec1;    // n-th sub diagonal
vec1 = mat1.diagonal<1>();     mat1.diagonal<1>() = vec1;   // first super diagonal
vec1 = mat1.diagonal<-2>();    mat1.diagonal<-2>() = vec1;  // second sub diagonal

优化的乘积和求逆

1
2
3
4
mat3  = scalar * diag1 * mat1;
mat3 += scalar * mat1 * vec1.asDiagonal();
mat3 = vec1.asDiagonal().inverse() * mat1
mat3 = mat1 * diag1.inverse()

三角矩阵

TriangularView 给出了密集矩阵三角部分的方法, 并允许对其执行优化操作. 相对的三角部分永远不会被引用, 可以用来存储其他信息.

注意: 如果 .triangularView() 模板成员函数用于依赖模板参数的类型的对象, 则需要 template 关键字; 详细信息请参见C++中的模板和类型名关键字.

引用带有可选单位或空对角线的三角矩阵 (读/写):

1
m.triangularView<Xxx>()

Xxx = Upper, Lower, StrictlyUpper, StrictlyLower, UnitUpper, UnitLower.

写入一个特定的三角部分 (只计算参考的三角部分):

1
m1.triangularView<Eigen::Lower>() = m2 + m3

转换为一个稠密矩阵, 使对三角形部分为零:

1
m2 = m1.triangularView<Eigen::UnitUpper>()

乘积:

1
2
m3 += s1 * m1.adjoint().triangularView<Eigen::UnitUpper>() * m2
m3 -= s1 * m2.conjugate() * m1.adjoint().triangularView<Eigen::Lower>()

求解线性方程:

, , 

1
2
3
L1.triangularView<Eigen::UnitLower>().solveInPlace(M2)
L1.triangularView<Eigen::Lower>().adjoint().solveInPlace(M3)
U1.triangularView<Eigen::Upper>().solveInPlace<OnTheRight>(M4)

对称/自伴矩阵

就像三角矩阵一样, 你可以参考一个方阵的任何三角部分, 将其视为一个自伴矩阵, 并进行特殊的优化操作. 同样, 相对的三角部分永远不会被引用, 可以用来存储其他信息.

注意: 如果 .selfadjointView() 模板成员函数用于依赖模板参数的类型的对象, 则需要 template 关键字; 详细信息请参见C++中的模板和类型名关键字.

转换为密集矩阵:

1
m2 = m.selfadjointView<Eigen::Lower>();

与另一个一般矩阵或向量的乘积:

1
2
m3  = s1 * m1.conjugate().selfadjointView<Eigen::Upper>() * m3;
m3 -= s1 * m3.adjoint() * m1.selfadjointView<Eigen::Lower>();

Rank 1 和 Rank K 更新:

1
2
M1.selfadjointView<Eigen::Upper>().rankUpdate(M2,s1);
M1.selfadjointView<Eigen::Lower>().rankUpdate(M2.adjoint(),-1);

Rank 2 更新:

1
M.selfadjointView<Eigen::Upper>().rankUpdate(u,v,s);

求解线性方程:

1
2
3
4
// via a standard Cholesky factorization
m2 = m1.selfadjointView<Eigen::Upper>().llt().solve(m2);
// via a Cholesky factorization with pivoting
m2 = m1.selfadjointView<Eigen::Lower>().ldlt().solve(m2);