列子集的特征矩阵的矩阵乘法

Question

将Eigen::Matrix 与一组随机列索引进行矩阵乘法的最快方法是什么？

Eigen::MatrixXd mat = Eigen::MatrixXd::Random(100, 1000);
// vector of random indices (linspaced here for brevity)
Eigen::VectorXi idx = VectorXi::LinSpaced(8,1000,9);

我正在使用 RcppEigen 和 R，它仍然在 Eigen 的 3.x 版本上（不支持带有索引数组的 ()），无论如何，我的理解是 () 运算符仍然执行深层复制。

现在我正在做一个深拷贝并生成一个新矩阵，其中只包含idx 中的列的数据：

template <typename T>
inline Eigen::Matrix<T, -1, -1> subset_cols(const Eigen::Matrix<T, -1, -1>& x, const std::vector<size_t>& cols) {
    Eigen::Matrix<T, -1, -1> y(x.rows(), cols.size());
    for (size_t i = 0; i < cols.size(); ++i)
        y.col(i) = x.col(cols[i]);
    return y;
}

然后做矩阵乘法：

Eigen::MatrixXd sub_mat = subset_cols(mat, idx);
Eigen::MatrixXd a = sub_mat * sub_mat.transpose();

a 是我想要的。

一定有办法避免深拷贝，而是使用Eigen::Map?

22 年 5 月 9 日编辑：作为对@Markus 的回复，他提出了一种使用原始数据访问和Eigen::Map 的方法。所提出的解决方案比深拷贝的矩阵乘法要慢一些。这里的基准测试是使用 Rcpp 代码和 R 完成的：

//[[Rcpp::depends(RcppClock)]]
#include <RcppClock.h>

//[[Rcpp::export]]
void bench(Eigen::MatrixXd mat, Eigen::VectorXi idx){
  Rcpp::Clock clock;
  size_t reps = 100;
  while(reps-- > 0){
    clock.tick("copy");
    Eigen::MatrixXd sub_mat = subset_cols(mat, idx);
    Eigen::MatrixXd a = sub_mat * sub_mat.transpose();
    clock.tock("copy");
    
    clock.tick("map");
    double *b_raw = new double[mat.rows() * mat.rows()];
    Eigen::Map<Eigen::MatrixXd> b(b_raw, mat.rows(), mat.rows());
    subset_AAt(b_raw, mat, idx);
    clock.tock("map");
  }
  clock.stop("clock");
}

这是一个 100,000 列的 100 行矩阵的 3 次运行。我们正在对 (1) 10 列的子集、(2) 1000 列的子集和 (3) 10000 列的子集进行矩阵乘法。

回复：

bench(
  matrix(runif(100000 * 100), 100, 100000), 
  sample(100000, 10) - 1)

# Unit: microseconds 
# ticker   mean     sd   min    max neval
#    copy  31.65  4.376 30.15  69.46   100
#     map 113.46 21.355 68.54 166.29   100

bench(
  matrix(runif(100000 * 100), 100, 100000), 
  sample(100000, 1000) - 1)

#  Unit: milliseconds 
#  ticker  mean     sd   min   max neval
#    copy 2.361 0.5789 1.972  4.86   100
#     map 9.495 2.4201 7.962 19.90   100

bench(
  matrix(runif(100000 * 100), 100, 100000), 
  sample(100000, 10000) - 1)

#  Unit: milliseconds 
#  ticker   mean     sd    min   max neval
#    copy  23.04  2.774  20.95  42.4   100
#     map 378.14 19.424 351.56 492.0   100

我在几台机器上进行了基准测试，结果相似。以上结果来自一个好的 HPC 节点。

编辑：2022 年 5 月 10 日这是一个代码 sn-p，它对列子集执行矩阵乘法的速度与不直接使用 Eigen BLAS 的任何代码一样快：

template <typename T>
Eigen::Matrix<T, -1, -1> subset_AAt(const Eigen::Matrix<T, -1, -1>& A, const Eigen::VectorXi& cols) {
  const size_t n = A.rows();
  Eigen::Matrix<T, -1, -1> AAt(n, n);
  for (size_t k = 0; k < cols.size(); ++k) {
    const T* A_data = A.data() + cols(k) * n;
    for (size_t i = 0; i < n; ++i) {
      T tmp_i = A_data[i];
      for (size_t j = 0; j <= i; ++j) {
        AAt(i * n + j) += tmp_i * A_data[j];
      }
    }
  }
  return AAt;
}

Answer 1

利用对称性

您可以利用生成的矩阵将是对称的，如下所示：

Mat sub_mat = subset_cols(mat, idx); // From your original post
Mat a = Mat::Zero(numRows, numRows);
a.selfadjointView<Eigen::Lower>().rankUpdate(sub_mat); // (1)
a.triangularView<Eigen::Upper>() = a.transpose(); // (2)

(1) 行将仅计算下部的 a += sub_mat * sub_mat.transpose()。 (2) 然后将下部写入上部。另请参阅文档（here 和 here）。 当然，如果你只能忍受下半部分，步骤（2）可以省略。

对于 100x100000 矩阵 mat，我得到了大约 1 倍的加速

• ~1.1x 时取 10 列，
• ~1.5x 时取 100 列，
• ~1.7x 时取 1000 列

在 Windows 上使用 MSVC 和在 Linux 上使用完全优化和 AVX 的 clang。

启用并行化

加快计算速度的另一种方法是通过使用 OpenMP 编译来启用parallelization。 Eigen 负责其余的工作。上面利用对称性的代码确实不是然而，从中受益。但是原来的代码

Eigen::MatrixXd sub_mat = subset_cols(mat, idx);
Eigen::MatrixXd a = sub_mat * sub_mat.transpose();

做。

对于 100x100000 矩阵 mat，在 Linux 上使用 clang，运行 4 个线程（在 4 个真实内核上）并与单线程相比，我得到了大约 1 倍的加速

• ~1.0x 时占用 10 列，即根本没有加速
• ~1.8x 时取 100 列
• ~2.0x 时取 1000 列

换句话说，除了极少数列外，4 个或更多内核的性能优于上面显示的对称方法。仅使用 2 个内核总是较慢。请注意，使用SMT 会损害我的测试中的性能，有时很明显。

其他注意事项

我已经在评论中写了这个，但为了完整起见： Eigen::Map 将不起作用，因为步幅是非等距的。使用slicing 给我的性能比在 Linux 上使用 clang 和 gcc 的复制方法好约 10%，但在 MSVC 上稍差。此外，正如您所指出的，它在 Eigen 的 3.3 分支上不可用。有一个 custom way 可以模仿它，但在我的测试中它的表现总是更差。 此外，在我的测试中，与复制方法相比，它并没有节省任何内存。

我认为在性能方面很难超越复制方法本身，因为特征矩阵默认为column major，这意味着复制几列相当便宜。此外，在没有真正了解细节的情况下，我怀疑 Eigen 可以将其优化的全部力量投入到完整矩阵上来计算乘积和转置，而无需处理视图或类似的事情。这可能会给 Eigen 更多的机会进行矢量化或缓存局部性。

除此之外，不仅应该打开优化，而且应该使用最高可能的指令集。在我的测试中打开 AVX 将性能提高了约 1.5 倍。不幸的是，我无法测试 AVX512。

Answer 2

万一有人在以后发现这很有帮助，我可以使用 OpenMP 和三角索引在已接受的问题中击败 Eigen 代码的性能。在这种情况下，我使用Rcpp::NumericMatrix，但您可以将Eigen::MatrixXd 插入：

    Rcpp::NumericMatrix Rcpp_AAt(const Rcpp::NumericMatrix& mat) {
    const size_t n = mat.cols();
    const size_t n_vals = n / 2 * (1 + n) - n;
    Rcpp::NumericMatrix res(n, n);
    #pragma omp parallel for
    for (size_t k = 0; k < (n_vals + n); ++k) {
        // k is linear index
        if (k >= n_vals) {
            size_t i = k - n_vals;
            double tmp = 0;
            for (size_t row = 0; row < mat.rows(); ++row)
                tmp += mat(row, i) * mat(row, i);
            res(i, i) = tmp;
        } else {
            size_t i = n - 2 - std::floor(std::sqrt(-8 * k + 4 * n * (n - 1) - 7) / 2.0 - 0.5);
            size_t j = k + i + 1 - n * (n - 1) / 2 + (n - i) * ((n - i) - 1) / 2;
            double tmp = 0;
            for (size_t row = 0; row < mat.rows(); ++row)
                tmp += mat(row, i) * mat(row, j);
            res(i, j) = tmp;
            res(j, i) = tmp;
        }
    }
    return res;
}

通过使用三角索引，我们允许 OpenMP 为所有列组合生成线程，这比一次仅在一列上并行（出于显而易见的原因）更有效。 Eigen 使用多线程，所以我认为这是公平的游戏。