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?
1 Like
Check https://github.com/JuliaNonconvex/NonconvexSemidefinite.jl. You can constrain matrix-valued functions of the decision variables to be positive semidefinite. There is an example in the tests https://github.com/JuliaNonconvex/NonconvexSemidefinite.jl/blob/master/test/runtests.jl.
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 https://github.com/TuringLang/DistributionsAD.jl by just calling using Distributions, DistributionsAD. If neither of these solutions work, please report here with the output of ]st Distributions.
1 Like
Here are some other packages that may be useful:
2 Likes
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.
See here.
@mohamed82008 Do you have any plans or thoughts about adding a global optimization capability for Nonlinear SDPs? For instance, a Branch and Bound Solver similar to YALMIP’s BMIBNB SDP cones in BMIBNB - YALMIP which at least attempts to solve Nonlinear SDPs to global optimality even though it might only succeed on relatively small and/or “easy” problems. The local Nonlinear SDP solver could be used as an upper solver for the Branch and Bound global solver. Even if not solved to global optimality, a global solver might produce a better solution than what you would find using a local optimizer, together with a guarantee on max objective difference from global optimum, which of course is helped by generating good lower solver problems to tighten the lower bound.
Yes! If it has a Julia wrapper, I will be happy to wrap it in NonconvexSemidefinite.
No Julia wrapper. The BMIBNB solver is written by YALMIP developer Johan Lofberg, and included in the YALMIP distribution. So it is basically MATLAB code (although much is YALMIP code, which is really MATLAB code), but calls various solvers to solve subproblems. The YALMIP code, to include the BMIBNB code is freely downloadable. I think if you were to incorporate stuff very directly from that code, you would want to get Johan’s concurrence.
As for whether a wrapper would be possible, I have no idea.
One other thing which might be kind of cool would be incorporating Quantum Relative Entropy into Nonlinear SDP. So Quantum Relative Entropy would be a convex cone as part of an overall non-convex problem. For example, non-convex quadratic objective with a Quantum Relative Entropy constraint. And maybe (or maybe not) there would be some nonlinear SDP constraints in there as well.
I can try implementing the algorithm but I have little scope for that in the next few months since NonconvexSemidefinite does the job for now albeit not in the best way possible. Combine it with deflation from https://www.researchgate.net/publication/358151220_Simplifying_deflation_for_non-convex_optimization_with_applications_in_Bayesian_inference_and_topology_optimization/stats and you can get multiple local minima. Either way, it’s helpful if you open an issue in Nonconvex.jl and post all the relevant references for anyone willing to give it a try.