Eigen （C++ 线性代数库）笔记

贡献者： addis

1. 安装

　　 Eigen 的主页。

1. 首先下载 Eigen, 不需要任何安装。
2. 把解压后的文件夹加入 include 路径。
3. 在使用 Eigen 的代码中加入 #include <Eigen/Dense> 等头文件即可。

2. 基础

• 与 Matlab 的语法对照表。
• ascii 快速参考
• 所有的 Eigen 名字都有 namespace 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;

3. operations

　　 获取信息

• 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.
• 复制 vector 给矩阵的每一列 mat.colwise() = v; 类似地，+=, -= 等也可以使用。

Aliasing

　　 如果矩阵同时出现在等号两边，就有可能出现 Aliasing. 例如 MatrixXi mat(3,3); mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2); 这时，要在等号右边的加上 .eval(). 特殊地，对于转置等，可以直接用 mat = mat.transposeInPlace(). 类似地有 adjointInPlace().

SVD (支持 complex)

　　 要用 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();

4. 外部接口

　　 这里介绍如何直接操作 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 object 如 Map<MatrixXf> a(duuble* pa, rows,cols);

5. 其他标量类型

　　 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

6. 稀疏矩阵（Sparse Matrix）

　　 参考：稀疏矩阵文档主页，SparseMatrix 类文档，解稀疏线性方程文档。

　　 Eigen::SparseMatrix<Scalar_, Options_, StorageIndex_> 模板中，Scalar_ 是矩阵元的类型，Options_ 只能是 Eigen::ColMajor 或 Eigen::RowMajor，默认是 ColMajor。StorageIndex_ 是内部数组的 index 的类型，必须是有符号整数，默认 int。注意这不是成员函数中行标和列表的 index 类型（Eigen::Index）。

　　 除非特殊说明，我们下面都以 RowMajor 为例。

内部结构

• Eigen 中 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 数组的指针

解线性方程组

• ConjugateGradient<SparseMatrix<Comp>, Lower|Upper> solver(H);
• solver.make(H)
• 要验证精度，可以用 (A*x - b).norm()/b.norm()。