Eigen: C++开源矩阵计算库
2020-12-13 04:33
标签:lse nta efault 常用 initial abs cts asi require Eigen库被分为一个Core模块和几个附加的模块,每个模块有一个相关的头文件,使用该模块时需要包含该头文件,为了能便利的使用eigen的几个模块,Eigen提供了Dense和Eigen两个头文件,各个头文件和模块如下表 Eigen库分为核心模块和额外模块两个部分, Eigen和Dense头文件方便的同时包含了几个头文件可以供使用 --Core 有关矩阵和数组的类,由基本的线性代数(包含 三角形和自伴乘积相关),还有相应对数组的操作。 --Geometry 几何学的类,有关转换、平移、进位制、2d旋转、3D旋转(四元组和角轴相关) --LU 逻辑单元的类,有关求逆,求行列式,LU分解解算器(FullPivLU,PartialPivLU) --Cholesky 包含LLT和LDLT的Cholesky因式分解法(Cholesky分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解) --Householder Householder变换,这个模块供几何线性代数模块使用 --SVD 奇异值分解,最小二乘解算器解决奇异值分解 --QR QR分解求解,三种方法:HouseholderQR、ColPivHouseholderQR,FullPivhouseholderQR --Eigenvalues 特征值和特征向量分解的方法:EigenSolver、SelfAdjointEigenSolver、ComplexEigenSolver --Sparse 稀疏矩阵相关类,对于系数矩阵的存储及相关线性代数 --Dense 包含:core、Geometry、LU、Cholesky、SVD、QR和Eigenvalues模块(头文件) --Eigen 包含上述所有模块(头文件) 1、矩阵的定义 Eigen中关于矩阵类的模板函数中,共有6个参数,但是常用的只有前三个,如下所示: 前三个参数分别表示矩阵元素的类型,行数和列数。 矩阵定义时可以使用Dynamic 来表示矩阵的行列数是未知,例如: 在Eigen中也提供了很多常见的简化定义形式,例如 注意: (1)Eigen中无论是矩阵还是数组、向量,无论是静态矩阵还是动态矩阵都提供默认构造函数,也就是你定义这些数据结构时都可以不用提供任何参数,其大小均由运行时来确定。 (2)矩阵的构造函数中只提供行列数、元素类型的构造参数,而不提供元素值的构造,对于比较小的、固定长度的向量提供初始化元素的定义,例如: 2、动态矩阵和静态矩阵 动态矩阵是指其大小在运行时确定,静态矩阵是指其大小在编译时确定,在Eigen中并未这样称呼矩阵。具体可见如下两段代码: 代码段1: 代码段2: 说明: 1)代码段1中MatrixXd表示任意大小的元素类型为double的矩阵变量,其大小只有在运行时被赋值之后才能知道; MatrixXd::Random(3,3)表示产生一个元素类型为double的3*3的临时矩阵对象。 2) 代码段2中Matrix3d表示元素类型为double大小为3*3的矩阵变量,其大小在编译时就知道; 3)上例中向量的定义也是类似,不过这里的向量时列优先,在Eigen中行优先的矩阵会在其名字中包含有row,否则就是列优先。 4)向量只是一个特殊的矩阵,其一个维度为1而已,如:typedef Matrix Vector3d 3、矩阵元素的访问 在矩阵的访问中,行索引总是作为第一个参数,需注意Eigen中遵循大家的习惯让矩阵、数组、向量的下标都是从0开始。矩阵元素的访问可以通过()操作符完成,例如m(2,3)即是获取矩阵m的第2行第3列元素(注意行列数从0开始)。可参看如下代码: 其输出结果为: 针对向量还提供[]操作符,注意矩阵则不可如此使用,原因为:在C++中m[i, j]中逗号表达式 “i, j”的值始终都是“j”的值,即m[i, j]对于C++来讲就是m[j]; 4、设置矩阵的元素 在Eigen中重载了"
代码段1: 输出结果为: 代码段二(使用下标进行复制): 5、重置矩阵大小 输出结果为: 6、如何选择动态矩阵和静态矩阵 Eigen对于这问题的答案是:对于小矩阵(一般大小小于16)的使用固定大小的静态矩阵,它可以带来比较高的效率,对于大矩阵(一般大小大于32)建议使用动态矩阵。 还需特别注意的是:如果特别大的矩阵使用了固定大小的静态矩阵则可能造成栈溢出的问题. 7、矩阵的运算 Eigen重载了基础的+ - * / += -= *= /= *可以表示标量和矩阵或者矩阵和矩阵 矩阵与矩阵的运算 Eigen提供+、-、一元操作符“-”、+=、-=,例如: 二元操作符+/-表示两矩阵相加(矩阵中对应元素相加/减,返回一个临时矩阵): B+C 或 B-C; 一元操作符-表示对矩阵取负(矩阵中对应元素取负,返回一个临时矩阵): -C; 组合操作法+=或者-=表示(对应每隔元素都做相应操作):A += B 或者 A-=B 代码段1为矩阵的加减操作,代码如下: 输出结果为: 矩阵与标量的运算 矩阵还提供与标量(单一个数字)的乘除操作,表示每个元素都与该标量进行乘除操作。例如: 二元操作符*在:A*a中表示矩阵A中的每个元素都与数字a相乘,结果放在一个临时矩阵中,矩阵的值不会改变。 对于a*A、A/a、A*=a、A /=a也是一样,例如下面的代码: 输出结果为: 需要注意: 在Eigen中,算术操作例如 “操作符+”并不会自己执行计算操作,他们只是返回一个“算术表达式对象”,而实际的计算则会延迟到后面的赋值时才进行。这些不影响你的使用,它只是为了方便Eigen的优化 8、求矩阵的转置、共轭矩阵、伴随矩阵、行列式、逆矩阵。 可以通过 成员函数transpose(), conjugate(),determinant(), adjoint()和inverse()来完成,注意这些函数返回操作后的结果,而不会对原矩阵的元素进行直接操作,如果要让原矩阵的进行转换,则需要使用响应的InPlace函数,例如:transposeInPlace() 、 adjointInPlace() 之类。 例如下面的代码所示: 输出结果为: 矩阵的特征值: m.eigenvalues(); 特征向量: m.eigenvector(); 9、矩阵相乘、矩阵向量相乘 矩阵的相乘,矩阵与向量的相乘也是使用操作符*,共有*和*=两种操作符,其用法可以参考如下代码: 输出结果为: 10、点积和叉积 11、矩阵的块操作 1)矩阵的块操作有两种使用方法,其定义形式为: (i,j); (2) 定义(1)表示返回从矩阵的(i, j)开始,每行取p个元素,每列取q个元素所组成的临时新矩阵对象,原矩阵的元素不变。 定义(2)中block(p, q)可理解为一个p行q列的子矩阵,该定义表示从原矩阵中第(i, j)开始,获取一个p行q列的子矩阵,返回该子矩阵组成的临时 矩阵对象,原矩阵的元素不变。 详细使用情况,可参考下面的代码段: 输出的结果为: 通过上述方式获取的子矩阵即可以作为左值也可以作为右值,也就是即可以用这个子矩阵给其他矩阵赋值,也可以给这个子矩阵对象赋值。 2)矩阵也提供了获取其指定行/列的函数,其实获取某行/列也是一种特殊的获取子块。可以通过 .col()和 .row()来完成获取指定列/行的操作,参数为列/行的索引。 输出结果为: 3)向量的块操作,其实向量只是一个特殊的矩阵,但是Eigen也为它单独提供了一些简化的块操作,如下三种形式: 输出结果为: 12、计算特征值和特征向量 13、矩阵分解 矩阵分解 (decomposition, factorization)是将矩阵拆解为数个矩阵的乘积,可分为三角分解、满秩分解、QR分解、Jordan分解和SVD(奇异值)分解等,常见的有三种:1)三角分解法 (Triangular Factorization),2)QR 分解法 (QR Factorization),3)奇异值分解法 (Singular Value Decompostion)。 Eigen总共提供了下面这些矩阵的分解方式: QR分解 QR分解法是将矩阵分解成一个正规正交矩阵与上三角形矩阵,所以称为QR分解法,与此正规正交矩阵的通用符号Q有关。 householderQR的分解方式的演示代码: 输出为: 矩阵的LU三角分解 矩阵的SVD(奇异值分解)分解 奇异值分解 (singular value decomposition,SVD) 是另一种正交矩阵分解法;SVD是最可靠的分解法,但是它比QR 分解法要花上近十倍的计算时间。[U,S,V]=svd(A),其中U和V分别代表两个正交矩阵,而S代表一对角矩阵。 和QR分解法相同, 原矩阵A不必为正方矩阵。使用SVD分解法的用途是解最小平方误差法和数据压缩。 矩阵的LLT分解 Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解。它要求矩阵的所有特征值必须大于零,故分解的下三角的对角元也是大于零的(LU三角分解法的变形)。 LDLT分解 对于对称正定矩阵A,可以将A分解为A=LDLT,并且该分解是唯一的,其中L为一下三角形单位矩阵(即主对角线元素皆为1),D为一对角矩阵(只在主对角线上有元素且不为0,其余皆为零),LT为L的转置矩阵。LDLT 是Cholesky分解的变形: LDLT分解法实际上是Cholesky分解法的改进,因为Cholesky分解法虽然不需要选主元,但其运算过程中涉及到开方问题,而LDLT分解法则避免了这一问题,可用于求解线性方程组。 = A.ldlt().solve(b)); // A sym. p.s.d. #include = A.llt().solve(b)); // A sym. p.d. #include = A.lu().solve(b)); // Stable and fast. #include = A.qr().solve(b)); // No pivoting. #include = A.svd().solve(b)); // Stable, slowest. #include // .ldlt() -> .matrixL() and .matrixD() 输出为: 结果不一致,这是为什么? 14、特征值分解 最小二乘法求解 最小二乘求解有两种方式,jacobiSvd或者colPivHouseholderQr,4*4以下的小矩阵速度没有区别,jacobiSvd可能更快,大矩阵最好用colPivHouseholderQr 输出为: 15、稀疏矩阵 稀疏矩阵的头文件包括: 初始化有两种方式: 1.使用三元组插入 2.直接将已知的非0值插入 稀疏矩阵支持大部分一元和二元运算: 二元运算中,稀疏矩阵和普通矩阵可以混合使用 16、C++数组和矩阵转换 使用Map函数,可以实现Eigen的矩阵和c++中的数组直接转换,语法如下: 矩阵的应用实例: 矩阵的定义 Eigen的基础使用 Eigen特殊矩阵生成 Eigen矩阵分块
Module
Header file
Contents
Core
Matrix and Array classes, basic linear algebra (including triangular and selfadjoint products), array manipulation
Geometry
Transform, Translation, Scaling, Rotation2D and 3D rotations (Quaternion, AngleAxis)
LU
Inverse, determinant, LU decompositions with solver (FullPivLU, PartialPivLU)
Cholesky
LLT and LDLT Cholesky factorization with solver
Householder
Householder transformations; this module is used by several linear algebra modules
SVD
SVD decompositions with least-squares solver (JacobiSVD, BDCSVD)
QR
QR decomposition with solver (HouseholderQR, ColPivHouseholderQR, FullPivHouseholderQR)
Eigenvalues
Eigenvalue, eigenvector decompositions (EigenSolver, SelfAdjointEigenSolver, ComplexEigenSolver)
Sparse
Sparse matrix storage and related basic linear algebra (SparseMatrix, SparseVector)
(see Quick reference guide for sparse matrices for details on sparse modules)
Includes Core, Geometry, LU, Cholesky, SVD, QR, and Eigenvalues header files
Includes Dense and Sparse header files (the whole Eigen library)
template
typedef matrixdouble, Dynamic, Dynamic> MarixXd;
typedef Matrixdouble , 3 , 1> Vector3d
Vector2d a(5.0, 6.0);
Vector3d b(5.0, 6.0, 7.0);
Vector4d c(5.0, 6.0, 7.0, 8.0);
#include
#include
#include
Here is the matrix m:
3 -1
2.5 1.5
Here is the vector v:
4
3
Matrix3f m;
m 1, 2, 3,
4, 5, 6,
7, 8, 9;
std::cout
1 2 3
4 5 6
7 8 9
VectorXf m_Vector_A;
MatrixXf m_matrix_B;
int m_iN =-1;
bool InitData(int pSrc[100][100], int iWidth, int iHeight)
{
if (NULL == pSrc || iWidth 0 || iHeight 0)
return false;
m_iN = iWidth*iHeight;
VectorXf tmp_A(m_iN);
MatrixXf tmp_B(m_iN, 5);
int i =0, j=0, iPos =0;
while(iiWidth)
{
j=0;
while(jiHeight)
{
tmp_A(iPos) = pSrc[i][j] * log((float)pSrc[i][j]);
tmp_B(iPos,0) = pSrc[i][j] ;
tmp_B(iPos,1) = pSrc[i][j] * i;
tmp_B(iPos,2) = pSrc[i][j] * j;
tmp_B(iPos,3) = pSrc[i][j] * i * i;
tmp_B(iPos,4) = pSrc[i][j] * j * j;
++iPos;
++j;
}
++i;
}
m_Vector_A = tmp_A;
m_matrix_B = tmp_B;
}
(3) 使用“=”操作符操作动态矩阵时,如果左右边的矩阵大小不等,则左边的动态矩阵的大小会被修改为右边的大小。例如下面的代码段:MatrixXf a(2,2);
std::cout "a is of size " "x" std::endl;
MatrixXf b(3,3);
a = b;
std::cout "a is now of size " "x"
a is of size 2x2
a is now of size 3x3
#include
a + b =
3 5
4 8
a - b =
-1 -1
2 0
Doing a += b;
Now a =
3 5
4 8
-v + w - v =
-1
-4
-6
#include
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
MatrixXcf a = MatrixXcf::Random(2,2);
cout "Here is the matrix a\n" endl;
cout "Here is the matrix a^T\n" endl;
cout "Here is the conjugate of a\n" endl;
cout "Here is the matrix a^*\n" endl;
cout "determinant: " endl;
cout "inverse:"
Here is the matrix a
(-0.211,0.68) (-0.605,0.823)
(0.597,0.566) (0.536,-0.33)
Here is the matrix a^T
(-0.211,0.68) (0.597,0.566)
(-0.605,0.823) (0.536,-0.33)
Here is the conjugate of a
(-0.211,-0.68) (-0.605,-0.823)
(0.597,-0.566) (0.536,0.33)
Here is the matrix a^*
(-0.211,-0.68) (0.597,-0.566)
(-0.605,-0.823) (0.536,0.33)
......#include
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 -0
2 0
Let‘s multiply mat by itself
Now mat is mat:
7 10
15 22
#include
matrix.block(i,j,p,q); (1)
matrix.block
#include
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
注意:
(1)需与获取矩阵的行数/列数的函数( rows(), cols() )的进行区别,不要弄混淆。
(2)函数参数为响应行/列的索引,需注意矩阵的行列均以0开始。
下面的代码段用于演示获取矩阵的指定行列:#include
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
获取向量的前n个元素:vector.head(n);
获取向量尾部的n个元素:vector.tail(n);
获取从向量的第i个元素开始的n个元素:vector.segment(i,n);
其用法可参考如下代码段:#include
v.head(3) =
1
2
3
v.tail3>() =
4
5
6
after ‘v.segment(1,4) *= 2‘, v =
1
4
6
8
10
6
#include
Decomposition
Method
Requirements on the matrix
Speed
Accuracy
PartialPivLU
partialPivLu()
Invertible(可逆)
++
+
FullPivLU
fullPivLu()
None
-
+++
HouseholderQR
householderQr()
None
++
+
ColPivHouseholderQR
colPivHouseholderQr()
None
+
++
FullPivHouseholderQR
fullPivHouseholderQr()
None
-
+++
LLT
llt()
Positive definite(正定)
+++
+
LDLT
ldlt()
Positive or negative semidefinite(正或负半定)
+++
++
void QR2()
{
Matrix3d A;
A1,1,1,
2,-1,-1,
2,-4,5;
HouseholderQRQR2(): HouseholderQR---------------------------------------------
A
1 1 1
2 -1 -1
2 -4 5
qr.matrixQR()
-3 3 -3
0.5 -3 3
0.5 -1 -3
R
-3 3 -3
0 -3 3
0 0 -3
Q
-0.333333 -0.666667 -0.666667
-0.666667 -0.333333 0.666667
-0.666667 0.666667 -0.333333
Q*R
1 1 1
2 -1 -1
2 -4 5
#include
Vector4d x2
Vector4d x3
Vector4d x4
Vector4d x5
Vector4d x6
// .llt() -> .matrixL()
// .lu() -> .matrixL() and .matrixU()
// .qr() -> .matrixQ() and .matrixR()
// .svd() -> .matrixU(), .singularValues(), and .matrixV() std::cout "The solution is:\n" "\n\n""\n\n"std::endl;
}
0
-1
-4
-3
-289.143
448.714
29.9082
3.97959
1.52903
0.1758
-0.340206
0.0423223
#include
Here is the matrix A:
0.680375 0.59688
-0.211234 0.823295
0.566198 -0.604897
Here is the right hand side b:
-0.329554
0.536459
-0.444451
The least-squares solution is:
-0.669626
0.314253
The least-squares solution is:
-0.669626
0.314253
#include
typedef Eigen::Tripletdouble> T;
std::vector
SparseMatrixdouble> mat(rows,cols);
mat.reserve(VectorXi::Constant(cols,6));
for(...)
{
// i,j 个非零值 v_ij != 0
mat.insert(i,j) = v_ij;
}
mat.makeCompressed(); // optional
sm1.real() sm1.imag() -sm1 0.5*sm1
sm1+sm2 sm1-sm2 sm1.cwiseProduct(sm2)
//dm表示普通矩阵
dm2 = sm1 + dm1;
//@param MatrixType 矩阵类型
//@param MapOptions 可选参数,指的是指针是否对齐,Aligned, or Unaligned. The default is Unaligned.
//@param StrideType 可选参数,步长
/*
Map
#include
// Basic usage
// Eigen // Matlab // comments
x.size() // length(x) // vector size
C.rows() // size(C,1) // number of rows
C.cols() // size(C,2) // number of columns
x(i) // x(i+1) // Matlab is 1-based
C(i, j) // C(i+1,j+1) //
A.resize(4, 4); // Runtime error if assertions are on.
B.resize(4, 9); // Runtime error if assertions are on.
A.resize(3, 3); // Ok; size didn‘t change.
B.resize(3, 9); // Ok; only dynamic cols changed.
A 1, 2, 3, // Initialize A. The elements can also be
4, 5, 6, // matrices, which are stacked along cols
7, 8, 9; // and then the rows are stacked.
B // B is three horizontally stacked A‘s.
A.fill(10); // Fill A with all 10‘s.
// Eigen // Matlab
MatrixXd::Identity(rows,cols) // eye(rows,cols)
C.setIdentity(rows,cols) // C = eye(rows,cols)
MatrixXd::Zero(rows,cols) // zeros(rows,cols)
C.setZero(rows,cols) // C = ones(rows,cols)
MatrixXd::Ones(rows,cols) // ones(rows,cols)
C.setOnes(rows,cols) // C = ones(rows,cols)
MatrixXd::Random(rows,cols) // rand(rows,cols)*2-1 // MatrixXd::Random returns uniform random numbers in (-1, 1).
C.setRandom(rows,cols) // C = rand(rows,cols)*2-1
VectorXd::LinSpaced(size,low,high) // linspace(low,high,size)‘
v.setLinSpaced(size,low,high) // v = linspace(low,high,size)‘
// Matrix slicing and blocks. All expressions listed here are read/write.
// Templated size versions are faster. Note that Matlab is 1-based (a size N // vector is x(1)...x(N)).
// Eigen // Matlab
x.head(n) // x(1:n)&