c++矩阵运算库Eigen

Eigen 的简介和下载安装

最近因为要写机械臂的运动学控制程序,因此难免就要在C++中进行矩阵运算。然而C++本身不支持矩阵运算,Qt库中的矩阵运算能力也相当有限,因此考虑寻求矩阵运算库Eigen的帮助。

Eigen是一个C++线性运算的模板库:他可以用来完成矩阵,向量,数值解等相关的运算。(Eigen is a C++ template library for linear algebra: matrices, vectors, numerical solvers, and related algorithms.)

Eigen库下载: http://eigen.tuxfamily.org/index.php?title=Main_Page

Eigen库的使用相当方便,压缩包中的Eigen文件夹拷贝到项目目录下,直接包含其中的头文件即可使用,省去了使用Cmake进行编译的烦恼。

Eigen官网有快速入门的参考文档:http://eigen.tuxfamily.org/dox/group__QuickRefPage.html

Eigen简单上手使用

要实现相应的功能只需要包含头相应的头文件即可:

[Core](http://eigen.tuxfamily.org/dox/group__Core__Module.html) #include [Matrix](http://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html)and[Array](http://eigen.tuxfamily.org/dox/classEigen_1_1Array.html)classes, basic linear algebra (including triangular and selfadjoint products), array manipulation
[Geometry](http://eigen.tuxfamily.org/dox/group__Geometry__Module.html) #include [Transform](http://eigen.tuxfamily.org/dox/classEigen_1_1Transform.html),[Translation](http://eigen.tuxfamily.org/dox/classEigen_1_1Translation.html),[Scaling](http://eigen.tuxfamily.org/dox/classScaling.html),[Rotation2D](http://eigen.tuxfamily.org/dox/classEigen_1_1Rotation2D.html)and 3D rotations ([Quaternion](http://eigen.tuxfamily.org/dox/classEigen_1_1Quaternion.html),[AngleAxis](http://eigen.tuxfamily.org/dox/classEigen_1_1AngleAxis.html))
[LU](http://eigen.tuxfamily.org/dox/group__LU__Module.html) #include Inverse, determinant, LU decompositions with solver ([FullPivLU](http://eigen.tuxfamily.org/dox/classEigen_1_1FullPivLU.html),[PartialPivLU](http://eigen.tuxfamily.org/dox/classEigen_1_1PartialPivLU.html))
[Cholesky](http://eigen.tuxfamily.org/dox/group__Cholesky__Module.html) #include [LLT](http://eigen.tuxfamily.org/dox/classEigen_1_1LLT.html)and[LDLT](http://eigen.tuxfamily.org/dox/classEigen_1_1LDLT.html)Cholesky factorization with solver
[Householder](http://eigen.tuxfamily.org/dox/group__Householder__Module.html) #include Householder transformations; this module is used by several linear algebra modules
[SVD](http://eigen.tuxfamily.org/dox/group__SVD__Module.html) #include SVD decomposition with least-squares solver ([JacobiSVD](http://eigen.tuxfamily.org/dox/classEigen_1_1JacobiSVD.html))
[QR](http://eigen.tuxfamily.org/dox/group__QR__Module.html) #include QR decomposition with solver ([HouseholderQR](http://eigen.tuxfamily.org/dox/classEigen_1_1HouseholderQR.html),[ColPivHouseholderQR](http://eigen.tuxfamily.org/dox/classEigen_1_1ColPivHouseholderQR.html),[FullPivHouseholderQR](http://eigen.tuxfamily.org/dox/classEigen_1_1FullPivHouseholderQR.html))
[Eigenvalues](http://eigen.tuxfamily.org/dox/group__Eigenvalues__Module.html) #include Eigenvalue, eigenvector decompositions ([EigenSolver](http://eigen.tuxfamily.org/dox/classEigen_1_1EigenSolver.html),[SelfAdjointEigenSolver](http://eigen.tuxfamily.org/dox/classEigen_1_1SelfAdjointEigenSolver.html),[ComplexEigenSolver](http://eigen.tuxfamily.org/dox/classEigen_1_1ComplexEigenSolver.html))
[Sparse](http://eigen.tuxfamily.org/dox/group__QuickRefPage.html) #include Sparse matrix storage and related basic linear algebra ([SparseMatrix](http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html), DynamicSparseMatrix,[SparseVector](http://eigen.tuxfamily.org/dox/classEigen_1_1SparseVector.html))
#include Includes Core, Geometry, LU, Cholesky, SVD, QR, and Eigenvalues header files
#include Includes Dense and Sparse header files (the whole[Eigen](http://eigen.tuxfamily.org/dox/namespaceEigen.html)library)

基本的矩阵运算只需要包含Dense即可

基本的数据类型:Array, matrix and vector

声明:

#include<Eigen/Dense>

...
//1D objects
Vector4d  v4;
Vector2f  v1(x, y);
Array3i   v2(x, y, z);
Vector4d  v3(x, y, z, w);
VectorXf  v5; // empty object
ArrayXf   v6(size);

//2D objects
atrix4f  m1;
MatrixXf  m5; // empty object
MatrixXf  m6(nb_rows, nb_columns);

赋值:

//    
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;
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

Resize矩阵:

//
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);

元素读取:

vector(i);    
vector[i];

matrix(i,j);

矩阵的预定义:

//vector x
x.setZero();
x.setOnes();
x.setConstant(value);
x.setRandom();
x.setLinSpaced(size, low, high);

VectorXf::Unit(4,1) == Vector4f(0,1,0,0)
                    == Vector4f::UnitY()

//matrix x
x.setZero(rows, cols);
x.setOnes(rows, cols);
x.setConstant(rows, cols, value);
x.setRandom(rows, cols);

x.setIdentity();//单位矩阵
x.setIdentity(rows, cols);

映射操作,可以将c语言类型的数组映射为矩阵或向量:

(注意:

1.映射只适用于对一维数组进行操作,如果希望将二维数组映射为对应的矩阵,可以借助"mat.row(i)=Map v(data[i],n)"进行循环来实现,其中data为二维数组名,n为数组第一维的维度。

2.应设置之后数组名和矩阵或向量名其实指向同一个地址,只是名字和用法不同,另外,对其中的任何一个进行修改都会修改另外一个)

//连续映射
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

//间隔映射
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|

算术运算:

add

subtract
mat3 = mat1 + mat2; mat3 += mat1;mat3 = mat1 - mat2; mat3 -= mat1;
scalar product mat3 = mat1 * s1; mat3 *= s1; mat3 = s1 * mat1;mat3 = mat1 / s1; mat3 /= s1;
matrix/vector

products[*](http://eigen.tuxfamily.org/dox/group__QuickRefPage.html#matrixonly)
col2 = mat1 * col1;row2 = row1 * mat1; row1 *= mat1;mat3 = mat1 * mat2; mat3 *= mat1;
transposition

adjoint[*](http://eigen.tuxfamily.org/dox/group__QuickRefPage.html#matrixonly)
mat1 = mat2.transpose(); mat1.transposeInPlace();mat1 = mat2.adjoint(); mat1.adjointInPlace();
[dot](http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html#adb71ddef4955ae7d353df12d05665191)product

inner product[*](http://eigen.tuxfamily.org/dox/group__QuickRefPage.html#matrixonly)
scalar = vec1.dot(vec2);scalar = col1.adjoint() * col2;scalar = (col1.adjoint() * col2).value();
outer product[*](http://eigen.tuxfamily.org/dox/group__QuickRefPage.html#matrixonly) mat = col1 * col2.transpose();

[norm](http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html#a0be1b433c65ce9d92c81a4718daf54e5)

[normalization](http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html#a8ed1fb2e792b1079639a74e3581fbc74)[*](http://eigen.tuxfamily.org/dox/group__QuickRefPage.html#matrixonly)
scalar = vec1.norm(); scalar = vec1.squaredNorm()vec2 = vec1.normalized(); vec1.normalize();// inplace

[cross product](http://eigen.tuxfamily.org/dox/group__QuickRefPage.html)[*](http://eigen.tuxfamily.org/dox/group__QuickRefPage.html#matrixonly) #include vec3 = vec1.cross(vec2);

Coefficient-wise & Array operators

Coefficient-wise operators for matrices and vectors:

[Matrix](http://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html)API[*](http://eigen.tuxfamily.org/dox/group__QuickRefPage.html#matrixonly) Via[Array](http://eigen.tuxfamily.org/dox/classEigen_1_1Array.html)conversions
mat1.cwiseMin(mat2)mat1.cwiseMax(mat2)mat1.cwiseAbs2()mat1.cwiseAbs()mat1.cwiseSqrt()mat1.cwiseProduct(mat2)mat1.cwiseQuotient(mat2) mat1.array().min(mat2.array())mat1.array().max(mat2.array())mat1.array().abs2()mat1.array().abs()mat1.array().sqrt()mat1.array() * mat2.array()mat1.array() / mat2.array()

Arrayoperators:*

Arithmetic operators array1 * array2 array1 / array2 array1 *= array2 array1 /= array2array1 + scalar array1 - scalar array1 += scalar array1 -= scalar
Comparisons array1 < array2 array1 > array2 array1 < scalar array1 > scalararray1 <= array2 array1 >= array2 array1 <= scalar array1 >= scalararray1 == array2 array1 != array2 array1 == scalar array1 != scalar
Trigo, power, and

misc functions

and the STL variants
array1.min(array2)array1.max(array2)array1.abs2()array1.abs() abs(array1)array1.sqrt() sqrt(array1)array1.log() log(array1)array1.exp() exp(array1)array1.pow(exponent) pow(array1,exponent)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)

Reductions

Eigenprovides several reduction methods such as:minCoeff(),maxCoeff(),sum(),prod(),trace()*,norm()*,squaredNorm()*,all(), andany(). All reduction operations can be done matrix-wise,column-wiseorrow-wise. Usage example:

5 3 1mat = 2 7 89 4 6 mat.minCoeff(); 1
mat.colwise().minCoeff(); 2 3 1
mat.rowwise().minCoeff(); 124

Typical use cases of all() and any():

//Typical use cases of all() and any():

if((array1 > 0).all()) ...      // if all coefficients of array1 are greater than 0 ...
if((array1 < array2).any()) ... // if there exist a pair i,j such that array1(i,j) < array2(i,j) ...

Sub-matrices

Read-write access to acolumnor arowof a matrix (or array):
mat1.row(i) = mat2.col(j);mat1.col(j1).swap(mat1.col(j2));
Read-write access to sub-vectors:

Default versions Optimized versions when the size

is known at compile time
vec1.head(n) vec1.head() the first`n`coeffs
vec1.tail(n) vec1.tail() the last`n`coeffs
vec1.segment(pos,n) vec1.segment(pos) the`n`coeffs in the

range [`pos`:`pos`+`n`- 1]

Read-write access to sub-matrices:

mat1.block(i,j,rows,cols)[(more)](http://eigen.tuxfamily.org/dox/classEigen_1_1DenseBase.html#a1dbaa2fc7b809720407130f48dfacf8f)mat1.block(i,j)[(more)](http://eigen.tuxfamily.org/dox/classEigen_1_1DenseBase.html#a3e433315822db2811a65e88c70672743)the`rows`x`cols`sub-matrix

starting from position (`i`,`j`)
mat1.topLeftCorner(rows,cols)mat1.topRightCorner(rows,cols)mat1.bottomLeftCorner(rows,cols)mat1.bottomRightCorner(rows,cols)mat1.topLeftCorner()mat1.topRightCorner()mat1.bottomLeftCorner()mat1.bottomRightCorner()the`rows`x`cols`sub-matrix

taken in one of the four corners
mat1.topRows(rows)mat1.bottomRows(rows)mat1.leftCols(cols)mat1.rightCols(cols)mat1.topRows()mat1.bottomRows()mat1.leftCols()mat1.rightCols()specialized versions of block()

when the block fit two corners

top

Miscellaneous operations

Reverse

Vectors, rows, and/or columns of a matrix can be reversed (seeDenseBase::reverse(),DenseBase::reverseInPlace(),VectorwiseOp::reverse()).
vec.reverse() mat.colwise().reverse() mat.rowwise().reverse()vec.reverseInPlace()

Replicate

Vectors, matrices, rows, and/or columns can be replicated in any direction (seeDenseBase::replicate(),VectorwiseOp::replicate())
vec.replicate(times) vec.replicatemat.replicate(vertical_times, horizontal_times) mat.replicate()mat.colwise().replicate(vertical_times, horizontal_times) mat.colwise().replicate()mat.rowwise().replicate(vertical_times, horizontal_times) mat.rowwise().replicate()
top

Diagonal, Triangular, and Self-adjoint matrices

(matrix world*)

Diagonal matrices

Operation Code
view a vector[as a diagonal matrix](http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html#adaf22d3a2069ec2c0df912cb87329e9c)

mat1 = vec1.asDiagonal();
Declare a diagonal matrix DiagonalMatrix diag1(size);diag1.diagonal() = vector;
Access the[diagonal](http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html#a0a45dd0ed5a44ec3f8f43239f2e4ac25)and[super/sub diagonals](http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html#ac062b66ec8aa4eb1f12cbf4ea3ff4f17)of a matrix as a vector (read/write) vec1 = mat1.diagonal(); mat1.diagonal() = vec1;// main diagonalvec1 = mat1.diagonal(+n); mat1.diagonal(+n) = vec1;// n-th super diagonalvec1 = mat1.diagonal(-n); mat1.diagonal(-n) = vec1;// n-th sub diagonalvec1 = mat1.diagonal<1>(); mat1.diagonal<1>() = vec1;// first super diagonalvec1 = mat1.diagonal<-2>(); mat1.diagonal<-2>() = vec1;// second sub diagonal

Optimized products and inverse mat3 = scalar * diag1 * mat1;mat3 += scalar * mat1 * vec1.asDiagonal();mat3 = vec1.asDiagonal().inverse() * mat1mat3 = mat1 * diag1.inverse()

Triangular views

TriangularViewgives a view on a triangular part of a dense matrix and allows to perform optimized operations on it. The opposite triangular part is never referenced and can be used to store other information.

Note
The .triangularView() template member function requires the`template`keyword if it is used on an object of a type that depends on a template parameter; see[The template and typename keywords in C++](http://eigen.tuxfamily.org/dox/TopicTemplateKeyword.html)for details.
Operation Code
Reference to a triangular with optional

unit or null diagonal (read/write):
m.triangularView()

`Xxx`=[Upper](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564ae70afef0d3ff7aca74e17e85ff6c9f2e),[Lower](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564af886b397626076218462d53d50eb96bc),[StrictlyUpper](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564abf11791f004a059cfdd9b941c76f3703),[StrictlyLower](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564a29eb98bd08096415c55f37ed4ac2af11),[UnitUpper](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564a65b1d67b2bb2e4a85b5f6a8863cd7109),[UnitLower](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564a0dc6c411b3fc7ae6e32860a7872b7d18)
Writing to a specific triangular part:

(only the referenced triangular part is evaluated)
m1.triangularView<[Eigen::Lower](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564af886b397626076218462d53d50eb96bc)>() = m2 + m3
Conversion to a dense matrix setting the opposite triangular part to zero: m2 = m1.triangularView<[Eigen::UnitUpper](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564a65b1d67b2bb2e4a85b5f6a8863cd7109)>()
Products: m3 += s1 * m1.adjoint().triangularView<[Eigen::UnitUpper](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564a65b1d67b2bb2e4a85b5f6a8863cd7109)>() * m2m3 -= s1 * m2.conjugate() * m1.adjoint().triangularView<[Eigen::Lower](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564af886b397626076218462d53d50eb96bc)>()
Solving linear equations:

![$ M_2 := L_1^{-1} M_2 $](http://eigen.tuxfamily.org/dox/form_170.png)

![$ M_3 := {L_1^*}^{-1} M_3 $](http://eigen.tuxfamily.org/dox/form_171.png)

![$ M_4 := M_4 U_1^{-1} $](http://eigen.tuxfamily.org/dox/form_172.png)

L1.triangularView<[Eigen::UnitLower](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564a0dc6c411b3fc7ae6e32860a7872b7d18)>().solveInPlace(M2)L1.triangularView().adjoint().solveInPlace(M3)U1.triangularView().solveInPlace<[OnTheRight](http://eigen.tuxfamily.org/dox/group__enums.html#gga1ee4b3b0c0d87990a374f7ef5b551fcfaeda0d7b1859ec757de18ee3b7c6c541c)>(M4)

Symmetric/selfadjoint views

Just as for triangular matrix, you can reference any triangular part of a square matrix to see it as a selfadjoint matrix and perform special and optimized operations. Again the opposite triangular part is never referenced and can be used to store other information.

Note
The .selfadjointView() template member function requires the`template`keyword if it is used on an object of a type that depends on a template parameter; see[The template and typename keywords in C++](http://eigen.tuxfamily.org/dox/TopicTemplateKeyword.html)for details.
Operation Code
Conversion to a dense matrix: m2 = m.selfadjointView<[Eigen::Lower](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564af886b397626076218462d53d50eb96bc)>();
Product with another general matrix or vector: m3 = s1 * m1.conjugate().selfadjointView<[Eigen::Upper](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564ae70afef0d3ff7aca74e17e85ff6c9f2e)>() * m3;m3 -= s1 * m3.adjoint() * m1.selfadjointView<[Eigen::Lower](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564af886b397626076218462d53d50eb96bc)>();
Rank 1 and rank K update:

![$ upper(M_1) \mathrel{{+}{=}} s_1 M_2 M_2^* $](http://eigen.tuxfamily.org/dox/form_173.png)

![$ lower(M_1) \mathbin{{-}{=}} M_2^* M_2 $](http://eigen.tuxfamily.org/dox/form_174.png)

M1.selfadjointView<[Eigen::Upper](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564ae70afef0d3ff7aca74e17e85ff6c9f2e)>().rankUpdate(M2,s1);M1.selfadjointView<[Eigen::Lower](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564af886b397626076218462d53d50eb96bc)>().rankUpdate(M2.adjoint(),-1);
Rank 2 update: (![$ M \mathrel{{+}{=}} s u v^* + s v u^* $](http://eigen.tuxfamily.org/dox/form_175.png)) M.selfadjointView<[Eigen::Upper](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564ae70afef0d3ff7aca74e17e85ff6c9f2e)>().rankUpdate(u,v,s);
Solving linear equations:

(![$ M_2 := M_1^{-1} M_2 $](http://eigen.tuxfamily.org/dox/form_176.png))
// via a standard Cholesky factorizationm2 = m1.selfadjointView<[Eigen::Upper](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564ae70afef0d3ff7aca74e17e85ff6c9f2e)>().llt().solve(m2);// via a Cholesky factorization with pivotingm2 = m1.selfadjointView<[Eigen::Lower](http://eigen.tuxfamily.org/dox/group__enums.html#gga7a42ebf69230517a2a10b2b8aa785564af886b397626076218462d53d50eb96bc)>().ldlt().solve(m2);

更多内容请关注我的博客:http://markma.tk/
原文链接: https://www.cnblogs.com/goingupeveryday/p/5699053.html

欢迎关注

微信关注下方公众号,第一时间获取干货硬货;公众号内回复【pdf】免费获取数百本计算机经典书籍

原创文章受到原创版权保护。转载请注明出处:https://www.ccppcoding.com/archives/237670

非原创文章文中已经注明原地址,如有侵权,联系删除

关注公众号【高性能架构探索】,第一时间获取最新文章

转载文章受原作者版权保护。转载请注明原作者出处!

(0)
上一篇 2023年2月13日 下午5:29
下一篇 2023年2月13日 下午5:30

相关推荐