#julia #mathematical-optimization
#julia #математическая оптимизация
Вопрос:
У меня есть этот код для оптимизации функции с помощью IPNewton
метода ( error.jl
):
import Optim
"""
Generate a matrix of constants used in computation
"""
function get_const(x::Vector{Float64}, sigma::Vector{Float64})::Array{Float64, 2}
exp.(-x'.^2 ./ (2 .* sigma.^2)) ./ (sigma .* sqrt(2 * π))
end
# Log likelihood for mixture model
log_likelihood(p, C::Array{Float64, 2}) = sum(log.(p' * C))
"""
Constraint: all probabilities (ps) must sum to 1
"""
function constraint!(c, ps)::typeof(c)
c[1] = sum(ps)
c
end
N = 100
x = range(-1, 1, length=1000) |> collect
sigma = range(0.001, 2, length=N) |> collect
C = get_const(x, sigma)
constraints = Optim.TwiceDifferentiableConstraints(
constraint!,
fill(0, N), fill(1, N), # 0 <= (each probability) <= 1
fill(1, N), fill(1, N) # 1 <= constraint(p) <= 1 (probabilities sum to 1)
)
p0 = fill(1, N) / N # initial guess == equal probabilities
res = Optim.optimize(
ps -> -log_likelihood(ps, C), # want to MAXIMIZE, so negate
constraints, p0,
Optim.IPNewton()
)
Project.toml
:
[deps]
Optim = "429524aa-4258-5aef-a3af-852621145aeb"
Версия Julia:
forcebru@thing ~/test> julia --version
julia version 1.5.3
Сообщение об ошибке:
forcebru@thing ~/test> julia error.jl
ERROR: LoadError: InexactError: Int64(0.01)
Stacktrace:
[1] Int64 at ./float.jl:710 [inlined]
[2] convert at ./number.jl:7 [inlined]
[3] setindex! at ./array.jl:847 [inlined]
[4] _unsafe_copyto!(::Array{Int64,1}, ::Int64, ::Array{Float64,1}, ::Int64, ::Int64) at ./array.jl:257
[5] unsafe_copyto! at ./array.jl:311 [inlined]
[6] _copyto_impl! at ./array.jl:335 [inlined]
[7] copyto! at ./array.jl:321 [inlined]
[8] copyto! at ./array.jl:347 [inlined]
[9] finite_difference_jacobian!(::Array{Float64,2}, ::typeof(constraint!), ::Array{Float64,1}, ::FiniteDiff.JacobianCache{Array{Int64,1},Array{Int64,1},Array{Int64,1},UnitRange{Int64},Nothing,Val{:central}(),Int64}, ::Nothing; relstep::Float64, absstep::Float64, colorvec::UnitRange{Int64}, sparsity::Nothing, dir::Bool) at /Users/forcebru/.julia/packages/FiniteDiff/jLwWI/src/jacobians.jl:338
[10] finite_difference_jacobian!(::Array{Float64,2}, ::Function, ::Array{Float64,1}, ::FiniteDiff.JacobianCache{Array{Int64,1},Array{Int64,1},Array{Int64,1},UnitRange{Int64},Nothing,Val{:central}(),Int64}, ::Nothing) at /Users/forcebru/.julia/packages/FiniteDiff/jLwWI/src/jacobians.jl:334 (repeats 2 times)
[11] jac! at /Users/forcebru/.julia/packages/NLSolversBase/QPnui/src/objective_types/constraints.jl:298 [inlined]
[12] initial_state(::Optim.IPNewton{typeof(Optim.backtrack_constrained_grad),Symbol}, ::Optim.Options{Float64,Nothing}, ::NLSolversBase.TwiceDifferentiable{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::NLSolversBase.TwiceDifferentiableConstraints{typeof(constraint!),NLSolversBase.var"#jac!#126"{typeof(constraint!),FiniteDiff.JacobianCache{Array{Int64,1},Array{Int64,1},Array{Int64,1},UnitRange{Int64},Nothing,Val{:central}(),Int64}},NLSolversBase.var"#con_hess!#130"{Int64,Array{Int64,2},Array{Int64,3},NLSolversBase.var"#jac_vec!#129"{Int64,Int64},FiniteDiff.JacobianCache{Array{Int64,1},Array{Int64,1},Array{Int64,1},UnitRange{Int64},Nothing,Val{:central}(),Int64}},Int64}, ::Array{Float64,1}) at /Users/forcebru/.julia/packages/Optim/D7azp/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl:135
[13] optimize(::NLSolversBase.TwiceDifferentiable{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::NLSolversBase.TwiceDifferentiableConstraints{typeof(constraint!),NLSolversBase.var"#jac!#126"{typeof(constraint!),FiniteDiff.JacobianCache{Array{Int64,1},Array{Int64,1},Array{Int64,1},UnitRange{Int64},Nothing,Val{:central}(),Int64}},NLSolversBase.var"#con_hess!#130"{Int64,Array{Int64,2},Array{Int64,3},NLSolversBase.var"#jac_vec!#129"{Int64,Int64},FiniteDiff.JacobianCache{Array{Int64,1},Array{Int64,1},Array{Int64,1},UnitRange{Int64},Nothing,Val{:central}(),Int64}},Int64}, ::Array{Float64,1}, ::Optim.IPNewton{typeof(Optim.backtrack_constrained_grad),Symbol}, ::Optim.Options{Float64,Nothing}) at /Users/forcebru/.julia/packages/Optim/D7azp/src/multivariate/solvers/constrained/ipnewton/interior.jl:228
[14] optimize(::Function, ::NLSolversBase.TwiceDifferentiableConstraints{typeof(constraint!),NLSolversBase.var"#jac!#126"{typeof(constraint!),FiniteDiff.JacobianCache{Array{Int64,1},Array{Int64,1},Array{Int64,1},UnitRange{Int64},Nothing,Val{:central}(),Int64}},NLSolversBase.var"#con_hess!#130"{Int64,Array{Int64,2},Array{Int64,3},NLSolversBase.var"#jac_vec!#129"{Int64,Int64},FiniteDiff.JacobianCache{Array{Int64,1},Array{Int64,1},Array{Int64,1},UnitRange{Int64},Nothing,Val{:central}(),Int64}},Int64}, ::Array{Float64,1}, ::Optim.IPNewton{typeof(Optim.backtrack_constrained_grad),Symbol}, ::Optim.Options{Float64,Nothing}; inplace::Bool, autodiff::Symbol) at /Users/forcebru/.julia/packages/Optim/D7azp/src/multivariate/optimize/interface.jl:148
[15] optimize(::Function, ::NLSolversBase.TwiceDifferentiableConstraints{typeof(constraint!),NLSolversBase.var"#jac!#126"{typeof(constraint!),FiniteDiff.JacobianCache{Array{Int64,1},Array{Int64,1},Array{Int64,1},UnitRange{Int64},Nothing,Val{:central}(),Int64}},NLSolversBase.var"#con_hess!#130"{Int64,Array{Int64,2},Array{Int64,3},NLSolversBase.var"#jac_vec!#129"{Int64,Int64},FiniteDiff.JacobianCache{Array{Int64,1},Array{Int64,1},Array{Int64,1},UnitRange{Int64},Nothing,Val{:central}(),Int64}},Int64}, ::Array{Float64,1}, ::Optim.IPNewton{typeof(Optim.backtrack_constrained_grad),Symbol}, ::Optim.Options{Float64,Nothing}) at /Users/forcebru/.julia/packages/Optim/D7azp/src/multivariate/optimize/interface.jl:147 (repeats 2 times)
[16] top-level scope at /Users/forcebru/test/error.jl:27
[17] include(::Function, ::Module, ::String) at ./Base.jl:380
[18] include(::Module, ::String) at ./Base.jl:368
[19] exec_options(::Base.JLOptions) at ./client.jl:296
[20] _start() at ./client.jl:506
in expression starting at /Users/forcebru/test/error.jl:27
forcebru@thing ~/test [1]>
Итак … InexactError: Int64(0.01)
? И, похоже, это также происходит изнутри Optim
?
Я понимаю, что InexactError
здесь означает, что Джулия не смогла преобразовать 0.01
в целое число, что имеет смысл. Но я понятия не имею, откуда это 0.01
вообще взялось! Как узнать, откуда происходит? Что не так с этим кодом и что можно сделать, чтобы это исправить?
РЕДАКТИРОВАТЬ: я заметил, что 0.01
должен быть элементом p0 = fill(1, N) / N
, потому что, если я установил N = 50
, ошибка становится InexactError: Int64(0.02)
, где 0.02 == 1/N
. Но почему он пытается преобразовать его в целое число??
Ответ №1:
После некоторого интенсивного просмотра этих частей сообщения об ошибке:
[8] copyto! at ./array.jl:347 [inlined]
[9] finite_difference_jacobian!(::Array{Float64,2}, ::typeof(constraint!), ::Array{Float64,1}, ::FiniteDiff.JacobianCache{Array{Int64,1},Array{Int64,1},Array{Int64,1},UnitRange{Int64},Nothing,Val{:central}(),Int64}, ::Nothing; relstep::Float64, absstep::Float64, colorvec::UnitRange{Int64}, sparsity::Nothing, dir::Bool) at /Users/forcebru/.julia/packages/FiniteDiff/jLwWI/src/jacobians.jl:338
...
[15] optimize(::Function, ::NLSolversBase.TwiceDifferentiableConstraints{typeof(constraint!),NLSolversBase.var"#jac!#126"{typeof(constraint!),FiniteDiff.JacobianCache{Array{Int64,1},Array{Int64,1},Array{Int64,1},UnitRange{Int64},Nothing,Val{:central}(),Int64}},NLSolversBase.var"#con_hess!#130"{Int64,Array{Int64,2},Array{Int64,3},NLSolversBase.var"#jac_vec!#129"{Int64,Int64},FiniteDiff.JacobianCache{Array{Int64,1},Array{Int64,1},Array{Int64,1},UnitRange{Int64},Nothing,Val{:central}(),Int64}},Int64}, ::Array{Float64,1}, ::Optim.IPNewton{typeof(Optim.backtrack_constrained_grad),Symbol}, ::Optim.Options{Float64,Nothing}) at /Users/forcebru/.julia/packages/Optim/D7azp/src/multivariate/optimize/interface.jl:147 (repeats 2 times)
[16] top-level scope at /Users/forcebru/test/error.jl:27
…Я видел, что FiniteDiff.JacobianCache
для ограничений было выведено, что они должны быть параметризованы по Int64
:
FiniteDiff.JacobianCache{
Array{Int64,1},
Array{Int64,1},
Array{Int64,1},
UnitRange{Int64},
Nothing,
Val{:central}(),
Int64
}
… что довольно странно, потому что я явно хочу оптимизировать по сравнению с реалами.
Оказывается, что в этой части кода:
constraints = Optim.TwiceDifferentiableConstraints(
constraint!,
fill(0, N), fill(1, N), # 0 <= (each probability) <= 1
fill(1, N), fill(1, N) # 1 <= constraint(p) <= 1 (probabilities sum to 1)
)
все fill(0, N)
друзья и являются целыми числами, потому 0
что это целое число. Похоже, это привело к попытке преобразования из float в integer .
Я изменил этот код на чтение:
constraints = Optim.TwiceDifferentiableConstraints(
constraint!,
fill(0., N), fill(1., N), # 0 <= (each probability) <= 1
fill(1., N), fill(1., N) # 1 <= constraint(p) <= 1 (probabilities sum to 1)
)
… и теперь ошибки нет (хотя алгоритм не сходится, но это другая проблема).
Комментарии:
1. Было бы неплохо открыть проблему, отмечая такие вещи. Вероятно, Optim следует
float(x)
куда-нибудь позвонить, чтобы очистить такой ввод.2. @mcabbott, на самом деле, я изначально начал писать это как проблему, но потом подумал, что мне просто не хватает чего-то простого, потому что я новичок в Julia, поэтому вместо этого спросил здесь. Действительно, ошибка, возникающая в недрах библиотеки, выглядит как ошибка в самой библиотеке.
3. Ошибка может быть слишком сильным словом, но было бы удобно перехватить и исправить целочисленный ввод, а не сбой. Причина, по которой вы получаете сбой из глубины, заключается в том, что типы часто остаются общими, так что они могут с радостью специализироваться на других странных числах.