Eigen 学习笔记

Eigen | 学习笔记

本文将记录使用Eigen重构机器学习模型中所用到的trick

1. Eigen 矩阵与 C++ 数组之间相互转换

#include <Eigen/Core>
#include <Eigen/Dense>

typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> EigenMatrixCol;
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> EigenMatrixRow;

// 数组 -> matrix (使用map实现)
// matrix在存储时都需要连续的内存空间，因此多维矩阵需要序列化成类似一维数组的形式，因此两者之间在某种程度上是等价的
// 多维矩阵在序列化时需要考虑行优先(row-major-storage)/列优先(column-major-storage)
// 行优先矩阵从数组中读取数据时以填充行数据的方式分割数组/列优先矩阵从数组中读取数据时以填充列数据的方式分割数组
// Eigen默认列优先, 为了运行方便，我们可以指定m1 * m2中m1为行优先m2为列优先,以方便矩阵乘法的运算
// https://eigen.tuxfamily.org/dox/group__TopicStorageOrders.html
double feat[8] = {0.01, 0.02, 0.03, 0.04};
EigenMatrixRow mat1 = Eigen::Map<Eigen::Matrix<double, 1, 4, Eigen::RowMajor>>(feat); // 注意(1, M)矩阵必须行优先
// 0.01, 0.02, 0.03, 0.04
EigenMatrixCol mat2 = Eigen::Map<Eigen::Matrix<double, 4, 1, Eigen::ColMajor>>(feat); // 注意(M, 1)矩阵必须列优先
// 0.01
// 0.02
// 0.03
// 0.04
EigenMatrixRow mat3 = Eigen::Map<Eigen::Matrix<double, 2, 2, Eigen::RowMajor>>(feat); // (N, M)矩阵行列均可
// 0.01 0.02
// 0.03 0.04
EigenMatrixCol mat4 = Eigen::Map<Eigen::Matrix<double, 2, 2, Eigen::ColMajor>>(feat);
// 0.01 0.03
// 0.02 0.04


// matrix -> 数组 (由属性.data实现)
// 数组必须是指针形式，不能 double feat2[4] = mat.data()
// 注意mat5和mat4在数值上虽然相等，但是存储形式已经发生变化
double* feat1 = mat1.data();
// 0.01 0.02 0.03 0.04
double* feat2 = mat2.data();
// 0.01 0.02 0.03 0.04
double* feat3 = mat3.data();
// 0.01 0.02 0.03 0.04
double* feat4 = mat4.data();
// 0.01 0.02 0.03 0.04
EigenMatrixRow mat5 = mat4;
double* feat5 = mat5.data();
// 0.01 0.03 0.02 0.04


// 常量 matrix 的初始化可以不借助数组, 直接 << 数值
EigenMatrixRow mat6(1, 3);
mat6 << -1.7476362058666228, 0.20625733243917574, 1.5413788734504181;
EigenMatrixRow mat7(7, 3);
mat7 <<   0.07596570881290864,  -0.029883796842660855, -0.046081911995243884,
            0.03717304278052893,  -0.36639910351580013,   0.32922606073528116,
            0.6881423398669501,    0.18600424502126356,  -0.8741465848882715,
            1.5223520003302156,   -0.24061676239120314,  -1.2817352379390454,
            -0.19630854550622312,   0.2059763534745924,   -0.009667807968369978,
            0.049970871859524445, -0.17010060664191579,   0.12012973478220677,
            0.48210908504659583,   0.18456527076604617,  -0.6666743558145485;

2. Eigen 矩阵乘法和元素乘法

Eigen中为了区分矩阵乘法和元素乘法，特别定义了两种不同的矩阵类型 Eigen::Matrix 和 Eigen::Array, 两个matrix之间的*被重载为矩阵乘法，两个array之间的*被重载为元素乘法。

#include <Eigen/Core>
#include <Eigen/Dense>

typedef Eigen::Matrix<float, 4, 4, Eigen::ColMajor> EigenMatrix44f;
typedef Eigen::Array<float, 4, 4, Eigen::ColMajor> EigenArray44f;

EigenArray44f a1, a2;
EigenMatrix44f m1, m2;
m1 = a1 * a2;                     // coeffwise product, implicit conversion from array to matrix.
a1 = m1 * m2;                     // matrix product, implicit conversion from matrix to array.
a2 = a1 + m1.array();             // mixing array and matrix is forbidden
m2 = a1.matrix() + m1;            // and explicit conversion is required.
Eigen::ArrayWrapper<EigenMatrix44f> m1a(m1);   // m1a is an alias for m1.array(), they share the same coefficients
Eigen::MatrixWrapper<EigenArray44f> a1m(a1);  

// more info can be found at https://eigen.tuxfamily.org/dox/group__QuickRefPage.html

3. Eigen Softmax 实现

softmax标准公式: \(\sigma (\mathbf {z} )_{j}={\frac {e^{z_{j}}}{\sum _{k=1}^{K}e^{z_{k}}}}\)

一般为了数值的稳定性，上式可以改写为: \(log [\sigma (\mathbf {z} )_{j}] = log[{\frac {e^{z_{j}-max(z)}}{\sum _{k=1}^{K}e^{z_{k}-max(z)}}}] = z_{j}-max(z)-log [{\sum _{k=1}^{K}e^{z_{k}-max(z)}}]\)

void EMatrixSoftmax(const EigenMatrixRow& input, EigenMatrixRow& output) {
    // using arrays allows to call native vectorizable exp and log
    EigenArray wMinusMax = input.colwise() - input.rowwise().maxCoeff();
    output = (wMinusMax.colwise() - wMinusMax.exp().rowwise().sum().log()).exp();
}