Eigen教程:1 Eigen简介和矩阵常见操作
Eigen是可以用来进行线性代数、矩阵、向量操作等运算的C++库，它里面包含了很多算法。它的License是MPL2。它支持多平台。
Eigen采用源码的方式提供给用户使用，在使用时只需要包含Eigen的头文件即可进行使用。之所以采用这种方式，是因为Eigen采用模板方式实现，由于模板函数不支持分离编译，所以只能提供源码而不是动态库的方式供用户使用。
一. 模块和头文件
Eigen库被分为一个Core模块和其他一些模块，每个模块有一些相应的头文件。 为了便于引用，Dense模块整合了一系列模块；Eigen模块整合了所有模块。一般情况下，#include<Eigen/Dense> 就够了。
| Module | Header file | Contents |
|---|---|---|
| Core | #include<Eigen/Core> | Matrix和Array类，基础的线性代数运算和数组操作 |
| Geometry | #include<Eigen/Geometry> | 旋转、平移、缩放、2维和3维的各种变换 |
| LU | #include<Eigen/LU> | 求逆，行列式，LU分解 |
| Cholesky | #include <Eigen/Cholesky> | LLT和LDLT Cholesky分解 |
| Householder | #include<Eigen/Householder> | 豪斯霍尔德变换，用于线性代数运算 |
| SVD | #include<Eigen/SVD> | SVD分解 |
| QR | #include<Eigen/QR> | QR分解 |
| Eigenvalues | #include<Eigen/Eigenvalues> | 特征值，特征向量分解 |
| Sparse | #include<Eigen/Sparse> | 稀疏矩阵的存储和一些基本的线性运算 |
| 稠密矩阵 | #include<Eigen/Dense> | 包含了Core/Geometry/LU/Cholesky/SVD/QR/Eigenvalues模块 |
| 矩阵 | #include<Eigen/Eigen> | 包括Dense和Sparse(整合库) |
二. 矩阵
Eigen中关于矩阵类的模板函数中，共有六个模板参数，常用的只有前三个。其前三个参数分别表示矩阵元素的类型、行数和列数。
Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>
Scalar是表示元素的类型，RowsAtCompileTime为矩阵的行，ColsAtCompileTime为矩阵的列。
其他模板参数
Matrix<typename Scalar,
int RowsAtCompileTime,
int ColsAtCompileTime,
int Options = 0,
int MaxRowsAtCompileTime = RowsAtCompileTime,
int MaxColsAtCompileTime = ColsAtCompileTime>
Options是一个比特标志位，这里，我们只介绍一种RowMajor，它表明matrix使用按行存储，默认是按列存储。Matrix<float, 3, 3, RowMajor>
MaxRowsAtCompileTime和MaxColsAtCompileTime表示在编译阶段矩阵的上限。主要是避免动态内存分配，使用数组。例如：Matrix<float, Dynamic, Dynamic, 0, 3, 4> 等价于 float [12];
矩阵定义时可以使用Dynamic来表示矩阵的行列数为未知。
Eigen中无论是矩阵还是数组、向量，无论是静态矩阵还是动态矩阵都提供默认构造函数，也就是定义这些数据结构时都可以不用提供任何参数，其大小均由运行时来确定。矩阵的构造函数中只提供行列数、元素类型的构造参数，而不提供元素值的构造，对于比较小的、固定长度的向量提供初始化元素的定义。
1. 矩阵类型
Eigen中的矩阵类型一般都是用类似MatrixXXX来表示，可以根据该名字来判断其数据类型，比如”d”表示double类型，”f”表示float类型，”i”表示整数，”c”表示复数；Matrix2f，表示的是一个2*2维的，其每个元素都是float类型。
2. 数据存储
Matrix创建的矩阵默认是按列存储，Eigen在处理按列存储的矩阵时会更加高效。如果想修改可以在创建矩阵的时候加入参数，如:
Matrix<int,3, 4, ColMajor> Acolmajor;
Matrix<int,3, 4, RowMajor> Arowmajor;
3. 动态矩阵和静态矩阵
动态矩阵是指其大小在运行时确定，静态矩阵是指其大小在编译时确定。
MatrixXd：表示任意大小的元素类型为double的矩阵变量，其大小只有在运行时被赋值之后才能知道。
Matrix3d：表示元素类型为double大小为3*3的矩阵变量，其大小在编译时就知道。
在Eigen中行优先的矩阵会在其名字中包含有row，否则就是列优先。
Eigen中的向量只是一个特殊的矩阵，其维度为1而已。
4. 矩阵元素的访问
在矩阵的访问中，行索引总是作为第一个参数，Eigen中矩阵、数组、向量的下标都是从0开始。矩阵元素的访问可以通过”()”操作符完成。例如m(2, 3)既是获取矩阵m的第2行第3列元素。
针对向量还提供”[]”操作符，注意矩阵则不可如此使用。
5. 设置矩阵的元素
在Eigen中重载了”<<”操作符，通过该操作符即可以一个一个元素的进行赋值，也可以一块一块的赋值。另外也可以使用下标进行赋值。
6. 重置矩阵大小
当前矩阵的行数、列数、大小可以通过rows()、cols()和size()来获取，对于动态矩阵可以通过resize()函数来动态修改矩阵的大小。注意：
- 固定大小的矩阵是不能使用resize()来修改矩阵的大小；
- resize()函数会析构掉原来的数据，因此调用resize()函数之后将不能保证元素的值不改变；
- 使用”=”操作符操作动态矩阵时，如果左右两边的矩阵大小不等，则左边的动态矩阵的大小会被修改为右边的大小。
- 如何选择动态矩阵和静态矩阵：对于小矩阵(一般大小小于16)使用固定大小的静态矩阵，它可以带来比较高的效率；对于大矩阵(一般大小大于32)建议使用动态矩阵。注意：如果特别大的矩阵使用了固定大小的静态矩阵则可能会造成栈溢出的问题。
7. 矩阵和向量的算术运算
在Eigen中算术运算重载了C++的+、-、*
-
矩阵的运算：提供+、-、一元操作符”-”、+=、-=；二元操作符+/-，表示两矩阵相加(矩阵中对应元素相加/减，返回一个临时矩阵)；一元操作符-表示对矩阵取负(矩阵中对应元素取负，返回一个临时矩阵)；组合操作法+=或者-=表示(对应每个元素都做相应操作)；矩阵还提供与标量(单一数字)的乘除操作，表示每个元素都与该标量进行乘除操作；
-
求矩阵的转置、共轭矩阵、伴随矩阵：可以通过成员函数transpose()、conjugate()、adjoint()来完成。注意：这些函数返回操作后的结果，而不会对原矩阵的元素进行直接操作，如果要让原矩阵进行转换，则需要使用响应的InPlace函数，如transpoceInPlace()等；
-
矩阵相乘、矩阵向量相乘：使用操作符 ∗ * ∗，共有 ∗ * ∗和 ∗ = *= ∗=两种操作符；
-
矩阵的块操作：有两种使用方法：
4.1matrix.block(i,j, p, q): 表示返回从矩阵(i, j)开始，每行取p个元素，每列取q个元素所组成的临时新矩阵对象，原矩阵的元素不变；4.2
matrix.block<p,q>(i, j):<p, q>可理解为一个p行q列的子矩阵，该定义表示从原矩阵中第(i, j)开始，获取一个p行q列的子矩阵，返回该子矩阵组成的临时矩阵对象，原矩阵的元素不变；
| operation | 构建一个动态尺寸的block | 构建一个固定尺寸的block |
|---|---|---|
| 起点(i,j)块大小(p,q) | .block(i,j,p,q) | .block< p,q >(i,j) |
两个版本都可以用于固定尺寸和动态尺寸的matrix/array。功能是等价的，只是固定尺寸的版本在block较小时速度更快一些。
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 in the middle
6 7
10 11
Block of size 1x1
1
Block of size 2x2
1 2
5 6
Block of size 3x3
1 2 3
5 6 7
9 10 11
8. 行和列
| Operation | Method |
|---|---|
| i t h i^{th} ith row | .matrix.row(i) |
| i t h i^{th} ithcolum | .matrix.col(j) |
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;
}
输出
Here is the matrix m:
1 2 3
4 5 6
7 8 9
2nd Row: 4 5 6
After adding 3 times the first column into the third column, the matrix m is:
1 2 6
4 5 18
7 8 30
9. 角相关操作
| Module | dynamic-size block | fixed-size block |
|---|---|---|
| 左上角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§; | matrix.leftCols< p >(); |
| 右p列 | matrix.rightCols§; | matrix.rightCols< p >(); |
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;
}
输出
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
10. 归约操作
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;
}
输出
Here is mat.sum(): 10
Here is mat.prod(): 24
Here is mat.mean(): 2.5
Here is mat.minCoeff(): 1
Here is mat.maxCoeff(): 4
Here is mat.trace(): 5
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;
输出
Here is the matrix m:
0.68 0.597 -0.33
-0.211 0.823 0.536
0.566 -0.605 -0.444
Its minimum coefficient (-0.605) is at position (2,1)
Here is the vector v: 1 0 3 -3
Its maximum coefficient (3) is at position 2
11. 操作的有效性
Eigen会检测执行操作的有效性，在编译阶段Eigen会检测它们，错误信息是繁冗的，但错误信息会大写字母突出，比如:
Matrix3f m;
Vector4f v;
v = m*v; // Compile-time error: YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES
当然动态尺寸的错误要在运行时发现，如果在debug模式，assertions会触发后，程序将崩溃。
MatrixXf m(3,3);
VectorXf v(4);
v = m * v; // Run-time assertion failure here: "invalid matrix product"
三. Map类
Map类：在已经存在的矩阵或向量中，不必拷贝对象，而是直接在该对象的内存上进行运算操作；通常Eigen与原生raw C/C++ 数组混合编程时候的选择；
1. Map的定义
Map<Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime> >
默认情况下，Map只需要一个模板参数。
为了构建Map变量，我们需要其余的两个信息
- 一个指向元素数组的指针
- Matrix/vector的尺寸
定义一个float类型的矩阵： Map<MatrixXf> mf(pf,rows,columns); pf是一个数组指针float *。
固定尺寸的整形vector声明： Map<const Vector4i> mi(pi);
注意:Map没有默认的构造函数，你需要传递一个指针来初始化对象。
Mat是灵活地足够去容纳多种不同的数据表示，其他的两个模板参数：
Map<typename MatrixType,
int MapOptions,
typename StrideType>
MapOptions标识指针是否是对齐的（Aligned），默认是Unaligned。StrideType表示内存数组的组织方式：行列的步长。
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;
cout << "Row-major using stride:\n" <<
Map<Matrix<int,2,4>, Unaligned, Stride<1,4> >(array) << endl;
输出
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
2. Map变量和Eigen变量的混用
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); // get the address storing the data for m2
MapType m2map(p,m2.size()); // m2map shares data with m2
MapTypeConst m2mapconst(p,m2.size()); // a read-only accessor for 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; // this will change m2, since they share the same array
cout << "Updated m2: " << m2 << endl;
cout << "m2 coefficient 2, constant accessor: " << m2mapconst(2) << endl;
/* m2mapconst(2) = 5; */ // this yields a compile-time error
输出
m1: 0.68 -0.211 0.566 0.597 0.823
m2: -0.605 -0.33 0.536 -0.444 0.108
Squared euclidean distance: 3.26
Squared euclidean distance, using map: 3.26
Updated m2: -0.605 -0.33 0.536 7 0.108
m2 coefficient 2, constant accessor: 0.536
Eigen提供的函数都兼容Map对象。
3. 改变mapped数组
Map对象声明后，可以通过C++的placement new语法来改变Map的数组。
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";
The mapped vector v is: 1 2 3 4
Now v is: 5 6 7 8 9
四. 向量
1. 向量的点运算和叉运算
dot()执行点积，cross()执行叉积，点运算得到1*1的矩阵。当然，点运算也可以用u.adjoint()*v来代替。
Eigen::Vector3d v(1,2,3);
Eigen::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;
输出结果是
Dot product: 8
Dot product via a matrix product: 8
Cross product:
1
-2
1
叉乘结果是一个向量，符合右手定则。
注意：
点积只对三维vector有效。对于复数，Eigen的点积是第一个变量共轭和第二个变量的线性积。
2. 向量的块操作
| operation | dynamic-size block | fixed-size block |
|---|---|---|
| 前n个 | vector.head(n); | vector.head< n >(); |
| 后n个 | vector.tail(n); | vector.tail< n >(); |
| i起始的n个元素 | vector.segment(i,n); | vector.segment< n >(i); |
Eigen 示例代码
- 从Eigen 3.2.5下载最新稳定版本3.2.5，解压缩；
- 新建一个
vs2013 控制台工程，将Eigen文件所在目录加入到工程属性的C/C++附加包含目录中，这样就可以使用Eigen中的函数了； - Eigen 附带的测试代码：
estEigen.cpp文件中的内容为
#include "stdafx.h"
#include <iostream>
#include <Eigen/Dense>
template <typename T>
static void matrix_mul_matrix(T* p1, int iRow1, int iCol1, T* p2, int iRow2, int iCol2, T* p3)
{
if (iRow1 != iRow2) return;
//列优先
//Eigen::Map< Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > map1(p1, iRow1, iCol1);
//Eigen::Map< Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > map2(p2, iRow2, iCol2);
//Eigen::Map< Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > map3(p3, iCol1, iCol2);
//行优先
Eigen::Map< Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > map1(p1, iRow1, iCol1);
Eigen::Map< Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > map2(p2, iRow2, iCol2);
Eigen::Map< Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > map3(p3, iCol1, iCol2);
map3 = map1 * map2;
}
int main(int argc, char* argv[])
{
//1. 矩阵的定义
Eigen::MatrixXd m(2, 2);
Eigen::Vector3d vec3d;
Eigen::Vector4d vec4d(1.0, 2.0, 3.0, 4.0);
//2. 动态矩阵、静态矩阵
Eigen::MatrixXd matrixXd;
Eigen::Matrix3d matrix3d;
//3. 矩阵元素的访问
m(0, 0) = 1;
m(0, 1) = 2;
m(1, 0) = m(0, 0) + 3;
m(1, 1) = m(0, 0) * m(0, 1);
std::cout << m << std::endl << std::endl;
//4. 设置矩阵的元素
m << -1.5, 2.4,
6.7, 2.0;
std::cout << m << std::endl << std::endl;
int row = 4;
int col = 5;
Eigen::MatrixXf matrixXf(row, col);
matrixXf << 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20;
std::cout << matrixXf << std::endl << std::endl;
matrixXf << Eigen::MatrixXf::Identity(row, col);
std::cout << matrixXf << std::endl << std::endl;
//5. 重置矩阵大小
Eigen::MatrixXd matrixXd1(3, 3);
m = matrixXd1;
std::cout << m.rows() << " " << m.cols() << std::endl << std::endl;
//6. 矩阵运算
m << 1, 2, 7,
3, 4, 8,
5, 6, 9;
std::cout << m << std::endl;
matrixXd1 = Eigen::Matrix3d::Random();
m += matrixXd1;
std::cout << m << std::endl << std::endl;
m *= 2;
std::cout << m << std::endl << std::endl;
std::cout << -m << std::endl << std::endl;
std::cout << m << std::endl << std::endl;
//7. 求矩阵的转置、共轭矩阵、伴随矩阵
std::cout << m.transpose() << std::endl << std::endl;
std::cout << m.conjugate() << std::endl << std::endl;
std::cout << m.adjoint() << std::endl << std::endl;
std::cout << m << std::endl << std::endl;
m.transposeInPlace();
std::cout << m << std::endl << std::endl;
//8. 矩阵相乘、矩阵向量相乘
std::cout << m*m << std::endl << std::endl;
vec3d = Eigen::Vector3d(1, 2, 3);
std::cout << m * vec3d << std::endl << std::endl;
std::cout << vec3d.transpose()*m << std::endl << std::endl;
//9. 矩阵的块操作
std::cout << m << std::endl << std::endl;
std::cout << m.block(1, 1, 2, 2) << std::endl << std::endl;
std::cout << m.block<1, 2>(0, 0) << std::endl << std::endl;
std::cout << m.col(1) << std::endl << std::endl;
std::cout << m.row(0) << std::endl << std::endl;
//10. 向量的块操作
Eigen::ArrayXf arrayXf(10);
arrayXf << 1, 2, 3, 4, 5, 6, 7, 8, 9, 10;
std::cout << vec3d << std::endl << std::endl;
std::cout << arrayXf << std::endl << std::endl;
std::cout << arrayXf.head(5) << std::endl << std::endl;
std::cout << arrayXf.tail(4) * 2 << std::endl << std::endl;
//11. 求解矩阵的特征值和特征向量
Eigen::Matrix2f matrix2f;
matrix2f << 1, 2, 3, 4;
Eigen::SelfAdjointEigenSolver<Eigen::Matrix2f> eigenSolver(matrix2f);
if (eigenSolver.info() == Eigen::Success) {
std::cout << eigenSolver.eigenvalues() << std::endl << std::endl;
std::cout << eigenSolver.eigenvectors() << std::endl << std::endl;
}
//12. 类Map及动态矩阵的使用
int array1[4] = { 1, 2, 3, 4 };
int array2[4] = { 5, 6, 7, 8 };
int array3[4] = { 0, 0, 0, 0};
matrix_mul_matrix(array1, 2, 2, array2, 2, 2, array3);
for (int i = 0; i < 4; i++)
std::cout << array3[i] << std::endl;
return 0;
}
上面代码示例中涉及了转置矩阵 m T m^T mT,共轭矩阵 m ‾ \overline m m和伴随矩阵 m ∗ m^* m∗,对于实数矩阵，conjugate不执行任何操作，adjoint等价于transpose。
参考
- 测试代码Github地址
- 官方Eigen库的教程
- fengbingchun Eigen/OpenBLAS