Nonlinear semidefinite programming recommendation

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

Here are some other packages that may be useful:

* PolyJuMP
* Polyopt.jl

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.

See here.