列子集的特征矩阵的矩阵乘法
Matrix multiplication of an Eigen Matrix for a subset of columns
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，它仍然是 3.x 版本的 Eigen（不支持 () 索引数组），无论如何，我的理解是 () 运算符仍然执行深拷贝。
现在我正在做一个深拷贝并生成一个新矩阵，其中的数据只包含 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 列的子集进行矩阵乘法。
R:
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 日
这是一个代码片段，它对列的子集执行矩阵乘法的速度与不直接使用 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;
}
利用对称性
您可以利用生成的矩阵将像这样对称：
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 and here)。
当然，如果你只接受下半部分，步骤（2）可以省略。
对于 100x100000 矩阵 mat，我的速度提高了大约
- 取10列时~1.1x，
- 取100列时~1.5x，
- 取 1000 列时~1.7 倍
在 Windows 使用 MSVC 和在 Linux 使用 clang 以及完全优化和 AVX。
启用并行化
另一种加快计算速度的方法是通过使用 OpenMP 进行编译来启用 parallelization。 Eigen 负责剩下的事情。但是，上面利用对称性的代码 不会 从中受益。但是原码
Eigen::MatrixXd sub_mat = subset_cols(mat, idx);
Eigen::MatrixXd a = sub_mat * sub_mat.transpose();
会。
对于 100x100000 矩阵 mat，在 Linux、运行 4 个线程（在 4 个真正的内核上）上使用 clang 并与单个线程进行比较，我的速度提高了大约
的一个因数
- 取10列时~1.0x，即完全没有加速
- 取 100 列时~1.8x
- 取 1000 列时~2.0x
换句话说，除了极少数列外，4核或更多核的性能优于上面显示的对称方法。只使用 2 个内核总是比较慢。请注意，使用 SMT 会损害我的测试性能，有时尤其如此。
其他说明
我已经在评论中写了这个，但为了完整起见：
Eigen::Map 将不起作用，因为步幅为 non-equidistant。使用 slicing gives me ~10% better performance than your copying method on Linux with clang and gcc, but somewhat worse on MSVC. Also, as you noted, it is not available on the 3.3 branch of Eigen. There is a custom way 来模仿它，但它在我的测试中表现总是更差。
另外，在我的测试中，与复制方法相比，它没有节省任何内存。
我认为在性能方面很难击败复制方法本身，因为本征矩阵默认为 column major，这意味着复制几列的成本相当低。此外，在不真正了解细节的情况下，我怀疑 Eigen 可以将其优化的全部力量投入到完整矩阵上以计算乘积和转置，而无需处理视图或类似的东西。这可能会给 Eigen 更多机会进行矢量化或缓存局部性。
除此之外，不仅应该打开优化，还应该使用尽可能高的指令集。在我的测试中打开 AVX 将性能提高了约 1.5 倍。很遗憾，我无法测试 AVX512。
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，它仍然是 3.x 版本的 Eigen（不支持 () 索引数组），无论如何，我的理解是 () 运算符仍然执行深拷贝。
现在我正在做一个深拷贝并生成一个新矩阵，其中的数据只包含 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 列的子集进行矩阵乘法。
R:
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 日 这是一个代码片段，它对列的子集执行矩阵乘法的速度与不直接使用 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;
}
利用对称性
您可以利用生成的矩阵将像这样对称：
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 and here)。
当然，如果你只接受下半部分，步骤（2）可以省略。
对于 100x100000 矩阵 mat，我的速度提高了大约
- 取10列时~1.1x，
- 取100列时~1.5x，
- 取 1000 列时~1.7 倍
在 Windows 使用 MSVC 和在 Linux 使用 clang 以及完全优化和 AVX。
启用并行化
另一种加快计算速度的方法是通过使用 OpenMP 进行编译来启用 parallelization。 Eigen 负责剩下的事情。但是，上面利用对称性的代码 不会 从中受益。但是原码
Eigen::MatrixXd sub_mat = subset_cols(mat, idx);
Eigen::MatrixXd a = sub_mat * sub_mat.transpose();
会。
对于 100x100000 矩阵 mat，在 Linux、运行 4 个线程（在 4 个真正的内核上）上使用 clang 并与单个线程进行比较，我的速度提高了大约
- 取10列时~1.0x，即完全没有加速
- 取 100 列时~1.8x
- 取 1000 列时~2.0x
换句话说，除了极少数列外，4核或更多核的性能优于上面显示的对称方法。只使用 2 个内核总是比较慢。请注意，使用 SMT 会损害我的测试性能，有时尤其如此。
其他说明
我已经在评论中写了这个，但为了完整起见：
Eigen::Map 将不起作用，因为步幅为 non-equidistant。使用 slicing gives me ~10% better performance than your copying method on Linux with clang and gcc, but somewhat worse on MSVC. Also, as you noted, it is not available on the 3.3 branch of Eigen. There is a custom way 来模仿它，但它在我的测试中表现总是更差。
另外，在我的测试中，与复制方法相比，它没有节省任何内存。
我认为在性能方面很难击败复制方法本身，因为本征矩阵默认为 column major，这意味着复制几列的成本相当低。此外，在不真正了解细节的情况下，我怀疑 Eigen 可以将其优化的全部力量投入到完整矩阵上以计算乘积和转置，而无需处理视图或类似的东西。这可能会给 Eigen 更多机会进行矢量化或缓存局部性。
除此之外，不仅应该打开优化，还应该使用尽可能高的指令集。在我的测试中打开 AVX 将性能提高了约 1.5 倍。很遗憾，我无法测试 AVX512。