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 | #include <Eigen/Core> | Matrix and Array classes, basic linear algebra (including triangular and selfadjoint products), array manipulation |
Geometry | #include <Eigen/Geometry> | Transform, Translation, Scaling, Rotation2D and 3D rotations (Quaternion, AngleAxis) |
LU | #include <Eigen/LU> | Inverse, determinant, LU decompositions with solver (FullPivLU, PartialPivLU) |
Cholesky | #include <Eigen/Cholesky> | LLT and LDLT Cholesky factorization with solver |
Householder | #include <Eigen/Householder> | Householder transformations; this module is used by several linear algebra modules |
SVD | #include <Eigen/SVD> | SVD decomposition with least-squares solver (JacobiSVD) |
QR | #include <Eigen/QR> | QR decomposition with solver (HouseholderQR, ColPivHouseholderQR, FullPivHouseholderQR) |
Eigenvalues | #include <Eigen/Eigenvalues> | Eigenvalue, eigenvector decompositions (EigenSolver, SelfAdjointEigenSolver,ComplexEigenSolver) |
Sparse | #include <Eigen/Sparse> | Sparse matrix storage and related basic linear algebra (SparseMatrix, DynamicSparseMatrix,SparseVector) |
#include <Eigen/Dense> | Includes Core, Geometry, LU, Cholesky, SVD, QR, and Eigenvalues header files | |
#include <Eigen/Eigen> | Includes Dense and Sparse header files (the whole Eigen 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<Vector> 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 * |
col2 = mat1 * col1; row2 = row1 * mat1; row1 *= mat1; mat3 = mat1 * mat2; mat3 *= mat1; |
transposition adjoint * |
mat1 = mat2.transpose(); mat1.transposeInPlace(); mat1 = mat2.adjoint(); mat1.adjointInPlace(); |
dot product inner product * |
scalar = vec1.dot(vec2); scalar = col1.adjoint() * col2; scalar = (col1.adjoint() * col2).value(); |
outer product * |
mat = col1 * col2.transpose();
|
norm normalization * |
scalar = vec1.norm(); scalar = vec1.squaredNorm()
vec2 = vec1.normalized(); vec1.normalize(); // inplace
|
cross product * | #include <Eigen/Geometry> vec3 = vec1.cross(vec2); |
Coefficient-wise & Array operators
Coefficient-wise operators for matrices and vectors:
Matrix API * | Via Array 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() |
Arithmetic operators | array1 * array2 array1 / array2 array1 *= array2 array1 /= array2 array1 + scalar array1 - scalar array1 += scalar array1 -= scalar |
Comparisons | array1 < array2 array1 > array2 array1 < scalar array1 > scalar array1 <= array2 array1 >= array2 array1 <= scalar array1 >= scalar array1 == 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
Eigen provides several reduction methods such as: minCoeff() , maxCoeff() , sum() , prod() , trace() *, norm()*, squaredNorm() *, all() , and any() . All reduction operations can be done matrix-wise, column-wise or row-wise. Usage example:
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 |
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 a column or a row of 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<n>() | the first ncoeffs |
vec1.tail(n) | vec1.tail<n>() | the last ncoeffs |
vec1.segment(pos,n) | vec1.segment<n>(pos) | the ncoeffs in the range [ pos: pos+ n- 1] |
Read-write access to sub-matrices: |
||
mat1.block(i,j,rows,cols) (more) | mat1.block<rows,cols>(i,j) (more) | the rowsx colssub-matrix starting from position ( i, j) |
mat1.topLeftCorner(rows,cols) mat1.topRightCorner(rows,cols) mat1.bottomLeftCorner(rows,cols) mat1.bottomRightCorner(rows,cols) | mat1.topLeftCorner<rows,cols>() mat1.topRightCorner<rows,cols>() mat1.bottomLeftCorner<rows,cols>() mat1.bottomRightCorner<rows,cols>() | the rowsx colssub-matrix taken in one of the four corners |
mat1.topRows(rows) mat1.bottomRows(rows) mat1.leftCols(cols) mat1.rightCols(cols) | mat1.topRows<rows>() mat1.bottomRows<rows>() mat1.leftCols<cols>() mat1.rightCols<cols>() | specialized versions of block() when the block fit two corners |
Miscellaneous operations
Reverse
Vectors, rows, and/or columns of a matrix can be reversed (see DenseBase::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 (see DenseBase::replicate(),VectorwiseOp::replicate())
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>()Diagonal, Triangular, and Self-adjoint matrices
(matrix world *)
Diagonal matrices
Operation | Code |
---|---|
view a vector as a diagonal matrix |
mat1 = vec1.asDiagonal(); |
Declare a diagonal matrix | DiagonalMatrix<Scalar,SizeAtCompileTime> diag1(size); diag1.diagonal() = vector; |
Access the diagonal and super/sub diagonals of a matrix as a vector (read/write) |
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
|
Optimized products and inverse |
mat3 = scalar * diag1 * mat1;
mat3 += scalar * mat1 * vec1.asDiagonal();
mat3 = vec1.asDiagonal().inverse() * mat1
mat3 = mat1 * diag1.inverse()
|
Triangular views
TriangularView gives 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++ for details.
Operation | Code |
---|---|
Reference to a triangular with optional unit or null diagonal (read/write): |
m.triangularView<Xxx>()
Xxx= Upper, Lower, StrictlyUpper, StrictlyLower, UnitUpper,UnitLower |
Writing to a specific triangular part: (only the referenced triangular part is evaluated) |
m1.triangularView<Eigen::Lower>() = m2 + m3 |
Conversion to a dense matrix setting the opposite triangular part to zero: | m2 = m1.triangularView<Eigen::UnitUpper>() |
Products: | m3 += s1 * m1.adjoint().triangularView<Eigen::UnitUpper>() * m2 m3 -= s1 * m2.conjugate() * m1.adjoint().triangularView<Eigen::Lower>() |
Solving linear equations:![]() ![]() ![]() |
L1.triangularView<Eigen::UnitLower>().solveInPlace(M2) L1.triangularView<Eigen::Lower>().adjoint().solveInPlace(M3) U1.triangularView<Eigen::Upper>().solveInPlace<OnTheRight>(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++ for details.
Operation | Code |
---|---|
Conversion to a dense matrix: | m2 = m.selfadjointView<Eigen::Lower>(); |
Product with another general matrix or vector: | m3 = s1 * m1.conjugate().selfadjointView<Eigen::Upper>() * m3; m3 -= s1 * m3.adjoint() * m1.selfadjointView<Eigen::Lower>(); |
Rank 1 and rank K update: ![]() ![]() |
M1.selfadjointView<Eigen::Upper>().rankUpdate(M2,s1); M1.selfadjointView<Eigen::Lower>().rankUpdate(M2.adjoint(),-1); |
Rank 2 update: (
![]() |
M.selfadjointView<Eigen::Upper>().rankUpdate(u,v,s); |
Solving linear equations: ( ![]() |
// 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); |
更多内容请关注我的博客:http://markma.tk/
- C++中四种类型转换
- C++多态性
- 设计模式之适配器模式(Adapter)C++实现
- 二维vector初始化
- 【C语言】关于斐波那契数的求法
- c++map的用法
- C语言范例学习03-上
- C++ string使用for改变为何必须使用引用
- 数据的压缩存储与解压缩算法实现(C语言)
- 数据的压缩存储与解压缩算法实现(C语言)
- 数据的压缩存储与解压缩算法实现(C语言)
- C++虚函数表分析
- C/C++中的预编译指令
- C++指针详解(复习一下)
- codeforces 495B - Modular Equations
- 编译预处理
- codeforces - 495A Digital Counter
- java socket解析和发送二进制报文工具(附java和C++转化问题)
- c++ primer学习(三)表达式
- SIFT算法C语言逐步实现详解(上)