JuliaStats / Klara.jl

MCMC inference in Julia
Other
166 stars 38 forks source link

make chol() happy on 0.5 #146

Closed omalled closed 7 years ago

omalled commented 7 years ago

Mads' tests fail ("import Mads; Mads.test()") inside Klara on julia 0.5 due to the changes in chol(). This PR fixes it, but I don't think this code will work on 0.4. I can add some "if VERSION < ... else ... end" type stuff, if that's the direction you guys would like to go.

papamarkou commented 7 years ago

Thanks a lot @omalled. The change from chol(C, Val{:L}) to ctranspose(chol(C)) should be cross-compatible between Julia v0.4 and v0.5 if I am not wrong, so it is the right way to fix this. The call to Hermitian() is not necessary though and could be dropped I think - what was your motivation behind calling Hermitian()?

omalled commented 7 years ago

Here's the error message I was getting that lead me to add Hermitian():

ERROR: LoadError: LoadError: ArgumentError: matrix is not symmetric/Hermitian. This error can be avoided by calling chol(Hermitian(A)) which will ignore either the upper or lower triangle of the matrix. in chol(::Array{Float64,2}) at ./linalg/cholesky.jl:166 in chol(::Array{Float64,2}, ::Type{Val{:L}}) at ./deprecated.jl:52 in ##_iterate#371(::Klara.BasicMCJob) at /homedir/.julia/v0.5/Klara/src/samplers/iterate/RAM.jl:107 in ##_run#372(::Klara.BasicMCJob) at /homedir/.julia/v0.5/Klara/src/jobs/BasicMCJob.jl:308 in #bayessampling#179(::Int64, ::Int64, ::Int64, ::Int64, ::Function, ::Dict{Any,Any}) at /homedir/.julia/v0.5/Mads/src/MadsMC.jl:48 in (::Mads.#kw##bayessampling)(::Array{Any,1}, ::Mads.#bayessampling, ::Dict{Any,Any}) at ./:0 in include_from_node1(::String) at ./loading.jl:426 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:? in macro expansion; at /homedir/.julia/v0.5/Mads/test/runtests.jl:17 [inlined] in anonymous at ./:? in include_from_node1(::String) at ./loading.jl:426 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:? in test(::String) at /homedir/.julia/v0.5/Mads/src/MadsTest.jl:24 in test() at /homedir/.julia/v0.5/Mads/src/MadsTest.jl:22 while loading /homedir/.julia/v0.5/Mads/examples/bayesian_sampling/runtests.jl, in expression starting on line 12 while loading /homedir/.julia/v0.5/Mads/test/runtests.jl, in expression starting on line 13

papamarkou commented 7 years ago

Thanks @omalled. There are two convoluted matters, which we may break into two standalone issues. One relates to fixing the chol(C, Val{:L}) breakage in v0.4, which you fixed already by using ctranspose(chol(C)). Would you like to keep this change in this PR, and remove for now the call to Hermitian() so that we investigate the latter in more detail?

If you agree, once you make the above change to this PR, I will merge it and will then open a new issue re the hermitian matter. A Cholesky decomposition usually fails because the matrix to be decomposed is not positive-definite. However, Vihola's RAM sampler constructs a positive-definite matrix, so it would be good to try and get to the bottom of this, i.e. to see why you get this error. Do you have a relevant example-script (or alternatively can you isolate the error via a simple example), which we can then further examine in a new issue, after merging the relevant aspects of the current PR?

omalled commented 7 years ago

I removed Hermitian() in the PR. Some pretty minimal code that reproduces the error we are seeing in Mads is below. Would you like me to create a separate issue for that?

import Klara
logtarget(x) = -sum(x .^ 2)
mcparams = Klara.BasicContMuvParameter(:p, logtarget=logtarget)
model = Klara.likelihood_model(mcparams, false)
initvals = zeros(3)
sampler = Klara.RAM(fill(1e-1, length(initvals)))
mcrange = Klara.BasicMCRange(nsteps=1000, burnin=100, thinning=1)
mcparams0 = Dict(:p=>initvals)
job = Klara.BasicMCJob(model, sampler, mcrange, mcparams0, tuner=Klara.VanillaMCTuner())
Klara.run(job)
papamarkou commented 7 years ago

Thanks @omalled, I merged the PR. Sure, you may open a new issue about the Hermitian matter.