Эффективный по памяти метод для создания объекта dist из матрицы расстояний - PullRequest
Новая
       6

Эффективный по памяти метод для создания объекта dist из матрицы расстояний

2 голосов
/ 25 января 2020

Я пытаюсь создать dist объекты из матриц большого расстояния. Я исчерпал память, используя stats::as.dist. Например, на текущей машине у меня доступно около 128 ГБ, но as.dist не хватает памяти при обработке матрицы 73 000 x 73 000 (около 42 ГБ). Учитывая, что конечный объект dist должен быть меньше половины размера матрицы (т. Е. Это нижний треугольник, хранящийся в виде вектора), мне кажется, что можно сделать этот расчет более эффективным с точки зрения памяти. путь - при условии, что мы будем осторожны, чтобы не создавать большие промежуточные объекты, а просто скопировать соответствующие элементы ввода непосредственно в вывод.

Глядя на код для getS3method('as.dist', 'default'), я вижу, что он выполняет вычисления с использованием ans <- m[row(m) > col(m)], что требует создания матриц row и col с той же размерностью, что и для ввода.

Я подумал, что смогу улучшить это, используя алгоритм из здесь сгенерировать индексы нижнего треугольника. Вот моя попытка использовать этот метод. 

as.dist.new = function(dm, diag = FALSE, upper = FALSE) {
  n = dim(dm)[1]
  stopifnot(is.matrix(dm))
  stopifnot(dim(dm)[2] == n)
  k = 1:((n^2 - n)/2)
  j <- floor(((2 * n + 1) - sqrt((2 * n - 1) ^ 2 - 8 * (k - 1))) / 2)
  i <- j + k - (2 * n - j) * (j - 1) / 2
  idx = cbind(i,j)
  remove(i,j,k)
  gc()
  d = dm[idx]

  class(d) <- "dist"
  attr(d, "Size") <- n
  attr(d, "call") <- match.call()
  attr(d, "Diag") <- diag
  attr(d, "Upper") <- upper
  d
}

Это работает отлично на меньших матрицах. Вот простой пример: 

N = 10
dm = matrix(runif(N*N), N, N) 
diag(dm) = 0

x = as.dist(dm)
y = as.dist.new(dm)

Однако, если мы создаем матрицу большего расстояния, она сталкивается с теми же проблемами с памятью, что и as.dist.

Например, эта версия дает сбой в моей системе: 

N = 73000
dm = matrix(runif(N*N), N, N)
gc() 
diag(dm) = 0
gc()

as.dist.new(dm)

У кого-нибудь есть предложения, как выполнить эту операцию более эффективно? R или R cpp решения приветствуются. NB, глядя на этот ответ на связанную проблему (генерирование матрицы полного расстояния из данных местоположения из 2 столбцов), кажется, что возможно сделать это, используя RcppArmadillo, но у меня нет опыта использования этого.

Ответы [ 2 ]

1 голос
/ 25 января 2020

Что ж, логика c для обхода нижнего трима angular относительно проста, и если вы делаете это в C ++, то она также может быть быстрой: 

library(Rcpp)

sourceCpp(code='
// [[Rcpp::plugins(cpp11)]]

#include <cstddef> // size_t

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector as_dist(const NumericMatrix& mat) {
  std::size_t nrow = mat.nrow();
  std::size_t ncol = mat.ncol();
  std::size_t size = nrow * (nrow - 1) / 2;
  NumericVector ans(size);

  if (nrow > 1) {
    std::size_t k = 0;
    for (std::size_t j = 0; j < ncol; j++) {
      for (std::size_t i = j + 1; i < nrow; i++) {
        ans[k++] = mat(i,j);
      }
    }
  }

  ans.attr("class") = "dist";
  ans.attr("Size") = nrow;
  ans.attr("Diag") = false;
  ans.attr("Upper") = false;
  return ans;
}
')

as_dist(matrix(1:100, 10, 10))
    1  2  3  4  5  6  7  8  9
2   2                        
3   3 13                     
4   4 14 24                  
5   5 15 25 35               
6   6 16 26 36 46            
7   7 17 27 37 47 57         
8   8 18 28 38 48 58 68      
9   9 19 29 39 49 59 69 79   
10 10 20 30 40 50 60 70 80 90
0 голосов
/ 26 января 2020

Мое текущее решение состоит в том, чтобы вычислять объект dist непосредственно по широту и долготе без генерации матрицы промежуточного расстояния вообще. На больших матрицах это в несколько сотен раз быстрее, чем «обычный» конвейер с geosphere::mdist(), за которым следует stats::as.dist(), и использует только столько памяти, сколько требуется для хранения конечного объекта dist.

Следующий R cpp source основан на использовании функций из здесь для вычисления расстояния боярышника в c ++ вместе с адаптацией метода @ Alexis для перебора элементов нижнего треугольника в c ++. 

#include <Rcpp.h>
using namespace Rcpp;

double inverseHaversine(double d){
  return 2 * atan2(sqrt(d), sqrt(1 - d)) * 6378137.0;
}

double distanceHaversine(double latf, double lonf, double latt, double lont, double tolerance){
  double d;
  double dlat = latt - latf;
  double dlon =  lont - lonf;

  d = (sin(dlat * 0.5) * sin(dlat * 0.5)) + (cos(latf) * cos(latt)) * (sin(dlon * 0.5) * sin(dlon * 0.5));
  if(d > 1 && d <= tolerance){
    d = 1;
  }
  return inverseHaversine(d);
}

double toRadians(double deg){
  return deg * 0.01745329251;  // PI / 180;
}

//-----------------------------------------------------------
// [[Rcpp::export]]
NumericVector calc_dist(Rcpp::NumericVector lat, 
                        Rcpp::NumericVector lon, 
                        double tolerance = 10000000000.0) {
  std::size_t nlat = lat.size();
  std::size_t nlon = lon.size();
  if (nlat != nlon) throw std::range_error("lat and lon different lengths");
  if (nlat < 2) throw std::range_error("Need at least 2 points");
  std::size_t size = nlat * (nlat - 1) / 2;
  NumericVector ans(size);
  std::size_t k = 0;
  double latf;
  double latt;
  double lonf;
  double lont;

  for (std::size_t j = 0; j < (nlat-1); j++) {
    for (std::size_t i = j + 1; i < nlat; i++) {
      latf = toRadians(lat[i]);
      lonf = toRadians(lon[i]);
      latt = toRadians(lat[j]);
      lont = toRadians(lon[j]);
      ans[k++] = distanceHaversine(latf, lonf, latt, lont, tolerance);
    }
  }

  return ans;
}

/*** R
as_dist = function(lat, lon, tolerance  = 10000000000.0) {
  dd = calc_dist(lat, lon, tolerance)
  attr(dd, "class") = "dist"
  attr(dd, "Size") = length(lat)
  attr(dd, "call") = match.call()
  attr(dd, "Diag") = FALSE
  attr(dd, "Upper") = FALSE
  return(dd)
}
*/
...