Optimization with SRBF results in SingularException

[christianhauschel]

Describe the bug 🐞

Optimization with SRBF (linear, quadratic, and cubic Radial) results in SingularException, optimization with EI + Kriging does not.

Minimal Reproducible Example 👇

using Surrogates 

x0 = [0.5, 0.2]
lb = 0.5 * x0
ub = 1.5 * x0

function obj(x)
    return x[1]^2 + x[2]^2
end

x = sample(5, lb, ub, SobolSample())

y = [obj(x[i]) for i = 1:length(x)]

surrogate_init = RadialBasis(x, y, lb, ub; rad=linearRadial())

x_opt, y_opt = surrogate_optimize(
    obj,
    SRBF(),
    lb,
    ub,
    surrogate_init,
    SobolSample();
)

# surrogate_init = Kriging(x, y, lb, ub)

# x_opt, y_opt = surrogate_optimize(
#     obj,
#     EI(),
#     lb,
#     ub,
#     surrogate_init,
#     SobolSample();
# )

Error & Stacktrace ⚠️

ERROR: LinearAlgebra.SingularException(22)
Stacktrace:
  [1] checknonsingular
    @ ~/.julia/juliaup/julia-1.8.5+0.x64.linux.gnu/share/julia/stdlib/v1.8/LinearAlgebra/src/factorization.jl:19 [inlined]
  [2] checknonsingular
    @ ~/.julia/juliaup/julia-1.8.5+0.x64.linux.gnu/share/julia/stdlib/v1.8/LinearAlgebra/src/factorization.jl:21 [inlined]
  [3] #bunchkaufman!#182
    @ ~/.julia/juliaup/julia-1.8.5+0.x64.linux.gnu/share/julia/stdlib/v1.8/LinearAlgebra/src/bunchkaufman.jl:102 [inlined]
  [4] #bunchkaufman#185
    @ ~/.julia/juliaup/julia-1.8.5+0.x64.linux.gnu/share/julia/stdlib/v1.8/LinearAlgebra/src/bunchkaufman.jl:199 [inlined]
  [5] #_factorize#124
    @ ~/.julia/juliaup/julia-1.8.5+0.x64.linux.gnu/share/julia/stdlib/v1.8/LinearAlgebra/src/symmetric.jl:615 [inlined]
  [6] _factorize
    @ ~/.julia/juliaup/julia-1.8.5+0.x64.linux.gnu/share/julia/stdlib/v1.8/LinearAlgebra/src/symmetric.jl:612 [inlined]
  [7] factorize
    @ ~/.julia/juliaup/julia-1.8.5+0.x64.linux.gnu/share/julia/stdlib/v1.8/LinearAlgebra/src/symmetric.jl:611 [inlined]
  [8] \
    @ ~/.julia/juliaup/julia-1.8.5+0.x64.linux.gnu/share/julia/stdlib/v1.8/LinearAlgebra/src/symmetric.jl:625 [inlined]
  [9] _calc_coeffs(x::Vector{Tuple{Float64, Float64}}, y::Vector{Float64}, lb::Vector{Float64}, ub::Vector{Float64}, phi::Function, q::Int64, scale_factor::Float64, sparse::Bool)
    @ Surrogates ~/.julia/packages/Surrogates/mN4Fb/src/Radials.jl:62
 [10] add_point!(rad::RadialBasis{Surrogates.var"#1#2", Int64, Vector{Tuple{Float64, Float64}}, Vector{Float64}, Vector{Float64}, Vector{Float64}, LinearAlgebra.Transpose{Float64, Vector{Float64}}, Float64, Bool}, new_x::Tuple{Float64, Float64}, new_y::Float64)
    @ Surrogates ~/.julia/packages/Surrogates/mN4Fb/src/Radials.jl:232
 [11] surrogate_optimize(obj::typeof(obj), ::SRBF, lb::Vector{Float64}, ub::Vector{Float64}, surr::RadialBasis{Surrogates.var"#1#2", Int64, Vector{Tuple{Float64, Float64}}, Vector{Float64}, Vector{Float64}, Vector{Float64}, LinearAlgebra.Transpose{Float64, Vector{Float64}}, Float64, Bool}, sample_type::SobolSample; maxiters::Int64, num_new_samples::Int64)
    @ Surrogates ~/.julia/packages/Surrogates/mN4Fb/src/Optimization.jl:176
 [12] surrogate_optimize(obj::Function, ::SRBF, lb::Vector{Float64}, ub::Vector{Float64}, surr::RadialBasis{Surrogates.var"#1#2", Int64, Vector{Tuple{Float64, Float64}}, Vector{Float64}, Vector{Float64}, Vector{Float64}, LinearAlgebra.Transpose{Float64, Vector{Float64}}, Float64, Bool}, sample_type::SobolSample)
    @ Surrogates ~/.julia/packages/Surrogates/mN4Fb/src/Optimization.jl:77
 [13] top-level scope
    @ ~/projects/conoise/CustomMDAO/test.jl:17

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

  [99985d1d] AbstractGPs v0.5.21
  [f67ccb44] HDF5 v0.17.2
  [91a5bcdd] Plots v1.40.4
  [f789e72d] PrettySections v0.1.0 `https://github.com/christianhauschel/PrettySections.jl.git#master`
  [5f89f4a4] PyFormattedStrings v0.1.11
⌅ [6fc51010] Surrogates v6.0.0
  [78aa1720] SurrogatesAbstractGPs v0.1.0
  [10745b16] Statistics
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

• Output of versioninfo()

Julia Version 1.8.5

[christianhauschel]

Same for:

  [b964fa9f] LaTeXStrings v1.3.1
  [91a5bcdd] Plots v1.40.4
  [8a4e6c94] QuasiMonteCarlo v0.3.3
  [6fc51010] Surrogates v6.9.0
  [50679fc6] SurrogatesPolyChaos v0.1.0

Julia Version 1.10.2

MVE

using Surrogates 
using QuasiMonteCarlo

x0 = [0.5, 0.2]
lb = 0.5 * x0
ub = 1.5 * x0

function obj(x)
    return x[1]^2 + x[2]^2
end

x = sample(5, lb, ub, SobolSample())

y = [obj(x[i]) for i = 1:length(x)]

surrogate_init = RadialBasis(x, y, lb, ub; rad=linearRadial())

x_opt, y_opt = surrogate_optimize(
    obj,
    SRBF(),
    lb,
    ub,
    surrogate_init,
    SobolSample();
)

# surrogate_init = Kriging(x, y, lb, ub)

# x_opt, y_opt = surrogate_optimize(
#     obj,
#     EI(),
#     lb,
#     ub,
#     surrogate_init,
#     SobolSample();
# )

[ChrisRackauckas]

That's only for certain choices of surrogates and certain choices of (x,y) right? It's definitely possible to get a singular exception with some cases which may need more points and have a degenerate solution if say everything falls on a line.