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

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 行矩阵的三个运行。我们对 (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 节点。

编辑：5/10/2022 这是一个代码片段，它对列的子集执行矩阵乘法，其速度与任何不直接使用 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)然后将下部写入上部。另请参阅文档（此处和此处）。当然，如果你只能接受下半部分，则步骤(2)可以省略。

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

• 取 10 列时约为 1.1 倍，
• 获取 100 列时约为 1.5 倍，
• 获取 1000 列时约为 1.7 倍

在 Windows 上使用 MSVC，在 Linux 上使用 clang 并进行全面优化和 AVX。

启用并行化

另一种加速计算的方法是通过使用 OpenMP 编译来实现并行化。Eigen 负责剩下的事情。然而，上面利用对称性的代码并没有从中受益。但原来的代码

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

做。

对于 100x100000 矩阵mat，在 Linux 上使用 clang，使用 4 个线程（在 4 个真实核心上）运行并与单个线程进行比较，我得到的速度大约为

• 获取 10 列时约为 1.0 倍，即根本没有加速
• 获取 100 列时约为 1.8 倍
• 获取 1000 列时约为 2.0 倍

换句话说，除了极少量的列之外，4 个或更多核心的性能优于上面所示的对称方法。仅使用 2 个核心总是较慢。请注意，在我的测试中，使用SMT会降低性能，有时甚至非常明显。

其他注意事项

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

我认为在性能方面很难击败复制方法本身，因为特征矩阵默认是列主矩阵，这意味着复制几列相当便宜。此外，在不真正了解细节的情况下，我怀疑 Eigen 可以在整个矩阵上发挥其优化的全部力量来计算乘积并转置，而无需处理视图或类似的东西。这可能会给 Eigen 更多矢量化或缓存局部性的机会。

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