[R-sig-Geo] 'LDLfactor' error in 'krige' function

Edzer Pebesma edzer.pebesma at uni-muenster.de
Tue Nov 17 13:00:54 CET 2009 

Tom, you can already do this:

library(gstat)
data(meuse)
coordinates(meuse)=~x+y
data(meuse.grid)
gridded(meuse.grid)=~x+y
meuse=meuse[c(1,1:155),] # replicate first observation
pr1 = predict(zinc~1,meuse,meuse.grid,vgm(1, "Exp", 300)) # will break
# the following will generate NA's for cells where the estimated
condition number of the covariance matrix exceeds 1e10:
pr1 = krige(zinc~1,meuse,meuse.grid,vgm(1, "Exp",
300),nmax=30,                 set=list(cn_max=1e10))
# generate a full interpolation as comparison:
pr1$idw = idw(zinc~1,meuse,meuse.grid)[[1]]
# show side by side:
spplot(pr1, c("var1.pred", "idw"), col.regions=bpy.colors())

... missing values are generated for those cells that result in a
singular covariance matrix, given their local neighbourhood of 30
nearest observations.

cn_max referers to the maximum allowed condition number, see
http://en.wikipedia.org/wiki/Condition_number

Two issues are (i) what singularity means when we use floating point
representations for real numbers, and (ii) that condition numbers of
matrices are expensive to evaluate, and an estimate based on the LU
decomposition is used. I threshold to 1e10 here, but that is purely for
illustrational purposes.

Note, for you of interest, that IIRC this thresholding is also done for
singularity of the X matrix, holding the predictors.

All this information didn't make it to the help pages of the krige
function. This help page would cover tens of pages otherwise. For full
documentation, still the "old", gstat stand-alone manual on gstat.org is
needed. I agree that this is not optimal.

With best wishes,
--
Edzer


Tomislav Hengl wrote:
> Hi Paul, Edzer,
>
> I understand why the singular matrix problem happens and I know that there is not really a
> mathematical solution around it:
>
>   
>> x <- matrix(runif(100), nrow=10)
>> x.i <- solve(x)
>> str(x.i)
>>     
>  num [1:10, 1:10] 0.8191 -1.0293 0.0826 1.068 -0.2106 ...
> # add a 'singular' column
>   
>> x[,1:10] <- rep(1, 10)
>> x.i <- solve(x)
>>     
> Error in solve.default(x) : 
>   Lapack routine dgesv: system is exactly singular
>
>
> However, I would very much support if you would integrate an "if" loop to check if it will happen,
> and then mask the prediction location using "NA".
>
> I mean, at the moment every time we want to run predictions, even if only at a single prediction
> location the model fails, we are not able to generate any output (this is sometimes really
> frustrating).
>
> Hence, I would support that you allow us to run predictions for all locations first, and then let
> the users judge if one should increase the search radius, remove some too-uniform predictors etc
> etc.
>
> I hope you agree with me,
>
> T. Hengl
> http://home.medewerker.uva.nl/t.hengl/
>
>
>   
>> -----Original Message-----
>> From: r-sig-geo-bounces at stat.math.ethz.ch [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf
>> Of Paul Hiemstra
>> Sent: Tuesday, November 17, 2009 9:53 AM
>> To: Edzer Pebesma
>> Cc: r-sig-geo
>> Subject: Re: [R-sig-Geo] 'LDLfactor' error in 'krige' function
>>
>> Edzer Pebesma wrote:
>>     
>>> My guess is that from constant data, the (co)variance is constant and
>>> zero, so the covariance matrix cannot be decomposed (hence: LDLfactor
>>> errors).
>>>
>>> Is this a case that autokrige should catch?
>>>
>>>       
>> I added a check in autoKrige. The output for the example below is now:
>>
>> kriging_result = autoKrige(zinc~1, meuse)
>> Error in autoKrige(zinc ~ 1, meuse) :
>>   All data in attribute 'zinc' is identical and equal to 0
>>    Can not interpolate this data
>>
>> A new version of automap with this feature has been uploaded to CRAN.
>>
>> cheers,
>> Paul
>>     
>>> --
>>> Edzer
>>>
>>> Mauricio Zambrano wrote:
>>>
>>>       
>>>> Dear List,
>>>>
>>>> During some OK interpolations of daily precipitation, with the
>>>> 'automap' library, I got the following error:
>>>>
>>>>
>>>> [using ordinary kriging]
>>>> "chfactor.c", line 130: singular matrix in function LDLfactor()
>>>> Error en predict.gstat(g, newdata = newdata, block = block, nsim = nsim,  :
>>>>   LDLfactor
>>>>
>>>>
>>>> The code I'm using works fine for other days, so I assume that the
>>>> distance between the measurement points is no the problem. When
>>>> looking at the data that rose the error, I realized that all the
>>>> measured values were equal to zero (I can not skip those days in which
>>>> all the measured points have the same value in advance, because the
>>>> measured value in those points change with time).
>>>>
>>>> According to a traceback that is given below, it seems that the cause
>>>> is in the 'predict.gstat' function of the 'gstat' package.
>>>>
>>>> The same error can be risen with:
>>>>
>>>> # Data preparation
>>>> data(meuse)
>>>> coordinates(meuse) =~ x+y
>>>> data(meuse.grid)
>>>> gridded(meuse.grid) =~ x+y
>>>>
>>>> meuse$zinc <- meuse$zinc*0
>>>> meuse$zinc
>>>>
>>>> # Ordinary kriging, no new_data object
>>>> kriging_result = autoKrige(zinc~1, meuse)
>>>>
>>>>
>>>> traceback()
>>>> 6: .Call("gstat_predict", as.integer(nrow(as.matrix(new.X))),
>>>> as.double(as.vector(raw$locations)),
>>>>        as.vector(new.X), as.integer(block.cols), as.vector(block),
>>>>        as.vector(bl_weights), as.integer(nsim), as.integer(BLUE))
>>>> 5: predict.gstat(g, newdata = newdata, block = block, nsim = nsim,
>>>>        indicators = indicators, na.action = na.action, debug.level =
>>>> debug.level)
>>>> 4: .local(formula, locations, ...)
>>>> 3: krige(formula, input_data, new_data, variogram_object$var_model,
>>>>        block = block, ...)
>>>> 2: krige(formula, input_data, new_data, variogram_object$var_model,
>>>>        block = block, ...)
>>>> 1: autoKrige(zinc ~ 1, meuse)
>>>>
>>>>
>>>> I would really appreciate any hint about the possible reason of this
>>>> error and if it could be overcome in some way.
>>>>
>>>>
>>>> Thanks in advance for any help.
>>>>
>>>>
>>>> Mauricio
>>>>
>>>>
>>>>         
>>>       
>> --
>> Drs. Paul Hiemstra
>> Department of Physical Geography
>> Faculty of Geosciences
>> University of Utrecht
>> Heidelberglaan 2
>> P.O. Box 80.115
>> 3508 TC Utrecht
>> Phone:  +3130 274 3113 Mon-Tue
>> Phone:  +3130 253 5773 Wed-Fri
>> http://intamap.geo.uu.nl/~paul
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>     
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>   

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de

More information about the R-sig-Geo mailing list