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

Posted 2023-02-19

【中文标题】使用 mclapply 或 %dopar% 从对角线切片组装矩阵，例如 Matrix::bandSparse

【英文标题】：assembling a matrix from diagonal slices with mclapply or %dopar%, like Matrix::bandSparse

【发布时间】：2018-10-14 02:43:08

【问题描述】：

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

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 版本。

【问题讨论】：

你可能想研究像 Eigen 或 Armadillo 这样的库。

为什么输入中的对角线重复？矩阵是稀疏的还是密集的？

您的绩效目标是什么？对于 4096 x 4096 矩阵，我得到了 2.7 秒（串行 R）和 0.3 秒（通过 Rcpp 串行 C++）的运行时间。

【参考方案1】：

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

我已经在 R（串行）、带有 Rcpp（串行）的 C++ 和带有 Rcpp 和 RcppParallel（并行）的 C++ 以及 OpenMP 中实现了波段到矩阵的转换。请注意，我使用的输入是一个没有重复对角线的向量列表。对于 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 的加速感到惊讶。

【讨论】：

哇...这是一个了不起的解决方案。我从没想过使用 OMP 会基于 C 语言。

@JamesDalgleish 我更新了 OpenMP 解决方案，因为原始版本有时会由于内存访问问题而崩溃。我应该比从线程代码访问 R 数据结构更清楚......