Время, затраченное на krige в пакете gstat в R

Question

Следующая программа R создает интерполированную поверхность, используя 470 точек данных, используя данные о методе Walker Lake в пакете gstat.

source("D:/kriging/allfunctions.r")   # Reads in all functions. 
source("D:/kriging/panel.gamma0.r")   # Reads in panel function for xyplot. 
library(lattice)       # Needed for "xyplot" function. 
library(geoR)        # Needed for "polygrid" function. 
library(akima) 
library(gstat); 
library(sp); 
walk470 <- read.table("D:/kriging/walk470.txt",header=T) 
attach(walk470) 
coordinates(walk470) = ~x+y 
walk.var1 <- variogram(v ~ x+y,data=walk470,width=10) #the width has to be tuned resulting different point pairs 
plot(walk.var1,xlab="Distance",ylab="Semivariance",main="Variogram for V, Lag Spacing = 5") 
model1.out <- fit.variogram(walk.var1,vgm(70000,"Sph",40,20000)) 
plot(walk.var1, model=model1.out,xlab="Distance",ylab="Semivariance",main="Variogram for V, Lag Spacing = 10") 
poly <- chull(coordinates(walk470)) 
plot(coordinates(walk470),type="n",xlab="X",ylab="Y",cex.lab=1.6,main="Plot of Sample and Prediction Sites",cex.axis=1.5,cex.main=1.6) 
lines(coordinates(walk470)[poly,]) 
poly.in <- polygrid(seq(2.5,247.5,5),seq(2.5,297.5,5),coordinates(walk470)[poly,]) 
points(poly.in) 
points(coordinates(walk470),pch=16) 
coordinates(poly.in) <- ~ x+y 
krige.out <- krige(v ~ 1, walk470,poly.in, model=model1.out) 
print(krige.out) 

Эта программа вычисляет следующую информацию по каждой точке 2688 точек

(470x470) matrix inversion 
(470x470) and (470x1) matrix multiplication 

Является GSTAT пакет использует некоторые умный способ расчета. Я знал из предыдущего запроса stackoverflow, что он использует cholesky-разложение для инверсии матрицы. Нормальная скорость для одной машины рассчитана так быстро.

Answer 1

Он использует декомпозицию LDL, которая похожа на Choleski. Когда вы используете глобальное кригинг, матрицу ковариации нужно разложить только один раз; то для каждой точки предсказания решается система, которая является O (n). Матрица 470x470 никогда не перевернута, ни решения, полученные путем ее умножения. Инверсы - это нотные устройства, но, по возможности, избегаются как вычислительная стратегия. В R, например, сравните время выполнения solve(A,b) с solve(A) %*% b.

Используйте source, Luke!