距离物体上的矩阵非常慢;如何使它更快?

做什么 dist 回来吗?

这个函数总是返回一个向量,包含整个矩阵的下三角部分(按列)。这个 diag 或 upper 参数只影响打印,即。, stats:::print.dist . 这个函数可以像打印矩阵一样打印这个向量,但实际上不是。

为什么 as.matrix 对“dist”对象有害?

在R核中很难有效地处理三角矩阵或进一步使其对称。 lower.tri 和 upper.tri 如果矩阵很大,可能会消耗内存: R: Convert upper triangular part of a matrix to symmetric matrix .

将“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 - How to get row & column subscripts of matched elements from a distance matrix . 如果我们试图避免形成完整的矩阵,那么操作是很棘手的。

一个 Rcpp公司 解决方案

我编写了一个Rcpp函数 dist2mat (请参阅本答案末尾的“dist2mat.cpp”源文件)。

函数有两个输入:“dist”对象 x 和(整数)缓存阻塞因子 bf . 该函数首先创建一个矩阵并填充其下三角部分,然后将下三角部分复制到上三角部分以使其对称。第二步是典型的换位操作,对于大型矩阵缓存阻塞是值得的。性能应该对缓存阻塞因子不敏感,只要它不太小或太大。128或256通常是一个不错的选择。

这是我第一次尝试Rcpp。我是一个使用R的传统C接口的C程序员。但我也在去熟悉Rcpp的路上。鉴于您不知道如何编写编译后的代码,您可能也不知道如何运行Rcpp函数。你需要

1. 安装 Rcpp 包装(不确定是否需要 Rtools 如果你在窗户上);
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

结果矩阵保留传递给的原始矩阵/数据帧的行名 距离 .

您可以在您的计算机上调整缓存阻塞因子。注意,缓存阻塞对小矩阵的影响并不明显。我试了一个10000 x 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 

我们可以设定基准 距离2英里 具有 as.矩阵 . 作为 as.矩阵 是内存消耗,我这里用一个小例子。

## 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>

注意如何 距离2英里 速度更快(参见“平均值”、“中值”),而且所需的RAM(参见“内存分配”)也更少。我已经准备好了 check = FALSE 禁用结果检查,因为 距离2英里 不返回“dimnames”属性(“dist”对象没有这样的信息),但是 as.矩阵 是的 1:2000 作为“dimnames”),所以它们并不完全相等。但你可以确认它们都是正确的。

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

“距离2mat.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;
  }