使用 mclapply 或 %dopar% 从对角线切片组装矩阵，例如 Matrix::bandSparse

Question

现在我正在 R 中处理一些巨大的矩阵，我需要能够使用对角带重新assemble它们。出于编程原因（为了避免必须对大小为 n 的矩阵执行 n*n 操作（数百万次计算），我只想进行 2n 次计算（数千次计算），因此选择运行我的功能矩阵的对角线带。现在，我得到了结果，但需要采用这些矩阵切片并 assemble 它们以允许我使用多个处理器的方式。

foreach 和 mclapply 都不允许我修改循环外的对象，所以我正在尝试考虑一个并行解决方案。如果有一些函数可以将非对角带分配给可以可靠地完成的矩阵的一部分，我完全赞成。

输入：

[1] 0.3503037

[1] 0.2851895 0.2851895

[1] 0.5233396 0.5233396 0.5233396

[1] 0.6250584 0.6250584 0.6250584 0.6250584

[1] 0.4300964 0.4300964 0.4300964 0.4300964 0.4300964

[1] 0.4300964 0.4300964 0.4300964 0.4300964 0.4300964

[1] 0.3949782 0.3949782 0.3949782 0.3949782

[1] 0.7852812 0.7852812 0.7852812

[1] 0.5309648 0.5309648

[1] 0.7718504

期望的输出（并行操作）：

          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.4300964 0.6250584 0.5233396 0.2851895 0.3503037

[2,] 0.3949782 0.4300964 0.6250584 0.5233396 0.2851895

[3,] 0.7852812 0.3949782 0.4300964 0.6250584 0.5233396

[4,] 0.5309648 0.7852812 0.3949782 0.4300964 0.6250584

[5,] 0.7718504 0.5309648 0.7852812 0.3949782 0.4300964

我越看这个，我需要一个并行化的 Matrix::bandSparse 版本。

Answer 1

如果要构建单个矩阵，您需要共享内存 并行性。 parallel 和 foreach 都实现了 分布式内存 并行。我知道一个实现共享内存的 R 包 (Rdsm)，但我没有使用过它。获得共享内存并行性的更自然的方法是使用 C++。

我已经在 R（串行）、C++ 和 Rcpp（串行）和 C++ 加上 OpenMP 和 Rcpp 和 RcppParallel（并行）中实现了波段到矩阵的转换。请注意，我使用的输入是一个没有重复对角线的向量列表。对于 OpenMP 解决方案，我将其转换为（参差不齐的）matrix，因为这可以轻松转换为线程安全的 RMatrix。即使在 R:

中，这种转换也非常快

#include <Rcpp.h>
#include <algorithm>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix diags2mtrCpp(int n, const ListOf<const NumericVector>& diags) {
  NumericMatrix mtr(n, n);
  int nDiags = diags.size();
  for (int i = 0; i < nDiags; ++i) {
    NumericVector diag(diags[i]);
    int nDiag = diag.size();
    int row = std::max(1, i - n + 2);
    int col = std::max(1, n - i);
    for (int j = 0; j < nDiag; ++j) {
      mtr(row + j - 1, col + j - 1) = diag(j);
    }
  }
  return mtr;
}

// [[Rcpp::plugins(openmp)]]
#include <omp.h>
// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
using namespace RcppParallel;

// [[Rcpp::export]]
NumericMatrix diags2mtrOmp(const NumericMatrix& diags_matrix, const IntegerVector& diags_length) {
  int nDiags = diags_matrix.cols();
  int n = diags_matrix.rows();
  NumericMatrix res(n, n);
  RMatrix<double> mtr(res);
  RMatrix<double> diags(diags_matrix);
  RVector<int> diagSize(diags_length);
  #pragma omp parallel for
  for (int i = 0; i < nDiags; ++i) {
    int nDiag = diagSize[i];
    int row = std::max(1, i - n + 2);
    int col = std::max(1, n - i);
    for (int j = 0; j < nDiag; ++j) {
      mtr(row + j - 1, col + j - 1) = diags(j, i);
    }
  }
  return res;
}


/*** R
set.seed(42)
n <- 2^12
n
diags <- vector(mode = "list", length = 2 * n - 1)
for (i in seq_len(n)) {
  diags[[i]] <- rep.int(runif(1), i)
  diags[[2 * n - i]] <- rep.int(runif(1), i)
}

diags_matrix <- matrix(0, nrow = n, ncol = length(diags))
diags_length <- integer(length(diags))
for (i in seq_along(diags)) {
  diags_length[i] <- length(diags[[i]])
  diags_matrix[ ,i] <- c(diags[[i]], rep.int(0, n - diags_length[i]))
}


diags2mtr <- function(n, diags) {
  mtr <- matrix(0, n, n)
  for (i in seq_along(diags)) {
    row <- max(1, i - n + 1)
    col <- max(1, n + 1 - i)
    for (j in seq_along(diags[[i]]))
      mtr[row + j - 1 , col + j - 1] <- diags[[i]][j]
  }
  mtr

}
system.time(mtr <- diags2mtr(n, diags))
system.time(mtrCpp <- diags2mtrCpp(n, diags))
system.time(mtrOmp <- diags2mtrOmp(diags_matrix, diags_length))
all.equal(mtr, mtrCpp)
all.equal(mtr, mtrOmp)
*/

在双核机器上对这些解决方案进行基准测试得到：

Unit: milliseconds
         expr        min        lq      mean    median        uq       max neval
    diags2mtr 2252.82538 2271.7221 2354.1251 2323.8221 2382.7958 2558.9282    10
 diags2mtrCpp  161.25920  190.9728  224.9094  226.2652  265.3675  279.3848    10
 diags2mtrOmp   95.50714  100.9555  105.8462  102.4064  105.7645  127.5200    10

我对 diags2mtrOmp 的加速感到惊讶。

使用 mclapply 或 %dopar% 从对角线切片组装矩阵，例如 Matrix::bandSparse

assembling a matrix from diagonal slices with mclapply or %dopar%, like Matrix::bandSparse

r

gpu

doparallel

mclapply

domc