Find gradient of a posterior

Specific Domains Statistics
1

Please can anyone help me with this.

I have been trying to use ForwardDiff or Flux.Tracker or Zygote to get the gradient of this posterior but has not had any success yet:

using AltDistributions, Distributions 
using LinearAlgebra, PositiveFactorizations


### number of off-diagonals of covmat
function getnumElem(n)
        return(Int64.(n * (n - 1) / 2))
end


function model(V::Vector)      ### figure out dimensions to put all params in V

    parnum = getnumElem(d); parnum = Int64.(parnum); udo = parnum + d;
    par2 = 2*udo;
    bn = size(X,2)*d;


    M = zeros(d,d)
    M[tril!(trues(size(M)), -1)] .= V[1:parnum]
    M[diagind(M)] .= 1.
    sigmae = exp.(V[(parnum+1):udo])
    Ve = Diagonal(sigmae) * M * Diagonal(sigmae)
 

    Mu = zeros(d,d)
    Mu[tril!(trues(size(Mu)), -1)] .= V[(udo+1):udo+parnum]
    Mu[diagind(Mu)] .= 1.
    sigmau = exp.(V[(udo+parnum+1):par2])
    Vg = Diagonal(sigmau) * Mu * Diagonal(sigmau)
    Vgchol = cholesky(Positive, Vg).L

    beta = V[(par2 + 1):(par2 + bn)];
    beta = reshape(beta, :, d)

    In = one(rand(Float64, size(Achol)))  #### change
    gcov = kron(Vg, Achol) 
    rcov = kron(Ve, In)

    mu = vec(X*beta);

    Rcov = cholesky(Positive, gcov+rcov).L

    a = AltMvNormal(Val(:L), mu, Rcov)
    llk = logpdf(a, vec(Y))

    loglik = llk
    logpri = sum(logpdf.(Cauchy(2.5), sigmau)) + sum(logpdf.(Cauchy(2.5), sigmae)) + 
    logpdf(LKJL(d), LowerTriangular(Mu)) + logpdf(LKJL(d), LowerTriangular(M)) + loglikelihood(Normal(0.,1000.), beta)

    logprob = loglik + logpri

    return(-logprob)
 
end

Thanks.

Uche

2

What have you tried and how did it fail?

3

I have tried ForwardDiff.gradient

julia> ForwardDiff.gradient(vars -> model(vars[1:nvar]), init)
ERROR: TypeError: in setindex!, in typeassert, expected Float64, got ForwardDiff.Dual{Nothing,Float64,12}
Stacktrace:
 [1] setindex!(::Array{Float64,1}, ::ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##19#20")),Float64},Float64,12}, ::Int64) at ./array.jl:769
 [2] materialize!(::SubArray{Float64,1,Array{Float64,1},Tuple{Array{Int64,1}},false}, ::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(identity),Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##19#20")),Float64},Float64,12},1}}}) at ./subarray.jl:243
 [3] MT_logprob(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##19#20")),Float64},Float64,12},1}) at ./REPL[11]:9
 [4] (::getfield(Main, Symbol("##19#20")))(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##19#20")),Float64},Float64,12},1}) at ./REPL[47]:1
 [5] chunk_mode_gradient(::getfield(Main, Symbol("##19#20")), ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(Main, Symbol("##19#20")),Float64},Float64,12,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##19#20")),Float64},Float64,12},1}}) at /home/ubuntu/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:140
 [6] gradient(::Function, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(Main, Symbol("##19#20")),Float64},Float64,12,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##19#20")),Float64},Float64,12},1}}, ::Val{true}) at /home/ubuntu/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:19
 [7] gradient(::Function, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(Main, Symbol("##19#20")),Float64},Float64,12,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##19#20")),Float64},Float64,12},1}}) at /home/ubuntu/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:15 (repeats 2 times)
 [8] top-level scope at none:0

and Flux.Tracker with same error.

I don’t think that is the right way to do it.

I did Zygote and it crashed instantly.

4

Also using optimize and autodiff works

using Optim
opt = optimize(vars -> model(vars[1:nvar]), init; autodiff = :forward)

I relied on this before but I am trying to use Flux optimizers now.

5

Notice that the error says it failed to convert a Float64 to a ForwardDiff.Dual.
ForwardDiff users these dual numbers to differentiate, while Flux uses a tracker type.
In both cases, trying to assign one of these into an Array{Float64} will cause that error.
Now, you keep creating intermediate arrays like this:
M = zeros(d,d)
If you instead did this:
M = zeros(eltype(V),d,d)
That should help. I haven’t looked at your code thoroughly, so there may be other instances, but I think making those and similar changes should at least get ForwardDiff working.

3 Likes
6

Thanks Elrod. I tried it but it didn’t work.

function model(V::Vector)      ### figure out dimensions to put all params in V

           parnum = getnumElem(d); parnum = Int64.(parnum); udo = parnum + d;
           par2 = 2*udo;
           bn = size(X,2)*d;


           M = zeros(eltype(V), d,d)
           M[tril!(trues(size(M)), -1)] .= V[1:parnum]
           M[diagind(M)] .= 1.
           sigmae = exp.(V[(parnum+1):udo])
           Ve = Diagonal(sigmae) * M * Diagonal(sigmae)
        

           Mu = zeros(eltype(V), d,d)
           Mu[tril!(trues(size(Mu)), -1)] .= V[(udo+1):udo+parnum]
           Mu[diagind(Mu)] .= 1.
           sigmau = exp.(V[(udo+parnum+1):par2])
           Vg = Diagonal(sigmau) * Mu * Diagonal(sigmau)
           Vgchol = cholesky(Positive, Vg).L

           beta = V[(par2 + 1):(par2 + bn)];
           beta = reshape(beta, :, d)

           In = one(rand(Float64, size(Achol)))  #### change
           gcov = kron(Vg, Achol) 
           rcov = kron(Ve, In)

           mu = vec(X*beta);

           Rcov = cholesky(Positive, gcov+rcov).L

           a = AltMvNormal(Val(:L), mu, Rcov)
           llk = logpdf(a, vec(Y))

           loglik = llk
           logpri = sum(logpdf.(Cauchy(2.5), sigmau)) + sum(logpdf.(Cauchy(2.5), sigmae)) + 
           logpdf(LKJL(d), LowerTriangular(Mu)) + logpdf(LKJL(d), LowerTriangular(M)) + loglikelihood(Normal(0.,1000.), beta)

           logprob = loglik + logpri

           return(-logprob)
        
       end




julia> using ForwardDiff

julia> ForwardDiff.gradient(vars -> model(vars[1:nvar]), init)
ERROR: MethodError: no method matching floattype(::Type{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12}})
Closest candidates are:
  floattype(::Type{T<:AbstractFloat}) where T<:AbstractFloat at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:264
  floattype(::Type{T<:Integer}) where T<:Integer at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:265
Stacktrace:
 [1] default_δ(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},2}) at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:269
 [2] default_tol(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},2}) at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:270
 [3] cholesky(::Type{Positive}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},2}, ::Type) at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:7 (repeats 2 times)
 [4] model(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},1}) at ./REPL[11]:20
 [5] (::getfield(Main, Symbol("##3#4")))(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},1}) at ./REPL[18]:1
 [6] chunk_mode_gradient(::getfield(Main, Symbol("##3#4")), ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},1}}) at /home/ubuntu/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:140
 [7] gradient(::Function, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},1}}, ::Val{true}) at /home/ubuntu/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:19
 [8] gradient(::Function, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},1}}) at /home/ubuntu/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:15 (repeats 2 times)
 [9] top-level scope at none:0
7

Well, it did work in the sense that it fixed your setindex! error you reported earlier.
Now you have:

ERROR: MethodError: no method matching floattype(::Type{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12}})
Closest candidates are:
  floattype(::Type{T<:AbstractFloat}) where T<:AbstractFloat at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:264
  floattype(::Type{T<:Integer}) where T<:Integer at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:265
Stacktrace:
 [1] default_δ(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},2}) at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:269
 [2] default_tol(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},2}) at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:270
 [3] cholesky(::Type{Positive}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},2}, ::Type) at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:7 (repeats 2 times)

Try replacing cholesky(Positive, gcov+rcov).L with cholesky(Symmetric(S,:L)).L.

8

I use cholesky(Positive, S) because I need a pivoted cholesky. The LinearAlgebra cholesky fails some times for me.

9 
julia> using ForwardDiff, LinearAlgebra, PositiveFactorizations

julia> function foo(S)
           sum(cholesky(Symmetric(S,:L)).L)
       end
foo (generic function with 1 method)

julia> function posfact(S)
           sum(cholesky(Positive, S).L)
       end
posfact (generic function with 1 method)

julia> S = randn(20,30) |> x -> x * x';

julia> foo(S)
66.65333177508367

julia> posfact(S)
66.65333177508364

julia> ForwardDiff.gradient(foo, S)
20×20 Array{Float64,2}:
 0.19737    0.0         0.0       0.0       0.0         0.0          0.0        …  0.0       0.0       0.0        0.0       0.0       0.0       0.0       0.0     
 0.347832   0.2902      0.0       0.0       0.0         0.0          0.0           0.0       0.0       0.0        0.0       0.0       0.0       0.0       0.0     
 0.182511  -0.38832     0.424377  0.0       0.0         0.0          0.0           0.0       0.0       0.0        0.0       0.0       0.0       0.0       0.0     
 0.164754   0.346277   -0.357542  0.516962  0.0         0.0          0.0           0.0       0.0       0.0        0.0       0.0       0.0       0.0       0.0     
 0.289372   0.565278   -0.485473  0.761614  0.333315    0.0          0.0           0.0       0.0       0.0        0.0       0.0       0.0       0.0       0.0     
 0.175999  -0.079724    0.497135  0.374754  0.0220011   0.176914     0.0        …  0.0       0.0       0.0        0.0       0.0       0.0       0.0       0.0     
 0.195705   0.178183    0.111643  0.506256  0.26643     0.196079     0.119574      0.0       0.0       0.0        0.0       0.0       0.0       0.0       0.0     
 0.211179   0.405277   -0.26541   0.632575  0.457406    0.0378999    0.168757      0.0       0.0       0.0        0.0       0.0       0.0       0.0       0.0     
 0.13156   -0.0814023   0.345046  0.398375  0.0363314   0.381171     0.235919      0.0       0.0       0.0        0.0       0.0       0.0       0.0       0.0     
 0.34357    1.02695    -0.906925  0.760491  0.975419   -0.471576     0.108528      0.0       0.0       0.0        0.0       0.0       0.0       0.0       0.0     
 0.268457   0.589819   -0.310423  0.618467  0.554997   -0.134947     0.230623   …  0.0       0.0       0.0        0.0       0.0       0.0       0.0       0.0     
 0.303612   0.672425   -0.523016  0.598369  0.656003   -0.218669     0.119476      0.0       0.0       0.0        0.0       0.0       0.0       0.0       0.0     
 0.206404   0.452179   -0.213944  0.570041  0.446251    0.0416389    0.277697      0.140257  0.0       0.0        0.0       0.0       0.0       0.0       0.0     
 0.472759   0.900328   -0.719251  0.682745  0.91941    -0.533732    -0.0132312     0.305508  0.57741   0.0        0.0       0.0       0.0       0.0       0.0     
 0.208625   0.427204   -0.216818  0.496521  0.362886    0.126911     0.202587      0.29715   0.566542  0.0906131  0.0       0.0       0.0       0.0       0.0     
 0.282203   0.504234   -0.287621  0.532574  0.480854   -0.00698481   0.148356   …  0.282713  0.678763  0.174499   0.141724  0.0       0.0       0.0       0.0     
 0.323329   0.534302   -0.340189  0.549029  0.561977   -0.0906971    0.130069      0.287799  0.754959  0.182195   0.296788  0.195021  0.0       0.0       0.0     
 0.343886   0.559048   -0.387695  0.583671  0.640261   -0.161942     0.129573      0.280975  0.769863  0.181305   0.301756  0.425594  0.266087  0.0       0.0     
 0.281491   0.519284   -0.368454  0.522619  0.508265   -0.0296958    0.133605      0.305698  0.828918  0.212788   0.300172  0.367075  0.450412  0.13574   0.0     
 0.278452   0.520246   -0.370156  0.526275  0.509413   -0.0287605    0.136445      0.304156  0.828044  0.212449   0.301309  0.363753  0.45114   0.267974  0.137388

julia> ForwardDiff.gradient(posfact, S)
ERROR: MethodError: no method matching floattype(::Type{ForwardDiff.Dual{ForwardDiff.Tag{typeof(posfact),Float64},Float64,12}})
Closest candidates are:
  floattype(::Type{T<:AbstractFloat}) where T<:AbstractFloat at /home/chriselrod/.julia/packages/PositiveFactorizations/vmuKm/src/cholesky.jl:264
  floattype(::Type{T<:Integer}) where T<:Integer at /home/chriselrod/.julia/packages/PositiveFactorizations/vmuKm/src/cholesky.jl:265
Stacktrace:
 [1] default_δ(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(posfact),Float64},Float64,12},2}) at /home/chriselrod/.julia/packages/PositiveFactorizations/vmuKm/src/cholesky.jl:269
 [2] default_tol(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(posfact),Float64},Float64,12},2}) at /home/chriselrod/.julia/packages/PositiveFactorizations/vmuKm/src/cholesky.jl:270
 [3] cholesky(::Type{Positive}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(posfact),Float64},Float64,12},2}, ::Type) at /home/chriselrod/.julia/packages/PositiveFactorizations/vmuKm/src/cholesky.jl:7 (repeats 2 times)
 [4] posfact(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(posfact),Float64},Float64,12},2}) at ./REPL[28]:2
 [5] chunk_mode_gradient(::typeof(posfact), ::Array{Float64,2}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(posfact),Float64},Float64,12,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(posfact),Float64},Float64,12},2}}) at /home/chriselrod/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:140
 [6] gradient(::Function, ::Array{Float64,2}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(posfact),Float64},Float64,12,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(posfact),Float64},Float64,12},2}}, ::Val{true}) at /home/chriselrod/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:19
 [7] gradient(::Function, ::Array{Float64,2}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(posfact),Float64},Float64,12,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(posfact),Float64},Float64,12},2}}) at /home/chriselrod/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:15 (repeats 2 times)
 [8] top-level scope at REPL[33]:1

What is the error ?
Both versions of cholesky allow for pivoting, and both default to not pivoting.

1 Like
10

So I changed this part of the code above and have same error:

     mu = vec(X*beta);
           S = gcov+rcov;

           #Rcov = cholesky(Positive, S).L
           Rcov = cholesky(Symmetric(S,:L)).L

julia> ForwardDiff.gradient(vars -> model(vars[1:nvar]), init)
ERROR: MethodError: no method matching floattype(::Type{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12}})
Closest candidates are:
  floattype(::Type{T<:AbstractFloat}) where T<:AbstractFloat at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:264
  floattype(::Type{T<:Integer}) where T<:Integer at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:265
Stacktrace:
 [1] default_δ(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},2}) at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:269
 [2] default_tol(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},2}) at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:270
 [3] cholesky(::Type{Positive}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},2}, ::Type) at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:7 (repeats 2 times)
 [4] model(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},1}) at ./REPL[10]:20
 [5] (::getfield(Main, Symbol("##3#4")))(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},1}) at ./REPL[19]:1
 [6] chunk_mode_gradient(::getfield(Main, Symbol("##3#4")), ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},1}}) at /home/ubuntu/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:140
 [7] gradient(::Function, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},1}}, ::Val{true}) at /home/ubuntu/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:19
 [8] gradient(::Function, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},1}}) at /home/ubuntu/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:15 (repeats 2 times)
 [9] top-level scope at none:0

Or is this not what you mean?

11 
julia> ForwardDiff.gradient(vars -> model(vars[1:nvar]), init)
ERROR: MethodError: no method matching floattype(::Type{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12}})
Closest candidates are:
  floattype(::Type{T<:AbstractFloat}) where T<:AbstractFloat at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:264
  floattype(::Type{T<:Integer}) where T<:Integer at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:265
Stacktrace:
 [1] default_δ(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},2}) at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:269
 [2] default_tol(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},2}) at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:270
 [3] cholesky(::Type{Positive}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##3#4")),Float64},Float64,12},2}, ::Type) at /home/ubuntu/.julia/packages/PositiveFactorizations/VybOm/src/cholesky.jl:7 (repeats 2 times)

Note the error was in cholesky(Positive,...), not in cholesky(Symmetric...).

1 Like
12

I see…

So cholesky(Positive,S) is the culprit.

I do not like using LinearAlgebra pivoted cholesky cos I don’t wanna choose the value for tol.

cholesky(Positive,S) does not fail with default.

I get this kind of errors when I use pivoted cholesky from LinearAlgebra

julia> model(init)
ERROR: RankDeficientException(1)
Stacktrace:
 [1] chkfullrank at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/cholesky.jl:498 [inlined]
 [2] #cholesky!#92 at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/cholesky.jl:195 [inlined]
 [3] #cholesky! at ./none:0 [inlined]
 [4] #cholesky#96 at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/cholesky.jl:296 [inlined]
 [5] (::getfield(LinearAlgebra, Symbol("#kw##cholesky")))(::NamedTuple{(:tol, :check),Tuple{Float64,Bool}}, ::typeof(cholesky), ::Symmetric{Float64,Array{Float64,2}}, ::Val{true}) at ./none:0
 [6] model(::Array{Float64,1}) at ./REPL[24]:33
 [7] top-level scope at none:0
13

Anyways, thanks Elrod.

I will figure out how to fix the cholesky problem.

Uche.