I want to do the following optimization:
\min E \ \ \text{s.t.} \ \ x \in \mathbb{R}, \ \ \text{$M$ is semidefinite},
where M_{ij} is a polynominal of x and E. Is there a Julia package suitable for this problem?
Check GitHub - JuliaNonconvex/NonconvexSemidefinite.jl: Nonlinear semi-definite programming algorithms for Nonconvex.jl. You can constrain matrix-valued functions of the decision variables to be positive semidefinite. There is an example in the tests NonconvexSemidefinite.jl/runtests.jl at master · JuliaNonconvex/NonconvexSemidefinite.jl · GitHub.
1 Like
I got an error when running your unittest … on the 26th line of the test, I got
Single semi-definite constraint: Error During Test at d:\Projects\numerical-boostrap\oscillator-simple-prototype\sdp-test.jl:26
Got exception outside of a @test
Need an adjoint for constructor Distributions.EachVariate{1,Array{Float64,2},Tuple{Base.OneTo{Int64}},SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},1}. Gradient is of type Array{Array{Float64,1},1}
Stacktrace:
[1] error(::String) at .\error.jl:33
[2] (::Zygote.Jnew{Distributions.EachVariate{1,Array{Float64,2},Tuple{Base.OneTo{Int64}},SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},1},Nothing,false})(::Array{Array{Float64,1},1}) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\lib\lib.jl:324
[3] (::Zygote.var"#1781#back#223"{Zygote.Jnew{Distributions.EachVariate{1,Array{Float64,2},Tuple{Base.OneTo{Int64}},SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},1},Nothing,false}})(::Array{Array{Float64,1},1}) at C:\Users\wujin\.julia\packages\ZygoteRules\AIbCs\src\adjoint.jl:67
[4] EachVariate at C:\Users\wujin\.julia\packages\Distributions\9Albf\src\eachvariate.jl:4 [inlined]
[5] (::typeof(∂(Distributions.EachVariate{1,Array{Float64,2},Tuple{Base.OneTo{Int64}},SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},1})))(::Array{Array{Float64,1},1}) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\compiler\interface2.jl:0
[6] EachVariate at C:\Users\wujin\.julia\packages\Distributions\9Albf\src\eachvariate.jl:11 [inlined]
[7] (::typeof(∂(Distributions.EachVariate{1,P,A,T,N} where N where T where A where P)))(::Array{Array{Float64,1},1}) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\compiler\interface2.jl:0
[8] eachvariate at C:\Users\wujin\.julia\packages\Distributions\9Albf\src\eachvariate.jl:31 [inlined]
[9] (::typeof(∂(eachvariate)))(::Array{Array{Float64,1},1}) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\compiler\interface2.jl:0
[10] loglikelihood at C:\Users\wujin\.julia\packages\Distributions\9Albf\src\common.jl:440 [inlined]
[11] f at d:\Projects\numerical-boostrap\oscillator-simple-prototype\sdp-test.jl:35 [inlined]
[12] (::typeof(∂(λ)))(::Float64) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\compiler\interface2.jl:0
[13] #207 at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\lib\lib.jl:203 [inlined]
[14] (::Zygote.var"#1747#back#209"{Zygote.var"#207#208"{typeof(∂(λ)),Tuple{Tuple{Nothing}}}})(::Float64) at C:\Users\wujin\.julia\packages\ZygoteRules\AIbCs\src\adjoint.jl:67
[15] #_#8 at C:\Users\wujin\.julia\packages\NonconvexCore\YlMJF\src\functions\functions.jl:170 [inlined] [16] (::typeof(∂(#_#8)))(::Float64) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\compiler\interface2.jl:0
[17] (::Zygote.var"#207#208"{typeof(∂(#_#8)),Tuple{Tuple{Nothing,Nothing},Tuple{Nothing}}})(::Float64) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\lib\lib.jl:203
[18] #1747#back at C:\Users\wujin\.julia\packages\ZygoteRules\AIbCs\src\adjoint.jl:67 [inlined]
[19] Objective at C:\Users\wujin\.julia\packages\NonconvexCore\YlMJF\src\functions\functions.jl:170 [inlined]
[20] (::typeof(∂(λ)))(::Float64) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\compiler\interface2.jl:0
[21] #133 at C:\Users\wujin\.julia\packages\NonconvexCore\YlMJF\src\models\vec_model.jl:89 [inlined]
[22] (::typeof(∂(λ)))(::Float64) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\compiler\interface2.jl:0
[23] #207 at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\lib\lib.jl:203 [inlined]
[24] (::Zygote.var"#1747#back#209"{Zygote.var"#207#208"{typeof(∂(λ)),Tuple{Tuple{Nothing}}}})(::Float64) at C:\Users\wujin\.julia\packages\ZygoteRules\AIbCs\src\adjoint.jl:67
[25] #_#8 at C:\Users\wujin\.julia\packages\NonconvexCore\YlMJF\src\functions\functions.jl:170 [inlined] [26] (::typeof(∂(#_#8)))(::Float64) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\compiler\interface2.jl:0
[27] (::Zygote.var"#207#208"{typeof(∂(#_#8)),Tuple{Tuple{Nothing,Nothing},Tuple{Nothing}}})(::Float64) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\lib\lib.jl:203
[28] #1747#back at C:\Users\wujin\.julia\packages\ZygoteRules\AIbCs\src\adjoint.jl:67 [inlined]
[29] Objective at C:\Users\wujin\.julia\packages\NonconvexCore\YlMJF\src\functions\functions.jl:170 [inlined]
[30] (::typeof(∂(λ)))(::Float64) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\compiler\interface2.jl:0
[31] CountingFunction at C:\Users\wujin\.julia\packages\NonconvexCore\YlMJF\src\functions\counting_function.jl:9 [inlined]
[32] (::typeof(∂(λ)))(::Float64) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\compiler\interface2.jl:0
[33] (::Zygote.var"#55#56"{typeof(∂(λ))})(::Float64) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\compiler\interface.jl:41
[34] gradient(::Function, ::Array{Float64,1}) at C:\Users\wujin\.julia\packages\Zygote\FPUm3\src\compiler\interface.jl:76
[35] (::NonconvexIpopt.var"#eval_grad_f#13"{NonconvexCore.CountingFunction{NonconvexCore.Objective{NonconvexCore.var"#133#140"{Model{Array{Any,1}},NonconvexCore.Unflatten{Array{Array{Float64,1},1},NonconvexCore.var"#Vector_from_vec#51"{Array{Array{Float64,1},1},Array{Array{Float64,1},1},Array{typeof(identity),1}}}},Base.RefValue{Float64},Set{Symbol}}}})(::Array{Float64,1}, ::Array{Float64,1}) at C:\Users\wujin\.julia\packages\NonconvexIpopt\dzR73\src\ipopt.jl:126
[36] eval_grad_f_wrapper(::Int32, ::Ptr{Float64}, ::Int32, ::Ptr{Float64}, ::Ptr{Nothing}) at C:\Users\wujin\.julia\packages\Ipopt\vtrOr\src\Ipopt.jl:163
[37] solveProblem(::Ipopt.IpoptProblem) at C:\Users\wujin\.julia\packages\Ipopt\vtrOr\src\Ipopt.jl:532
[38] optimize!(::NonconvexIpopt.IpoptWorkspace{NonconvexCore.VecModel{Array{Float64,1}},Ipopt.IpoptProblem,Array{Float64,1},IpoptOptions{NamedTuple{(:max_iter, :hessian_approximation, :jac_c_constant, :jac_d_constant),Tuple{Int64,String,String,String}}},Base.RefValue{Int64}}) at C:\Users\wujin\.julia\packages\NonconvexIpopt\dzR73\src\ipopt.jl:52
[39] #optimize#131 at C:\Users\wujin\.julia\packages\NonconvexCore\YlMJF\src\models\vec_model.jl:73 [inlined]
[40] optimize!(::NonconvexSemidefinite.SDPBarrierWorkspace{NonconvexCore.VecModel{Array{Float64,1}},Array{Float64,1},SDPBarrierOptions{Float64,Float64,Int64,IpoptOptions{NamedTuple{(:max_iter, :hessian_approximation, :jac_c_constant, :jac_d_constant),Tuple{Int64,String,String,String}}},Bool},IpoptAlg}) at C:\Users\wujin\.julia\packages\NonconvexSemidefinite\W8Ncn\src\NonconvexSemidefinite.jl:162
[41] #optimize#131 at C:\Users\wujin\.julia\packages\NonconvexCore\YlMJF\src\models\vec_model.jl:73 [inlined]
[42] optimize(::Model{Array{Any,1}}, ::SDPBarrierAlg{IpoptAlg}, ::Array{Array{Float64,1},1}; kwargs::Base.Iterators.Pairs{Symbol,SDPBarrierOptions{Float64,Float64,Int64,IpoptOptions{NamedTuple{(:max_iter, :hessian_approximation, :jac_c_constant, :jac_d_constant),Tuple{Int64,String,String,String}}},Bool},Tuple{Symbol},NamedTuple{(:options,),Tuple{SDPBarrierOptions{Float64,Float64,Int64,IpoptOptions{NamedTuple{(:max_iter, :hessian_approximation, :jac_c_constant, :jac_d_constant),Tuple{Int64,String,String,String}}},Bool}}}}) at C:\Users\wujin\.julia\packages\NonconvexCore\YlMJF\src\common.jl:233
[43] top-level scope at d:\Projects\numerical-boostrap\oscillator-simple-prototype\sdp-test.jl:50
[44] top-level scope at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\Test\src\Test.jl:1115
[45] top-level scope at d:\Projects\numerical-boostrap\oscillator-simple-prototype\sdp-test.jl:28
[46] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091
[47] invokelatest(::Any, ::Any, ::Vararg{Any,N} where N; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at .\essentials.jl:710
[48] invokelatest(::Any, ::Any, ::Vararg{Any,N} where N) at .\essentials.jl:709
[49] inlineeval(::Module, ::String, ::Int64, ::Int64, ::String; softscope::Bool) at c:\Users\wujin\.vscode\extensions\julialang.language-julia-1.5.10\scripts\packages\VSCodeServer\src\eval.jl:211
[50] (::VSCodeServer.var"#58#62"{String,Int64,Int64,String,Module,Bool,Bool,VSCodeServer.ReplRunCodeRequestParams})() at c:\Users\wujin\.vscode\extensions\julialang.language-julia-1.5.10\scripts\packages\VSCodeServer\src\eval.jl:155
[51] withpath(::VSCodeServer.var"#58#62"{String,Int64,Int64,String,Module,Bool,Bool,VSCodeServer.ReplRunCodeRequestParams}, ::String) at c:\Users\wujin\.vscode\extensions\julialang.language-julia-1.5.10\scripts\packages\VSCodeServer\src\repl.jl:185
[52] (::VSCodeServer.var"#57#61"{String,Int64,Int64,String,Module,Bool,Bool,Bool,VSCodeServer.ReplRunCodeRequestParams})() at c:\Users\wujin\.vscode\extensions\julialang.language-julia-1.5.10\scripts\packages\VSCodeServer\src\eval.jl:153
[53] hideprompt(::VSCodeServer.var"#57#61"{String,Int64,Int64,String,Module,Bool,Bool,Bool,VSCodeServer.ReplRunCodeRequestParams}) at c:\Users\wujin\.vscode\extensions\julialang.language-julia-1.5.10\scripts\packages\VSCodeServer\src\repl.jl:36
[54] (::VSCodeServer.var"#56#60"{String,Int64,Int64,String,Module,Bool,Bool,Bool,VSCodeServer.ReplRunCodeRequestParams})() at c:\Users\wujin\.vscode\extensions\julialang.language-julia-1.5.10\scripts\packages\VSCodeServer\src\eval.jl:124
[55] with_logstate(::Function, ::Any) at .\logging.jl:408
[56] with_logger at .\logging.jl:514 [inlined]
[57] (::VSCodeServer.var"#55#59"{VSCodeServer.ReplRunCodeRequestParams})() at c:\Users\wujin\.vscode\extensions\julialang.language-julia-1.5.10\scripts\packages\VSCodeServer\src\eval.jl:201
[58] #invokelatest#1 at .\essentials.jl:710 [inlined]
[59] invokelatest(::Any) at .\essentials.jl:709
[60] macro expansion at c:\Users\wujin\.vscode\extensions\julialang.language-julia-1.5.10\scripts\packages\VSCodeServer\src\eval.jl:34 [inlined]
Test Summary: | Error Total
Single semi-definite constraint | 1 1
ERROR: LoadError: Some tests did not pass: 0 passed, 0 failed, 1 errored, 0 broken.
in expression starting at d:\Projects\numerical-boostrap\oscillator-simple-prototype\sdp-test.jl:26
Here the main problem seems to be
Need an adjoint for constructor Distributions.EachVariate{1,Array{Float64,2},Tuple{Base.OneTo{Int64}},SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},1}. Gradient is of type Array{Array{Float64,1},1}
Tests passed with an older version of Distributions. It’s possible you are using a version that’s not Zygote compatible. Try an older version of Distributions, or try using GitHub - TuringLang/DistributionsAD.jl: Automatic differentiation of Distributions using Tracker, Zygote, ForwardDiff and ReverseDiff by just calling using Distributions, DistributionsAD. If neither of these solutions work, please report here with the output of ]st Distributions.
1 Like
Changing
using Distributions, ChainRulesTestUtils, Random
into
using Distributions, DistributionsAD, ChainRulesTestUtils, Random
works and all test pass. Here are versions of the packages involved, which may be useful for anyone else facing the same problem:
[cdddcdb0] ChainRulesTestUtils v1.5.2
[31c24e10] Distributions v0.25.45
[ced4e74d] DistributionsAD v0.6.37
[bf347577] NonconvexIpopt v0.1.2
[9e21ff56] NonconvexSemidefinite v0.1.2 `https://github.com/JuliaNonconvex/NonconvexSemidefinite.jl#master`
[91a5bcdd] Plots v1.15.2
[9a3f8284] Random
There are still some issues. It seems that in some problems the semidefinite constraint isn’t working:
# Taken from https://github.com/JuliaNonconvex/NonconvexSemidefinite.jl/blob/master/test/runtests.jl
# Test cases of semidefinite programming
using NonconvexSemidefinite, NonconvexIpopt, LinearAlgebra, Test
using Distributions, DistributionsAD, ChainRulesTestUtils, Random
sd_constraint((x⁴, x²)) = [
1.0 0 x² - 1 ;
0 x² 0 ;
x² - 1 0 x⁴
]
# Objective function
f((x⁴, x²)) = 2x² + 3x⁴
model = Model(f)
x0 = [4, 2]
lbs = [0.0, 0.0]
ubs = [Inf, Inf]
addvar!(model, lbs, ubs)
add_sd_constraint!(model, sd_constraint)
alg = SDPBarrierAlg(sub_alg=IpoptAlg())
options = SDPBarrierOptions(sub_options=IpoptOptions(max_iter=200))
result = optimize(model, alg, x0, options = options)
minimum, minimizer, optimal_ind = result.minimum, result.minimizer, result.optimal_ind
m_mat_minimizer = sd_constraint(minimizer)
@show minimum
@show eigen(m_mat_minimizer).values
Result:
minimum = 0.0
(eigen(m_mat_minimizer)).values = [-0.6180339887498936, 0.0, 1.618033988749895]
3-element Array{Float64,1}:
-0.6180339887498936
0.0
1.618033988749895
It can be seen that sd_constraint has a negative eigenvalue after optimization.
What’s the output of sd_constraint(x0)? Is it positive definite?
Found the bug and fixed it locally, check the latest release again in an hour or so.
Use NonconvexSemidefinite 0.1.3.
2 Likes
There seems to be yet another problem, though I’m not sure whether it is a bug or feature and whether it comes from something in the auto differentiation code. Consider the following demo, where we keep a matrix to be semidefinite during the optimization:
using NonconvexSemidefinite, NonconvexIpopt, LinearAlgebra, Test
using Distributions, DistributionsAD, ChainRulesTestUtils, Random
K = 3
g = 1
Mij(x⁴, x², n) = begin
if n == 0
return 1.0
end
if isodd(n)
return 0.0
end
if n == 2
return x²
end
if n == 4
return x⁴
end
Mij(x⁴, x², n - 2) * x² + Mij(x⁴, x², n - 4)
end
sd_constraint((x⁴, x²)) = [Mij(x⁴, x², i + j) for i in 0 : K, j in 0 : K]
# Objective function
f((x⁴, x²)) = 2x² + 3x⁴
model = Model(f)
x0 = [1.0, 0.8]
lbs = [0.0, 0.0]
ubs = [Inf, Inf]
addvar!(model, lbs, ubs)
add_sd_constraint!(model, sd_constraint)
alg = SDPBarrierAlg(sub_alg=IpoptAlg())
options = SDPBarrierOptions(sub_options=IpoptOptions(max_iter=400))
result = optimize(model, alg, x0, options = options)
minimum, minimizer, optimal_ind = result.minimum, result.minimizer, result.optimal_ind
m_mat_minimizer = sd_constraint(minimizer)
@show minimum
@show eigen(m_mat_minimizer).values
The demo works well and finds the expected solution (just letting x⁴ == 0.0 and x² == 0.0). However, the following code
using NonconvexSemidefinite, NonconvexIpopt, LinearAlgebra, Test
using Distributions, DistributionsAD, ChainRulesTestUtils, Random
K = 3
g = 1
Mij(x⁴, x², n) = begin
if n == 0
return 1.0
end
if isodd(n)
return 0.0
end
if n == 2
return x²
end
if n == 4
return x⁴
end
(n - 4) * Mij(x⁴, x², n - 2) * x² + Mij(x⁴, x², n - 4)
end
sd_constraint((x⁴, x²)) = [Mij(x⁴, x², i + j) for i in 0 : K, j in 0 : K]
# Objective function
f((x⁴, x²)) = 2x² + 3x⁴
model = Model(f)
x0 = [1.0, 0.8]
lbs = [0.0, 0.0]
ubs = [Inf, Inf]
addvar!(model, lbs, ubs)
add_sd_constraint!(model, sd_constraint)
alg = SDPBarrierAlg(sub_alg=IpoptAlg())
options = SDPBarrierOptions(sub_options=IpoptOptions(max_iter=400))
result = optimize(model, alg, x0, options = options)
minimum, minimizer, optimal_ind = result.minimum, result.minimizer, result.optimal_ind
m_mat_minimizer = sd_constraint(minimizer)
@show minimum
@show eigen(m_mat_minimizer).values
throws a strange error, saying
LoadError: MethodError: no method matching getindex(::Nothing, ::Int64)
The only difference between the two versions is that the second version has an additional (n - 4) factor at the end of the definition of Mij. It can be verified by Mathematica that the second problem also has a solution at x⁴ == 0.0 and x² == 0.0.
I’m not sure what is going on here …
Plus: should we continue the discussion in a Github issue of JuliaNonconvex/NonconvexSemidefinite.jl?
Looks like a Zygote automatic differentiation issue. Please open an issue on NonconvexSemidefinite and I will be happy to help.