Increase the Speed of Krig2D
QUESTION: I've heard on the grapevine that it is possible to increase the speed of Krig2D by an order of magnitude. Is this true?
ANSWER: Apparently so. Michele Cappellari discovered a mathematically equivalent expression in Krig2D that speeds this extremely slow function up by at least an order of magnitude. As of IDL 8.2.3, you should findthe following lines of code in the krig2d.pro file (found in the IDL lib directory).
for j=0,ny-1 do begin ;Each output point
y0 = bounds[1] + gs[1] * j
for i=0,nx-1 do begin
x0 = bounds[0] + gs[0] * i
d[0] = reform(sqrt((x-x0)^2 + (y-y0)^2),n) ;distance
d = call_function(fname, d, t) ;Get rhs
d[n] = 1.0 ;lagrange constr
d=LUSOL( a, indx, d) ;Solve using LUDC results.
r[i,j] = total(d * z)
endfor
endfor
Replace this code with the following lines, which essentially just moves the LUSOL call outside the double FOR loop and changes the order of the array operations. According to Michele, this works because array multiplication is cummutative.
az = LUSOL(a, indx, [z,0.0]) ;Solve using LUDC results.
for j=0,ny-1 do begin ;Each output point
y0 = bounds[1] + gs[1] * j
for i=0,nx-1 do begin
x0 = bounds[0] + gs[0] * i
d[0] = reform(sqrt((x-x0)^2 + (y-y0)^2),n) ;distance
d = call_function(fname, d, t) ;Get rhs
d[n] = 1.0 ;lagrange constr
r[i,j] = total(d * az)
endfor
endfor
Michele supplies the following example routine, which is more than 30 times faster than the unmodified code and more than seven times faster than the built-in (and presumably optimized) GridData code, which suggests that GridData is making the same mistake and can be made faster, too, although not by users.
; Create a dataset of N points.
n = 500 ;# of scattered points
seed = -121147L ;For consistency
x = RANDOMU(seed, n)
y = RANDOMU(seed, n)
; Create a dependent variable in the form a function of (x,y)
z = 3 * EXP(-((9*x-2)^2 + (7-9*y)^2)/4) + $
3 * EXP(-((9*x+1)^2)/49 - (1-0.9*y)) + $
2 * EXP(-((9*x-7)^2 + (6-9*y)^2)/4) - $
EXP(-(9*x-4)^2 - (2-9*y)^2)
s = [0.5,0]
tic
z2 = krig2d_modified(z, x, y, EXPONENTIAL=s, NX=200, NY=200)
toc
tic
z2 = GRIDDATA(x, y, z, /KRIGING, DIMENSION=[200,200], VARIOGRAM=[2,s])
toc
tic
z2 = krig2d(z, x, y, EXPONENTIAL=s, NX=200, NY=200)
toc
END
Further Enhancements
The folks at ExelisVis got right on this, but they discovered a further problem. The Krig_Sphere function code in Krig2D doesn't match the documentation, and there appears to be a problem with the code itself.
Eventually, things got sorted out and ExelisVis posted a new Krig2D_Fast program on the IDL newsgroup. Running Michele's example code with the new program was about 80 times faster than the original IDL Krig2D program on my computer! Similar code has been added to the GridData program and will make that program much faster at Kriging in IDL 8.3. Here are the results as reported by ExelisVis:
Old KRIG2D: 30 seconds New KRIG2D: 0.51 seconds New GRIDDATA: 1.3 seconds
Note, too, that Michele reports the "new covariance functions now looks right and corrects for some major numerical differences one can see against the old Krig2D, when the data extend beyond the region where the semivariogram was supposed to flatten out."
 |
| The 2D result of the fast Kriging algorithm. |
Coyote Kriging Program
Because I wanted to have a kriging function that I could use with all versions of IDL, and because I wanted to learn more about kriging in general, I decided to write my own Krig2D function, which I called cgKrig2D. I decided in my function to use semivariograms rather than covariance matrices, and I modified the kriging models appropriately. I also added additional error handling and I discovered and fixed a problem with the original Krig2D command that was always confusing to me. Namely, whenever I used a 2D grid as input to the function, I seemed to get nonsense values back. I never understood why.
It turns out that when kriging a 2D grid that you must match the size of the input grid with the size of the output grid. This is not mentioned in the documentation. Here is an example of what I mean.
inputGrid = Dist(11) outputGrid = Krig2D(inputGrid, Spherical=[5.0, 0.0]) !P.Multi=[0,2,1] cgDisplay, 1200, 450 cgSurf, inputGrid cgSurf, outputGrid !P.Multi=0
You see the result of this code in the figure below.
 |
| If the dimensions of the input and output grid are not matched, chaos ensues. |
inputGrid = Dist(11) outputGrid = cgKrig2D(inputGrid, Spherical=[5.0, 0.0]) !P.Multi=[0,2,1] cgDisplay, 1200, 450 cgSurf, inputGrid cgSurf, outputGrid !P.Multi=0
You see the result of this code in the figure below.
 |
| The cgKrig2D program makes sure dimensions of gridded data are matched. |
Version of IDL used to prepare this article: IDL 8.2.3.
