r - 在 R 中使用 gstat 进行克里金法时，多慢算太慢

Question

我正在尝试使用 R 的 gstat 包中的 krige 函数来插入 R 中的一些空间海洋深度数据。我发现超过约 1000点，该函数开始花费不合理的时间来完成(即，几小时到几天尚未完成)。这是正常现象还是我做错了什么？我特别担心，因为我的最终目标是对非常大的数据集(> 30,000 个数据点)进行时空克里金法，并且我担心考虑到这些运行时间，这是不可行的。

我正在运行 gstat-1.1-3 和 R-3.3.2。下面是我正在运行的代码:

library(sp); library(raster); library(gstat)
v.utm # SpatialPointsDataFrame with >30,000 points

# Remove points with identical positons
zd = zerodist(v.utm)
nzd = v.utm[-zd[,1],] # Layer with no identical positions

# Make a raster layer covering point layer
resolution=1e4
e = extent(as.matrix(v.utm@coords))+resolution
r = raster(e,resolution=resolution) 
proj4string(r) = proj4string(v.utm)

# r is a 181x157 raster

# Fit variogram
fv = fit.variogram(variogram(AVGDEPTH~1, nzd),model=vgm(6000,"Exp",1,5e5,1))

# Krige on random sample of 500 points - works fine
size=500
ss=nzd[sample.int(nrow(nzd),size),]
depth.krig = krige(AVGDEPTH~1,ss,as(r,"SpatialPixelsDataFrame"),
               model=depth.fit)

# Krige on random sample of 5000 points - never seems to end
size=5000
ss=nzd[sample.int(nrow(nzd),size),]
depth.krig = krige(AVGDEPTH~1,ss,as(r,"SpatialPixelsDataFrame"),
               model=depth.fit)

Answer 1

Choleski 分解(或类似分解)的复杂度为 O(n^3)，这意味着如果将点数乘以 10，则所需时间将增加 1000 倍。有两种方法可以解决这个问题，至少就 gstat 而言:

1. 安装 BLAS 的优化版本(例如 OpenBLAS 或 MKL) - 这不能解决 O(n^3) 问题，但可以最大程度地加快 n 倍的速度(其中 n 为可用内核数)
2. 避免通过选择局部邻域(参数 maxdist 和/或 nmax)来分解完整协方差矩阵