R Solve:System Is Exactly Singular

Error in solve.default(oout$hessian) : Lapack routine dgesv: system is exactly singular: U[1,1] = 0

tl;dr I'm not sure your objective function makes sense, my guess is that you have a typo. (Furthermore, if your objective function is working with mle, you don't need to set fixed explicitly: the profile method will automatically compute likelihood profiles for you ...)

Let's start with the full model, and let's use optim() rather than stats4::mle() (I know you want to get back to mle so that you can do likelihood profiling, but it's a little bit easier to debug optim() problems since there is one less layer of code to dig through.)

Since optim() wants an objective function that takes a vector rather than a list of arguments, write a wrapper (we could also use do.call(fr, as.list(p))):

fr0 <- function(p) {
  fr(x1=p[1], x2=p[2], x3=p[3])
}
opt1 <- optim(fn=fr0, par=c(1,2,3), method="Nelder-Mead")

Results:

$par
[1]    28.51486   812.09978 -7095.39630

$value
[1] -6238.881

$counts
function gradient 
     502       NA 

$convergence
[1] 1

Note that the value of x[3] is strongly negative, as is the objective function value, and the convergence code is non-zero: in particular means (from ?optim):

‘1’ indicates that the iteration limit ‘maxit’ had been reached.

If we set control=list(maxit=2000) and try again x3 and the objective function get even smaller and the convergence code is still 1!

Then we look more carefully at the objective function and notice that it goes to -Inf as x3 → Inf, so we'll never reach the answer. (Presumably at some point we'll get to a floating-point problem, but 10 million iterations only gets us as a far as -1e17 ...)

If I change x3 to x3^2 in your function everything seems to work OK ... maybe that's what you intended ... ???

In R, kriging gives system is exactly singular error

The problem is, as @mrip mentioned, that you have duplicate observations. In this context, points with exactly the same x and y coordinates. The following code reveals the duplicates:

coor = cbind(x, y)
duplicated(coor)
# [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#[13] FALSE  TRUE  TRUE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE FALSE
#[25]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE
coor[duplicated(coor),]
#     [,1] [,2]
#[1,]   45   40
#[2,]   25   60
#[3,]   25   55
#[4,]   40   35
#[5,]   25   50
#[6,]   55   25
#[7,]   40   40
#[8,]   25   55
#[9,]   55   20
#[10,]   40   35
#[11,]   25   50
#[12,]   50   20
#[13,]   25   45

This leads to an uninvertible covariance matrix in the kriging equations, which in turn leads to the Lapack error from dgesv. The solution is to remove the duplicates from the dataset.

The following dataset works:

x = runif(100)
y = runif(100)
z = runif(100)
krig = kriging(x, y, z, lag = 3)
image(krig)

Related Topics