贡献者: addis
Eigen 的主页。
#include <Eigen/Dense>
等头文件即可。
apt install libeigen3-dev
Eigen
// 矩阵
// 第一个维度固定在 3, 第二个维度可以动态变化
Matrix<double, 3, Dynamic> A;
// 两个动态维度, Aligned 可以把起始地址 padding 到一个 2^n 的整数倍的地址
// 以满足 SIMD 优化需求, 默认开启。 RowMajor 是行主序矩阵,默认是 ColMajor
Matrix<double, Dynamic, Dynamic, Aligned | RowMajor> A
typedef Matrix<int, 4, 4> Matrix4i;
typedef Matrix<double, 4, 4> Matrix4d;
typedef Matrix<double, Dynamic, Dynamic> MatrixXd;
typedef Matrix<double, Dynamic, 1> VectorXd; // 矢量区分行和列
typedef Matrix<double, 1, Dynamic> RowVectorXd;
// declaration & initialization
MatrixXd a(10, 15); // No initialization
a.setRandom(m,n); // 随机赋值
a.setZero(); 把矩阵赋值为 0;
a.setZero(m,n); 把矩阵赋值为 m*n 的 0 矩阵;
以下函数用法类似
setOnes(); setConstant(); setIdentity();(单位矩阵) setRandom(); setLinSpaced();(等间距赋值)
注意其中 setConstant()
输入 1 个变量时与 fill()
功能相同,输入 3 个变量时最后一个为常数。setLinSpaced(size, val1, val2)
只能对 Vector 或者 RowVector 使用。如果想赋值给 MatrixXd
, 可以用 MatrixXd a = VectorXd::LinSpaced(3, 3, 4)
; 对于整型,setLinSpaced
不保证最后一个值等于 val2, 而是保证间隔相等。这时候只能先生成 Xd, 然后 .cast<int>.
* 逗号赋值都是 RowMaor 的,无论 a 是什么 major.
a << 1,2,3,4;
获取信息
Matrix<>::size()
相当于 Matlab 的 numel()
, 另外,rows()
和 cols()
分别是行数和列数。
Matrix<>::data()
可以获得 Matrix 数据的指针,用于直接读写矩阵数据。注意 rowwise 需要专门声明。也可以用 &
来获取矩阵元的地址。
Matrix<>::transpose()
复数矩阵共轭用 conjugate()
, 共轭转置用 adjoint()
Matrix<>::sum(); Matrix<>::colwise().sum(); Matrix<>::rowwise().sum()
v.dot(w), v.cross(w)
类似的有 mean(), prod(), all(), any(), maxCoeff(), minCoeff()
获取子矩阵
对 Vector, 有 Matrix<>::head(n), tail(n), segment(n)
. 对 Matrix, 有 block(i,j,rows,cols)
, 有 row(i), col(j), topRows(n), middleRows(i, rows), bottomRows(n), leftCols(n), middleCols(j, cols), rightCols(n), topLeftCorner(rows,cols), bottomRightCorner(rows,cols)
等。
Matrix<>::cwiseAbs()
cwiseInverse(); cwiseMax(); cwiseMin(); cwiseSign(); cwiseSqrt(); Matrix<>::array()
用于把 Matrix 转换为 array, 逐个元素运算,或加减一个常数。如 Matrix<>::array().square()
round(); pow(); sqrt();
注意 array 可以直接赋值给 matrix.
mat.colwise() = v;
类似地,+=, -=
等也可以使用。
如果矩阵同时出现在等号两边,就有可能出现 Aliasing. 例如 MatrixXi mat(3,3); mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);
这时,要在等号右边的加上 .eval()
. 特殊地,对于转置等,可以直接用 mat = mat.transposeInPlace()
. 类似地有 adjointInPlace()
.
要用 SVD, 要 #include<Eigen/SVD>
, 例程:
MatrixXd A(3,3), U, V, X;
A << 1, 2, 3, 4, 5, 6, 7, 8, 9;
BDCSVD<MatrixXd> svd(A, ComputeThinU | ComputeThinV);
U = svd.matrixU(); V = svd.matrixV(); X = svd.singularValues();
这里介绍如何直接操作 Matrix<>
底层的数据,以及如何将内存中已有的数据直接给 Eigen 使用而无需复制。
Matrix<>::data()
可以返回第一个矩阵元的指针。然后就可以当做矩阵来用了。
Matrix<>::innerStride()
一律返回 1, outerStride()
返回 innerSize()
, 而 innerSize
对于 c-major 是 rows()
, 对于 r-major 是 cols()
. 这就说明矩阵的内存显然是连续的(无论怎么优化), 可以放心使用。
Matrix<>
到底是 col-major 还是 row-major, 默认是 col-major. 不推荐用 #define EIGEN_DEFAULT_TO_ROW_MAJOR
将默认改为 row-major. 建议用例如 typedef Matrix<double, Dynamic, Dynamic, Aligned | RowMajor> RMatrixXd.
Map<>
类即可,所有接受 Matrix<>
的函数同样也接受 Map<>
(除了自己写的函数)
Map<>
有三个模板参数,但是只需要用一个,即 Map<Matrix<>>
.
Map<MatrixXf> a(duuble* pa, rows,cols);
Eigen 使用模板变成的巨大优势就是可以使用自定义标量类型。代码已经默认支持 float,double,long double
以及对应的复数。要让矩阵元支持其他类型(如 g++ 的 __float128
,甚至 GMP 的任意精度类型),只需要告诉 Eigen 新增标量的特性即可:
namespace Eigen {
template<> struct NumTraits<Qdoub> : GenericNumTraits<Qdoub>
{
typedef Qdoub Real;
typedef Qdoub NonInteger;
typedef Qdoub Nested;
static inline Real epsilon() { return FLT128_EPSILON; }
static inline Real dummy_precision() { return 0; }
static inline int digits10() { return FLT128_DIG; }
enum {
IsInteger = 0,
IsSigned = 1,
IsComplex = 0,
RequireInitialization = 0,
ReadCost = 2,
AddCost = 2,
MulCost = 6
};
};
template<> struct NumTraits<Qcomp> : GenericNumTraits<Qcomp>
{
typedef Qdoub Real;
typedef Qdoub NonInteger;
typedef Qdoub Nested;
static inline Real epsilon() { return FLT128_EPSILON; }
static inline Real dummy_precision() { return 0; }
static inline int digits10() { return FLT128_DIG; }
enum {
IsInteger = 0,
IsSigned = 0,
IsComplex = 1,
RequireInitialization = 0,
ReadCost = 4,
AddCost = 4,
MulCost = 12
};
};
typedef Matrix<Qdoub, Dynamic, Dynamic> MatrixXq;
typedef Matrix<Qcomp, Dynamic, Dynamic> MatrixXqc;
typedef Matrix<Qdoub, Dynamic, 1> VectorXq;
typedef Matrix<Qcomp, Dynamic, 1> VectorXqc;
} // namespace Eigen
参考:稀疏矩阵文档主页,SparseMatrix 类文档,解稀疏线性方程文档。
Eigen::SparseMatrix<Scalar_, Options_, StorageIndex_>
模板中,Scalar_
是矩阵元的类型,Options_
只能是 Eigen::ColMajor
或 Eigen::RowMajor
,默认是 ColMajor
。StorageIndex_
是内部数组的 index 的类型,必须是有符号整数,默认 int
。注意这不是成员函数中行标和列表的 index 类型(Eigen::Index
)。
除非特殊说明,我们下面都以 RowMajor
为例。
RowMajor
的稀疏矩阵类似 CSR(子节 3 ),但允许储存矩阵元的数组 Values
中有一些预留的空位。
RowMajor
的 inner dimension 是行的方向,outer dimmension 是列。ColMajor 则相反。
InnerIndices
是非零矩阵元的列标,大小元素都和 Values
一一对应,空位处的值无意义。
OuterStarts
:Values[OuterStarts[i]]
是第 i
行的第一个非零元。如果该行为零,那就和下一个非零元的 OuterStarts
相同。OuterStarts
的尺寸比列数多 1,最后一个元素是 Values
的长度。
InnerNNZs
是每一列非零元素的个数(不包含空位),长度和列数相同。
Scalar coeff(Index row, Index col) const
Scalar &coeffRef(Index row, Index col)
Index cols() const
列数
Index rows() const
行数
void conservativeResize(Index rows, Index cols)
改变矩阵尺寸,保留原来的值。
StorageIndex *innerIndexPtr()
返回 InnerIndices
数组的指针
StorageIndex *innerNonZeroPtr()
返回 InnerNNZs
数组的指针
Scalar &insert(Index row, Index col)
返回一个引用,用于插入一个值。该矩阵元之前不能存在值
bool isCompressed() const
是否是 CSR 格式
void makeCompressed()
变为 CSR 格式
sqeeze()
类似 makeCompressed()
,但仍然会保留少量空位。
Index nonZeros() const
非零元的总个数
StorageIndex *outerIndexPtr()
返回 OuterStarts
数组的指针
void setFromTriplets(const InputIterators &begin, const InputIterators &end)
从 vector<Eigen::Triplet<Scalar>>
给稀疏矩阵赋值(不一定是 vector
)。
Scalar *valuePtr()
返回 Values
数组的指针
SLISC/tests/eigen.cpp
ConjugateGradient<SparseMatrix<Comp>, Lower|Upper> solver(H);
solver.make(H)
(A*x - b).norm()/b.norm()
。