Время, затрачиваемое на криг в пакете 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 использует какой-то умный способ для расчета. Из предыдущего запроса stackru я знал, что он использует разложение Холецкого для инверсии матрицы. Это нормальная скорость для одной машины, чтобы вычислить это так быстро.

Answer 1

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

Используйте источник, Люк!