using Surrogates, Stheno, Distributions
d = Normal()
dat = rand(d, 300)
function test(x)
x = collect(x)
μ = x[1]
logsigma = x[2]
σ = exp(logsigma)
d2 = Normal(μ, σ)
-log(exp(loglikelihood(d2, dat)))
end
#optimize(test, [1., 1.], BFGS(), autodiff = :forward) Works fine with Optim
ig = [1., 1.]
test2(x) = test(collect(x))
l = ig.*0.5
u = ig .* 2.
s = Surrogates.sample(20, l, u, SobolSample())
fs = test2.(s)
ind = isfinite.(fs)
σ² = 0.05
gp = Stheno.GP(σ² * stretch(matern52(), 1.), Stheno.GPC())
my_krig = SthenoKriging(s[ind], fs[ind], gp)
surrogate_optimize(test, LCBS(), l, u, my_krig, SobolSample())
I actually think there is a seperate issue here. When I run this example, I get
ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
[1] checkpositivedefinite at /home/andrew/Programs/julia-1.3.1-linux-x86_64/julia-1.3.1/share/julia/stdlib/v1.3/LinearAlgebra/src/factorization.jl:18 [inlined]
[2] #cholesky!#124(::Bool, ::typeof(LinearAlgebra.cholesky!), ::LinearAlgebra.Symmetric{Float64,Array{Float64,2}}, ::Val{false}) at /home/andrew/Programs/julia-1.3.1-linux-x86_64/julia-1.3.1/share/julia/stdlib/v1.3/LinearAlgebra/src/cholesky.jl:226
[3] #cholesky! at ./none:0 [inlined] (repeats 2 times)
[4] #cholesky#129 at /home/andrew/Programs/julia-1.3.1-linux-x86_64/julia-1.3.1/share/julia/stdlib/v1.3/LinearAlgebra/src/cholesky.jl:348 [inlined]
[5] cholesky at /home/andrew/Programs/julia-1.3.1-linux-x86_64/julia-1.3.1/share/julia/stdlib/v1.3/LinearAlgebra/src/cholesky.jl:348 [inlined] (repeats 2 times)
[6] |(::GP{Stheno.ZeroMean{Float64},Stheno.EQ}, ::Stheno.Observation{Stheno.FiniteGP{Stheno.CompositeGP{Tuple{typeof(LinearAlgebra.cross),Array{GP{Stheno.ZeroMean{Float64},Stheno.EQ},1}}},Stheno.BlockData{Array{Float64,1},ColVecs{Float64,Array{Float64,2}}},BlockArrays.BlockArray{Float64,2,LinearAlgebra.Diagonal{LinearAlgebra.Diagonal{Float64,Array{Float64,1}},Array{LinearAlgebra.Diagonal{Float64,Array{Float64,1}},1}},BlockArrays.BlockSizes{2,Tuple{Array{Int64,1},Array{Int64,1}}}}},Array{Float64,1}}) at /home/andrew/.julia/packages/Stheno/j89z0/src/composite/conditioning.jl:9
[7] map at /home/andrew/.julia/packages/Stheno/j89z0/src/composite/conditioning.jl:69 [inlined]
[8] |(::Tuple{GP{Stheno.ZeroMean{Float64},Stheno.EQ}}, ::Tuple{Stheno.Observation{Stheno.FiniteGP{GP{Stheno.ZeroMean{Float64},Stheno.EQ},ColVecs{Float64,Array{Float64,2}},LinearAlgebra.Diagonal{Float64,Array{Float64,1}}},Array{Float64,1}}}) at /home/andrew/.julia/packages/Stheno/j89z0/src/composite/conditioning.jl:71
[9] _condition_gps(::ColVecs{Float64,Array{Float64,2}}, ::Array{Array{Float64,1},1}, ::Tuple{GP{Stheno.ZeroMean{Float64},Stheno.EQ}}, ::Float64) at /home/andrew/.julia/packages/Surrogates/dfpEW/src/SthenoKriging.jl:95
[10] _prepare_gps at /home/andrew/.julia/packages/Surrogates/dfpEW/src/SthenoKriging.jl:82 [inlined]
[11] add_point!(::SthenoKriging{Array{Tuple{Float64,Float64},1},Array{Float64,1},Tuple{GP{Stheno.ZeroMean{Float64},Stheno.EQ}},Float64,Tuple{Stheno.CompositeGP{Tuple{typeof(|),GP{Stheno.ZeroMean{Float64},Stheno.EQ},LinearAlgebra.Cholesky{Float64,Array{Float64,2}},Array{Float64,1},Stheno.FiniteGP{Stheno.CompositeGP{Tuple{typeof(LinearAlgebra.cross),Array{GP{Stheno.ZeroMean{Float64},Stheno.EQ},1}}},Stheno.BlockData{Array{Float64,1},ColVecs{Float64,Array{Float64,2}}},BlockArrays.BlockArray{Float64,2,LinearAlgebra.Diagonal{LinearAlgebra.Diagonal{Float64,Array{Float64,1}},Array{LinearAlgebra.Diagonal{Float64,Array{Float64,1}},1}},BlockArrays.BlockSizes{2,Tuple{Array{Int64,1},Array{Int64,1}}}}},Array{Float64,1}}}}}, ::Tuple{Float64,Float64}, ::Float64) at /home/andrew/.julia/packages/Surrogates/dfpEW/src/SthenoKriging.jl:75
I get the same error even before the optimization step if I construct the Sobol grid with 100 points instead of 20. I'm not really sure why. It doesn't seem to have anything to do with the objective. I got the same thing when I used random function values.
Anyway, regardless, this objective function will be infinite valued at some points due to poor numerical conditioning. What should Surrogates do in that case? Decline to add that value and try something else?
On my actual problem I didn't get this positive definite issue. Instead, it ran the optimization step for about 10 points, then said