LinearAlgebra.LAPACKException(3) ERROR when calling `confidence_intervals`

[Datseris]

Hi,

I'm running code that guaranteed worked just a week ago, but since then I've done some package updates, and also moved to Julia 1.6. Julia 1.6 is not the problem, same error I get when falling back to Julia 1.5.1.

Here is a MWE:

using LsqFit
xs, ys = rand(100), rand(100)
p0 = [1.0, 1.0, 0.01]
regression(x, p) = @. p[1] + p[2]*x #+ p[3]*log(-x)
fit = LsqFit.curve_fit(regression, xs, ys, p0)
dfit = LsqFit.coef(fit)[2]
dci = LsqFit.confidence_interval(fit, 0.05)[2]

ERROR: LoadError: LinearAlgebra.LAPACKException(3)
Stacktrace:
  [1] chklapackerror(ret::Int64)
    @ LinearAlgebra.LAPACK C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\lapack.jl:38
  [2] trtrs!(uplo::Char, trans::Char, diag::Char, A::Matrix{Float64}, B::Matrix{Float64}) 
    @ LinearAlgebra.LAPACK C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\lapack.jl:3426
  [3] ldiv!
    @ C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\triangular.jl:739 [inlined]
  [4] inv(A::LinearAlgebra.UpperTriangular{Float64, Matrix{Float64}})
    @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\triangular.jl:821
  [5] inv(A::Matrix{Float64})
    @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\dense.jl:811
  [6] estimate_covar(fit::LsqFit.LsqFitResult{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Vector{Float64}})
    @ LsqFit ~\.julia\packages\LsqFit\wmSYK\src\curve_fit.jl:203
  [7] stderror(fit::LsqFit.LsqFitResult{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Vector{Float64}}; rtol::Float64, atol::Int64)
    @ LsqFit ~\.julia\packages\LsqFit\wmSYK\src\curve_fit.jl:217
  [8] confidence_interval(fit::LsqFit.LsqFitResult{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, alpha::Float64; rtol::Float64, atol::Int64)
    @ LsqFit ~\.julia\packages\LsqFit\wmSYK\src\curve_fit.jl:246
  [9] confidence_interval(fit::LsqFit.LsqFitResult{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Vector{Float64}}, alpha::Float64)
    @ LsqFit ~\.julia\packages\LsqFit\wmSYK\src\curve_fit.jl:246

Any ideas how to fix?

p.s.: While making a MWE, I've checked the Tests folder of this package, and couldn't find a call to confidence_interval anywhere... Where is this function tested?

[Datseris]

version: [2fda8390] LsqFit v0.12.0.

[pkofod]

I would guess singular or near-singular information matrix (or whatever equivalent is used here).

[pkofod]

When will I ever learn to read the provided example before answering. The problem is exactly that. Your third parameter gives you a zero in the hessian of the likelihood function (the interpretation we use when we think about standard errors). Your third gradient entry is exactly, constantly zero so your last row and column of your hessian is too.

G(p) = [g1, g2, 0]

H(p) = [h1  h2  0
        h2  h3  0
        0   0   0]

something like that. In other words, the third parameter is unidentified. It could be anything. https://en.wikipedia.org/wiki/Identifiability

[Datseris]

Thank you. Although a bit too late for when I was using this, nevertheless I appreciate the answer. So to conclude, I shouldn't use "null" parameters and always use as many parameters as I actually use in my model.

A suggestion for LsqFit: to throw a warning when this happens? If I read your message correctly, the package has all the information it needs to make such a statement?

[pkofod]

I knew it from your model. I guess it can be a bit hard to check. But if the inv (should be a factorization, but...) fails, we could certainly catch it and ask the user to check for identification, yes.