Error: LinearAlgebra.SingularException(4)

I am trying to fit a Pumas.FOCEI model. I get the error message: LinearAlgebra,SingleException(4). Any pointers for me to look for? The error does not show up when i use Pumas. FO method.

Hi Anita,

Could you please post the entire error message. Best if you run the command that you used on the REPL and then copy the error message from there and post it here

Vijay, here is the error message:

LinearAlgebra.SingularException(4)
_A_ldiv_B at triangular.jl:404 [inlined]
\ at triangular.jl:25 [inlined]
macro expansion at solve.jl:56 [inlined]
_solve at solve.jl:46 [inlined]
\ at solve.jl:1 [inlined]
_analytical_solve(::Central1Periph1MetaPeriph1, ::Float64, ::Float64, ::LabelledArrays.SLArray{Tuple{4},ForwardDiff.Dual{ForwardDiff.Tag{getfield(Pumas, Symbol("##190#191")){Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}},PumasModel{ParamSet{NamedTuple{(:tvcl, :tvQ, :tvv, :tvvt, :tvt, :tvcl_metab, :tvQ_metab, :tvv_metab, :tvvt_metab, :Ω, :σ²_add, :σ²_prop),Tuple{RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},PDiagDomain{PDMats.PDiagMat{Float64,Array{Float64,1}}},VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}},VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}}}}},getfield(Main, Symbol("##49#54")),getfield(Main, Symbol("##50#55")),getfield(Main, Symbol("##51#56")),Central1Periph1MetaPeriph1,getfield(Main, Symbol("##52#57")),getfield(Main, Symbol("##53#58"))},Subject{NamedTuple{(:Conc_B, :Conc_Glucu),Tuple{Array{Union{Missing, Float64},1},Array{Union{Missing, Float64},1}}},NamedTuple{(:Age, :BSA),Tuple{Float64,Float64}},Array{Pumas.Event{Float64,Float64,Float64,Float64,Float64,Float64,Int64},1},Array{Float64,1}},NamedTuple{(:tvcl, :tvQ, :tvv, :tvvt, :tvt, :tvcl_metab, :tvQ_metab, :tvv_metab, :tvvt_metab, :Ω, :σ²_add, :σ²_prop),Tuple{Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,PDMats.PDiagMat{Float64,Array{Float64,1}},Array{Float64,1},Array{Float64,1}}},Tuple{}},Float64},Float64,9},1,4,(:Central, :CPeripheral, :Metabolite, :MPeripheral)}, ::LabelledArrays.SLArray{Tuple{4},ForwardDiff.Dual{ForwardDiff.Tag{getfield(Pumas, Symbol("##190#191")){Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}},PumasModel{ParamSet{NamedTuple{(:tvcl, :tvQ, :tvv, :tvvt, :tvt, :tvcl_metab, :tvQ_metab, :tvv_metab, :tvvt_metab, :Ω, :σ²_add, :σ²_prop),Tuple{RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},PDiagDomain{PDMats.PDiagMat{Float64,Array{Float64,1}}},VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}},VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}}}}},getfield(Main, Symbol("##49#54")),getfield(Main, Symbol("##50#55")),getfield(Main, Symbol("##51#56")),Central1Periph1MetaPeriph1,getfield(Main, Symbol("##52#57")),getfield(Main, Symbol("##53#58"))},Subject{NamedTuple{(:Conc_B, :Conc_Glucu),Tuple{Array{Union{Missing, Float64},1},Array{Union{Missing, Float64},1}}},NamedTuple{(:Age, :BSA),Tuple{Float64,Float64}},Array{Pumas.Event{Float64,Float64,Float64,Float64,Float64,Float64,Int64},1},Array{Float64,1}},NamedTuple{(:tvcl, :tvQ, :tvv, :tvvt, :tvt, :tvcl_metab, :tvQ_metab, :tvv_metab, :tvvt_metab, :Ω, :σ²_add, :σ²_prop),Tuple{Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,PDMats.PDiagMat{Float64,Array{Float64,1}},Array{Float64,1},Array{Float64,1}}},Tuple{}},Float64},Float64,9},1,4,(:Central, :CPeripheral, :Metabolite, :MPeripheral)}, ::NamedTuple{(:CL1, :Q1, :V1, :Vp1, :T, :CL2, :Q2, :V2, :Vp2, :tvcl, :tvQ, :tvv, :tvvt, :tvt, :tvcl_metab, :tvQ_metab, :tvv_metab, :tvvt_metab, :Ω, :σ²_add, :σ²_prop, :η),Tuple{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Pumas, Symbol("##190#191")){Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}},PumasModel{ParamSet{NamedTuple{(:tvcl, :tvQ, :tvv, :tvvt, :tvt, :tvcl_metab, :tvQ_metab, :tvv_metab, :tvvt_metab, :Ω, :σ²_add, :σ²_prop),Tuple{RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Flo...

Below is the code. I am using FO method for now. FO method gives the output . But the infer function gives the error: PosDefException:matrix is not positive definite; Cholesky fractionation failed

parent_glu = @model begin
  @param   begin
     tvcl        ∈    RealDomain(lower=0)
     tvQ         ∈    RealDomain(lower=0)
     tvv         ∈    RealDomain(lower=0)
     tvvt        ∈    RealDomain(lower=0)
     tvt         ∈    RealDomain(lower=0)
     tvcl_metab  ∈    RealDomain(lower=0)
     tvQ_metab   ∈    RealDomain(lower=0)
     tvv_metab   ∈    RealDomain(lower=0)
     tvvt_metab  ∈    RealDomain(lower=0)
     Ω           ∈    PDiagDomain(9)
     σ²_add      ∈    VectorDomain(2, lower=[0.001, 0.001])
     σ²_prop     ∈    VectorDomain(2, lower=[0.001, 0.001])
  end

  @random begin
     η ~ MvNormal(Ω)
  end

  @pre begin
     CL1      =    tvcl  * exp(η[1]) # CLpar
     Q1       =    tvQ   * exp(η[2]) # Qpar
     V1       =    tvv   * exp(η[3]) # Cp_par
     Vp1      =    tvvt  * exp(η[4]) # Ct_par
     T        =    tvt  * exp(η[5]) # CLpm
     CL2      =    tvcl_metab * exp(η[6]) # CLglu
     Q2       =    tvQ_metab * exp(η[7])  # Qglu
     V2       =    tvv_metab * exp(η[8])  # Cp_glu
     Vp2      =    tvvt_metab * exp(η[9]) # Ct_glu
  end

  @dynamics Central1Periph1MetaPeriph1

  @derived begin
    Cparent     :=    @. Central / V1
    Conc_B      ~     @. Normal(Cparent, sqrt(Cparent^2*σ²_prop[1] + σ²_add[1]))
    CMetab      :=    @. Metabolite / V2
    Conc_Glucu   ~    @. Normal(CMetab, sqrt(CMetab^2*σ²_prop[2] + σ²_add[2]))
  end
end

param_normal = (
    tvcl       =     0.9,
    tvQ        =     8,
    tvv        =     6.6,
    tvvt       =     5,
    tvt        =     44,
    tvcl_metab =     4,
    tvQ_metab  =     5,
    tvv_metab  =     7.6,
    tvvt_metab =     6,
    Ω          =     Diagonal([0.09,0.09,0.09,0.09,0.09,0.09,0.09,0.09,0.09]),
    σ²_add     =     [0.01, 0.01],
    σ²_prop    =     [0.01, 0.01]
    )

result_parent_glu = @time fit(parent_glu,
                              df_belino_pk,
                              param_normal,
                              Pumas.FO()
                            )
infer_parent_glu = infer(result_parent_glu)
PosDefException: matrix is not positive definite; Cholesky factorization failed.
chkposdef at lapack.jl:50 [inlined]
sygvd!(::Int64, ::Char, ::Char, ::Array{Float64,2}, ::Array{Float64,2}) at lapack.jl:5075
#eigen!#85 at symmetric.jl:677 [inlined]
eigen! at symmetric.jl:677 [inlined]
#eigen#68(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}}, ::typeof(LinearAlgebra.eigen), ::LinearAlgebra.Symmetric{Float64,Array{Float64,2}}, ::LinearAlgebra.Symmetric{Float64,Array{Float64,2}}) at eigen.jl:403
eigen(::LinearAlgebra.Symmetric{Float64,Array{Float64,2}}, ::LinearAlgebra.Symmetric{Float64,Array{Float64,2}}) at eigen.jl:402
#vcov#253(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}}, ::typeof(vcov), ::Pumas.FittedPumasModel{PumasModel{ParamSet{NamedTuple{(:tvcl, :tvQ, :tvv, :tvvt, :tvt, :tvcl_metab, :tvQ_metab, :tvv_metab, :tvvt_metab, :Ω, :σ²_add, :σ²_prop),Tuple{RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},PDiagDomain{PDMats.PDiagMat{Float64,Array{Float64,1}}},VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}},VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}}}}},getfield(Main, Symbol("##99#104")),getfield(Main, Symbol("##100#105")),getfield(Main, Symbol("##101#106")),Central1Periph1MetaPeriph1,getfield(Main, Symbol("##102#107")),getfield(Main, Symbol("##103#108"))},Array{Subject{NamedTuple{(:Conc_B, :Conc_Glucu),Tuple{Array{Union{Missing, Float64},1},Array{Union{Missing, Float64},1}}},NamedTuple{(:Age, :BSA),Tuple{Float64,Float64}},Array{Pumas.Event{Float64,Float64,Float64,Float64,Float64,Float64,Int64},1},Array{Float64,1}},1},Optim.MultivariateOptimizationResults{Optim.BFGS{LineSearches.InitialStatic{Float64},LineSearches.BackTracking{Float64,Int64},getfield(Pumas, Symbol("##211#212")){NLSolversBase.OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}},Array{Float64,1}},Nothing,Optim.Flat},Float64,Array{Float64,1},Float64,Float64,Array{Optim.OptimizationState{Float64,Optim.BFGS{LineSearches.InitialStatic{Float64},LineSearches.BackTracking{Float64,Int64},getfield(Pumas, Symbol("##211#212")){NLSolversBase.OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}},Array{Float64,1}},Nothing,Optim.Flat}},1},Bool},Pumas.FO,Array{Array{Float64,1},1},Tuple{},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}},TransformVariables.TransformTuple{NamedTuple{(:tvcl, :tvQ, :tvv, :tvvt, :tvt, :tvcl_metab, :tvQ_metab, :tvv_metab, :tvvt_metab, :Ω, :σ²_add, :σ²_prop),Tuple{TransformVariables.ShiftedExp{true,Int64},TransformVariables.ShiftedExp{true,Int64},TransformVariables.ShiftedExp{true,Int64},TransformVariables.ShiftedExp{true,Int64},TransformVariables.ShiftedExp{true,Int64},TransformVariables.ShiftedExp{true,Int64},TransformVariables.ShiftedExp{true,Int64},TransformVariables.ShiftedExp{true,Int64},TransformVariables.ShiftedExp{true,Int64},Pumas.PDiagTransform,Pumas.ElementArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},Pumas.ElementArrayTransform{TransformVariables.ShiftedExp{true,Float64},1}}}}}) at likelihoods.jl:1364
vcov at likelihoods.jl:1361 [inlined]
#infer#272(::Float64, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}}, ::typeof(infer), ::Pumas.FittedPumasModel{PumasModel{ParamSet{NamedTuple{(:tvcl, :tvQ, :tvv, :tvvt, :tvt, :tvcl_metab, :tvQ_metab, :tvv_metab, :tvvt_metab, :Ω, :σ²_add, :σ²_prop),Tuple{RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},RealDomain{Int64,TransformVariables.Infinity{true},Float64},PDiagDomain{PDMats.PDiagMat{Float64,Array{Float64,1}}},VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}},VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}}}}},getfield(Main, Symbol("##99#104")),getfield(Main, Symbol("##100#105")),getfield(Main, Symbol("##101#106")),Central1Periph1MetaPeriph1,getfield(Main, Symbol("##102#107")),getfield(Main, Symbol("##103#108"))},Array{Subject{NamedTuple{(:Conc_B, :Conc_Glucu),Tuple{Array{Union{Mis...

this error comes when you try to fit the model using Pumas.FOCEI()?

The error message is for FO methhod at the infer step. (PosDefException: matrix is not positive definite; Cholesky factorization failed.)

what function triggered the error message for Error:LinearAlgebra.SingularException(4), the title of your post?

This was for Pumas. FOCEI ()

1 Like

These are two different errors. The error during the FOCEI fitting is because the eigenspace of the ODE coefficient matrix collapses completely. Likely because of some of the PK parameters either ending up at zero or infinity. Would you be able to add

optimize_fn=Pumas.DefaultOptimizeFN(show_trace=true)

as a last argument to the fit function. That will print the optimization trace. Then please report the values of the last step before it fails.

The error you get from the infer of the FO fit is because outer score estimate of the covariance matrix is singular. (Indeed we should provide a better error message). Again, this could happen if some of the parameters end up with extreme values. Would you be able to show the result of just the fit? Before seeing the estimates, a guess could be that this is related to the number of random effects in the model. Unless you have very rich data, the model is probably way overparameterized.

1 Like

Thanks for the technical insights. For FO method, my initial PK estimates were way off, so i had to tweak the initial estimates. For FOCEI, my error model estimates were way off, so i had to change the initial estimates. Now i am able to run both methods.

1 Like