сборка матрицы из диагональных срезов с помощью mclapply или% dopar%, например Matrix :: bandSparse - PullRequest
Новая
       34

сборка матрицы из диагональных срезов с помощью mclapply или% dopar%, например Matrix :: bandSparse

0 голосов
/ 03 мая 2018

Прямо сейчас я работаю с некоторыми огромными матрицами в 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, которая распараллелена.

1 Ответ

0 голосов
/ 04 мая 2018

Если вы хотите построить одну матрицу, вы ищете разделяемую память параллелизм. И 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.

...