对于一个随机的列索引集,对 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
?
编辑 5/9/22:
针对@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行、100,000列矩阵的三个运算实例。我们将对(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;
}
Eigen::Map
不起作用,因为步幅是非等距的。在Linux上使用slicling比你的subset_cols()
方式性能提高了约10%,在clang和gcc上表现更好,但在MSVC上表现较差。正如你所指出的,它在3.3分支上不可用。有一种自定义的方法可以模仿它,但在我的测试中始终表现较差。启用AVX可以获得最佳改进(快约1.5倍),甚至可以启用AVX512。 - Sedenion