原文地址:http://www.cnblogs.com/goingupeveryday/p/5699053.html
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(); |
dotproduct 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-wiseor 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 columnor a rowof 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 n coeffs |
vec1.tail(n) |
vec1.tail<n>() |
the last n coeffs |
vec1.segment(pos,n) |
vec1.segment<n>(pos) |
the n coeffs in the range [ pos : pos + n - 1] |
Read-write access to sub-matrices: |
||
mat1.block(i,j,rows,cols) |
mat1.block<rows,cols>(i,j) |
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<rows,cols>() mat1.topRightCorner<rows,cols>() mat1.bottomLeftCorner<rows,cols>() mat1.bottomRightCorner<rows,cols>() |
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<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 diagonaland super/sub diagonalsof 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>()
|
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/