as.matrix на объекте расстояния чрезвычайно медленный; как сделать это быстрее?

Я нашел пакет R Rlof который использует многопоточность для вычисления матрицы расстояний, и это делает замечательную работу.

Тем не менее, вывод функции distmc это вектор, а не матрица. применение as.matrix Для этого "дистального" объекта получается гораздо дороже, чем многопоточный расчет расстояний.

Глядя на справку по функции, есть варианты печати диагонали и верхнего треугольника, но я не понимаю, где они должны использоваться.

Есть ли способ сэкономить время as.matrix каким-то образом?

Воспроизводимый пример:

set.seed(42)
M1 = matrix(rnorm(15000*20), nrow = 15000, ncol =20)
system.time({dA = distmc(M1, method = "euclidean", diag = TRUE,
                         upper = TRUE, p = 2)})
system.time(A = as.matrix(dA))
Источник

1 ответ

Решение

Что значит dist вернуть?

Эта функция всегда возвращает вектор, содержащий нижнюю треугольную часть (по столбцу) полной матрицы. diag или же upper аргумент влияет только на печать, т. е. stats:::print.dist, Эта функция может печатать этот вектор как матрицу; но это действительно не так.

Почему as.matrix плохо для "dist" объекта?

Трудно эффективно работать с треугольными матрицами или сделать их симметричными в R ядре. lower.tri а также upper.tri может быть занимать память, если ваша матрица велика: R: преобразовать верхнюю треугольную часть матрицы в симметричную матрицу.

Преобразование объекта dist в матрицу хуже. Посмотрите на код stats:::as.matrix.dist:

function (x, ...) 
{
    size <- attr(x, "Size")
    df <- matrix(0, size, size)
    df[row(df) > col(df)] <- x
    df <- df + t(df)
    labels <- attr(x, "Labels")
    dimnames(df) <- if (is.null(labels)) 
    list(seq_len(size), seq_len(size))
    else list(labels, labels)
    df
}

Использование row, col а также t это кошмар. И финал "dimnames<-" генерирует еще один большой временный матричный объект. Когда память становится узким местом, все идет медленно.

Но нам все еще, возможно, нужна полная матрица, поскольку она проста в использовании.

Неловко то, что проще работать с полной матрицей, поэтому мы этого хотим. Рассмотрим этот пример: R - Как получить индексы строк и столбцов совпадающих элементов из матрицы расстояний. Операции сложны, если мы пытаемся избежать формирования полной матрицы.

Решение Rcpp

Я пишу функцию Rcpp dist2mat (см. исходный файл "dist2mat.cpp" в конце этого ответа).

Функция имеет два входа: объект dist x и (целочисленный) коэффициент блокировки кэша bf, Функция работает, сначала создав матрицу и заполнив ее нижнюю треугольную часть, затем скопировав нижнюю треугольную часть в верхнюю треугольную, чтобы сделать ее симметричной. Второй шаг - это типичная операция транспонирования, и для крупномасштабной матричной кэш-памяти стоит блокировать. Производительность должна быть нечувствительной к фактору блокировки кэша, если она не слишком мала или не слишком велика. 128 или 256, как правило, хороший выбор.

Это моя первая попытка с Rcpp. Я был программистом C, использующим традиционный интерфейс C. Но я уже собираюсь познакомиться с Rcpp. Учитывая, что вы не знаете, как писать скомпилированный код, вы, вероятно, также не знаете, как запускать функции Rcpp. Вам нужно

  1. устанавливать Rcpp пакет (не уверен, нужен ли вам Rtools, если вы работаете в Windows);
  2. скопируйте мой "dist2mat.cpp" в файл в текущем рабочем каталоге вашего R (вы можете получить его из getwd() в вашей R сессии). Файл.cpp - это простой текстовый файл, поэтому вы можете создавать, редактировать и сохранять его в любом текстовом редакторе.

Теперь давайте начнем витрину.

library(Rcpp)
sourceCpp("dist2mat.cpp")  ## this takes some time; be patient

## a simple test with `dist2mat`
set.seed(0)
x <- dist(matrix(runif(10), nrow = 5, dimnames = list(letters[1:5], NULL)))
A <- dist2mat(x, 128)  ## cache blocking factor = 128
A
#          a         b         c         d         e
#a 0.0000000 0.9401067 0.9095143 0.5618382 0.4275871
#b 0.9401067 0.0000000 0.1162289 0.3884722 0.6968296
#c 0.9095143 0.1162289 0.0000000 0.3476762 0.6220650
#d 0.5618382 0.3884722 0.3476762 0.0000000 0.3368478
#e 0.4275871 0.6968296 0.6220650 0.3368478 0.0000000

Полученная матрица сохраняет имена строк вашей исходной матрицы / фрейма данных, переданных в dist,

Вы можете настроить фактор блокировки кэша на вашем компьютере. Обратите внимание, что эффект блокировки кэша не проявляется для небольших матриц. Здесь я попробовал один 10000 х 10000.

## mimic a "dist" object without actually calling `dist`
n <- 10000
x <- structure(numeric(n * (n - 1) / 2), class = "dist", Size = n)

system.time(A <- dist2mat(x, 64))
#   user  system elapsed 
#  0.676   0.424   1.113 

system.time(A <- dist2mat(x, 128))
#   user  system elapsed 
#  0.532   0.140   0.672 

system.time(A <- dist2mat(x, 256))
#   user  system elapsed 
#  0.616   0.140   0.759

Мы можем измерять dist2mat с as.matrix, Как as.matrix занимает много оперативной памяти, здесь я приведу небольшой пример.

## mimic a "dist" object without actually calling `dist`
n <- 2000
x <- structure(numeric(n * (n - 1) / 2), class = "dist", Size = n)

library(bench)
bench::mark(dist2mat(x, 128), as.matrix(x), check = FALSE)
## A tibble: 2 x 14
#  expression         min   mean  median     max `itr/sec` mem_alloc  n_gc n_itr
#  <chr>         <bch:tm> <bch:> <bch:t> <bch:t>     <dbl> <bch:byt> <dbl> <int>
#1 dist2mat(x, …   24.6ms   26ms  25.8ms  37.1ms     38.4     30.5MB     0    20
#2 as.matrix(x)   154.5ms  155ms 154.8ms 154.9ms      6.46   160.3MB     0     4
## ... with 5 more variables: total_time <bch:tm>, result <list>, memory <list>,
##   time <list>, gc <list>

Обратите внимание, как dist2mat быстрее (см. "mean", "median"), а также насколько меньше оперативной памяти (см. "mem_alloc") ему нужно. Я поставил check = FALSE отключить проверку результатов, потому что dist2mat не возвращает атрибут "dimnames" (у объекта "dist" такой информации нет), но as.matrix делает (это устанавливает 1:2000 как "dimnames"), поэтому они не совсем равны. Но вы можете убедиться, что они оба верны.

A <- dist2mat(x, 128)
B <- as.matrix(x)
range(A - B)
#[1] 0 0

"Dist2mat.cpp"

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix dist2mat(NumericVector& x, int bf) {

  /* input validation */
  if (!x.inherits("dist")) stop("Input must be a 'dist' object");
  int n = x.attr("Size");
  if (n > 65536) stop("R cannot create a square matrix larger than 65536 x 65536");

  /* initialization */
  NumericMatrix A(n, n);

  /* use pointers */
  size_t j, i, jj, ni, nj; double *ptr_x = &x[0];
  double *A_jj, *A_ij, *A_ji, *col, *row, *end;

  /* fill in lower triangular part */
  for (j = 0; j < n; j++) {
    col = &A(j + 1, j); end = &A(n, j);
    while (col < end) *col++ = *ptr_x++;
    }

  /* cache blocking factor */
  size_t b = (size_t)bf;

  /* copy lower triangular to upper triangular; cache blocking applied */
  for (j = 0; j < n; j += b) {
    nj = n - j; if (nj > b) nj = b;
    /* diagonal block has size nj x nj */
    A_jj = &A(j, j);
    for (jj = nj - 1; jj > 0; jj--, A_jj += n + 1) {
      /* copy a column segment to a row segment */
      col = A_jj + 1; row = A_jj + n;
      for (end = col + jj; col < end; col++, row += n) *row = *col;
      }
    /* off-diagonal blocks */
    for (i = j + nj; i < n; i += b) {
      ni = n - i; if (ni > b) ni = b;
      /* off-diagonal block has size ni x nj */
      A_ij = &A(i, j); A_ji = &A(j, i);
      for (jj = 0; jj < nj; jj++) {
        /* copy a column segment to a row segment */
        col = A_ij + jj * n; row = A_ji + jj;
        for (end = col + ni; col < end; col++, row += n) *row = *col;
        }
      }
    }

  /* add row names and column names */
  A.attr("dimnames") = List::create(x.attr("Labels"), x.attr("Labels"));

  return A;
  }
Источник
Другие вопросы по тегам