mariohsouto / proxsdp.jl Goto Github PK
View Code? Open in Web Editor NEWSemidefinite programming optimization solver
License: MIT License
Semidefinite programming optimization solver
License: MIT License
No matter what we do: opencollab/arpack-ng#218
Originally: scipy/scipy#10374
Julia alternatives:
https://github.com/haampie/ArnoldiMethod.jl
however, there are no symmetric specialized methods (Lanczos):
https://github.com/haampie/ArnoldiMethod.jl/issues/92
https://github.com/Jutho/KrylovKit.jl
Has Lanczos!
This issue is used to trigger TagBot; feel free to unsubscribe.
If you haven't already, you should update your TagBot.yml
to include issue comment triggers.
Please see this post on Discourse for instructions and more details.
If you'd like for me to do this for you, comment TagBot fix
on this issue.
I'll open a PR within a few hours, please be patient!
I was trying ProxSDP in solving a convex SDP whose solution rank is 1, my optimization is:
min (trace(QZ))
subject to: some equality constraints on entries of Z.
Where Q is a given cost matrix and Z is my PSD decision variable.
This optimization gives rank 1 solution in CVX using SeDuMi (Matlab), and also gives rank 1 solution in Julia using Mosek solver.
However, ProxSDP does not give rank 1 solution, and the resulting Z^* is incorrect, although I did notice that ProxSDP is much faster than Mosek.
So I am wondering, is ProxSDP guaranteed to find the optimal solution of an SDP? (Just like conventional Interior Points Methods)?
julia> model = Model(with_optimizer(ProxSDP.Optimizer))
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: SolverName() attribute not implemented by the optimizer.
Check there is enough convergence lanczos:
https://github.com/JuliaLinearAlgebra/Arpack.jl/blob/a7cdb6d7f19f076f5fadd8b58a9c5a061c48322f/src/Arpack.jl#L191
using JuMP, ProxSDP
model = Model()
@variable(model, k[i=1:4])
@constraint(model, [k[1], k[2], k[3] + k[4]] in SecondOrderCone())
JuMP.optimize!(model, with_optimizer(ProxSDP.Optimizer, log_verbose=true))
The resulting error is
ERROR: MathOptInterface.UnsupportedConstraint{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.SecondOrderCone}: `MathOptInterface.VectorAffineFunction{Float64}`-in-`MathOptInterface.SecondOrderCone` constraints is not supported by the the model.
Stacktrace:
[1] bridge_type(::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{ProxSDP.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP.JuMPMOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}, ::Type{MathOptInterface.VectorAffineFunction{Float64}}, ::Type{MathOptInterface.SecondOrderCone}) at /home/guilhermebodin/.julia/packages/MathOptInterface/A2LvE/src/Bridges/lazybridgeoptimizer.jl:110
[2] concrete_bridge_type(::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{ProxSDP.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP.JuMPMOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}, ::Type{MathOptInterface.VectorAffineFunction{Float64}}, ::Type{MathOptInterface.SecondOrderCone}) at /home/guilhermebodin/.julia/packages/MathOptInterface/A2LvE/src/Bridges/bridgeoptimizer.jl:59
[3] add_constraint(::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{ProxSDP.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP.JuMPMOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}, ::MathOptInterface.VectorAffineFunction{Float64}, ::MathOptInterface.SecondOrderCone) at /home/guilhermebodin/.julia/packages/MathOptInterface/A2LvE/src/Bridges/bridgeoptimizer.jl:315
[4] copyconstraints!(::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{ProxSDP.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP.JuMPMOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}, ::MathOptInterface.Utilities.UniversalFallback{JuMP.JuMPMOIModel{Float64}}, ::Bool, ::MathOptInterface.Utilities.IndexMap, ::Type{MathOptInterface.VectorAffineFunction{Float64}}, ::Type{MathOptInterface.SecondOrderCone}) at /home/guilhermebodin/.julia/packages/MathOptInterface/A2LvE/src/Utilities/copy.jl:158
[5] default_copy_to(::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{ProxSDP.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP.JuMPMOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}, ::MathOptInterface.Utilities.UniversalFallback{JuMP.JuMPMOIModel{Float64}},::Bool) at /home/guilhermebodin/.julia/packages/MathOptInterface/A2LvE/src/Utilities/copy.jl:200
[6] #automatic_copy_to#55 at /home/guilhermebodin/.julia/packages/MathOptInterface/A2LvE/src/Utilities/copy.jl:15 [inlined]
[7] #automatic_copy_to at ./none:0 [inlined]
[8] #copy_to#1 at /home/guilhermebodin/.julia/packages/MathOptInterface/A2LvE/src/Bridges/bridgeoptimizer.jl:91 [inlined]
[9] (::getfield(MathOptInterface, Symbol("#kw##copy_to")))(::NamedTuple{(:copy_names,),Tuple{Bool}}, ::typeof(MathOptInterface.copy_to), ::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{ProxSDP.Optimizer,MathOptInterface.Utilities.UniversalFallback{JuMP.JuMPMOIModel{Float64}}},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Bridges.AllBridgedConstraints{Float64}}}, ::MathOptInterface.Utilities.UniversalFallback{JuMP.JuMPMOIModel{Float64}}) at ./none:0
[10] attach_optimizer(::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer,MathOptInterface.Utilities.UniversalFallback{JuMP.JuMPMOIModel{Float64}}}) at /home/guilhermebodin/.julia/packages/MathOptInterface/A2LvE/src/Utilities/cachingoptimizer.jl:125
[11] attach_optimizer(::Model) at /home/guilhermebodin/.julia/packages/JuMP/xQHBA/src/optimizer_interface.jl:32
[12] #optimize!#77(::Bool, ::Bool, ::Function, ::Model, ::JuMP.OptimizerFactory) at /home/guilhermebodin/.julia/packages/JuMP/xQHBA/src/optimizer_interface.jl:115
[13] optimize!(::Model, ::JuMP.OptimizerFactory) at /home/guilhermebodin/.julia/packages/JuMP/xQHBA/src/optimizer_interface.jl:101
[14] top-level scope at none:0
Use LAPACK LA_SPEVX (http://www.netlib.org/lapack95/DOC/la_spevx.txt) function to compute a truncated eigenvalue decomposition of a lower triangular matrix in packed format.
Run full eig decomps after convergence
That is, instead of assuming one PSD constraint X in PSDCone()
, we can have X_1 in PSDCone()
, X_2 in PSDCone()
, etc. Is this how ProxSDP is implemented? Or does it assume one matrix block X in PSDCone()
?
I expect this question to be related to the issue:
Chordal decomposition issue
I really want to try ProxSDP for solving low-rank SDPs in my research. But I am not able to run ProxSDP on my Mac OS with Julia 1.0, or my Linux PC with Julia 0.6.
I am wondering what are the correct (recommended) ways to be able to run your solver?
After solving a problem:
SDPproblem = # defined through JuMP
with_ProxSDP = with_optimizer(ProxSDP.Optimizer, tol_primal=1e-10, tol_dual=1e-10, log_verbose=true)
optimize!(SDPproblem, with_ProxSDP)
...
----------------------------------------------------------------------
ProxSDP failed to converge in 70.66 seconds
Primal objective = -0.0
Dual objective = 0.00026
Duality gap (%) = 3.133284148680549e13 %
======================================================================
objective_value(SDPproblem)
OptimizeNotCalled()
...
Additional info:
It is a SDP problem which should (for theoretical reasons) have a low-rank solution;
As I define it, it's a feasibility problem, but I add a dummy optimization variable to maximize and bound it by 0.0 from above.
ProxSDP is failing geomean2 tests that were added jump-dev/MathOptInterface.jl#962:
https://travis-ci.com/blegat/SolverTests/jobs/268551133
remove this sentence:
"The main caveat is that currently ProxSDP must have one and only one PSD variable, no other variables are allowed." because its not true anymore
remove MOI version from manual
I was wondering if as a pure-Julia SDP solver, would it be possible to use something like BigFloat
's with ProxSDP.jl
? I know JuMP automatically converts inputs to Float64
but I thought maybe it would be possible to use the library directly without going through JuMP. I thought I would ask first if it were possible first, before spending too much time on it though. I understand it would surely be slow, but there are some quite small SDPs that it would still be interesting to me to check some things at high precision.
jump-dev/MathOptInterface.jl#528
Will enable F
in S
for SOC and SDCone.
Consequently we get RootDet and GeoMean for free
Follow the steps of: jump-dev/SCS.jl#144
Just add an “else” in parse arg function
running moi_randsdp
with varbounds = false
leads to error on julia >= 1.6.
some values of v
in psd_vec_to_square
are +/-Inf
check by adding:
if isnan(v[cont])
@show v
elseif !isfinite(v[cont])
@show 1, v
end
in that function
From the MOI wrapper, I see that you are transforming the PSD set from "upper triangular column-wise" (the MOI definition) to "lower triangular column-wise" (the SCS definition). Is there any reason not to use the MOI definition ?
As noted here:
The advantage of the first format is the mapping between the (i, j) matrix indices and the k index of the vectorized form. It is simpler and does not depend on the side dimension of the matrix.
so there is a small advantage for the MOI version but I see no advantage for the SCS one.
See here: https://ericphanson.github.io/ConvexTests.jl/dev/ProxSDP.html (and see here for a discourse post introducing that repo, ConvexTests.jl). Some of the failures might just be settings and things like that; e.g., with COSMO most of the failures were max-iteration related, but it would be great to figure out what settings can solve a variety of problems to 1e-3
, and it's possible the tests are surfacing bugs as well.
ProxSDP should exclude relentr test:
https://travis-ci.com/blegat/SolverTests/jobs/281969930
Lanczos is specialized for symmetric matrices, however, it might be possible that Arnoldi, although less efficient might be more robust.
Running the example on the front page produces the following error:
julia> JuMP.optimize!(model)
ERROR: error compiling #optimize!#77: error compiling optimize!: error compiling chambolle_pock: error compiling primal_step!: error compiling sdp_cone_projection!: error compiling eig!: error compiling _AUPD!: error compiling saupd: could not load library "/home/affans/.julia/packages/Arpack/qof0w/deps/usr/lib/libarpack.so"
libopenblas64_.so.0: cannot open shared object file: No such file or directory
Stacktrace:
[1] optimize! at /home/affans/.julia/packages/JuMP/jnmGG/src/optimizer_interface.jl:105 [inlined] (repeats 2 times)
[2] top-level scope at none:0
I cannot reproduce the timings from the sdplib
instances of this repo's preprint. For example, ProxSDP
terminates after 100k iterations without solving gpp124-1
to 1e-3
tolerance.
Update: See comments below
Download gpp124-1.jld2 and run
using LinearAlgebra, SparseArrays
using ProxSDP, JuMP
using JLD2
function solve_sdpa_jump(c, F; solver=ProxSDP.Optimizer, kwargs...)
model = Model(with_optimizer(solver; kwargs...))
m = length(c); n = size(F[1], 1)
@variable(model, x[1:m])
A = hcat([F[i+1][:] for i = 1:m]...)
b = -Vector(F[1][:])
@constraint(model, A*x + b in MOI.PositiveSemidefiniteConeSquare(n))
@objective(model, Min, dot(c, x))
JuMP.optimize!(model)
@show JuMP.termination_status(model)
return value.(x), JuMP.objective_value(model)
end
@load "gpp124-1.jld2" c F
solution, solution_objective = solve_sdpa_jump(c, F, solver=ProxSDP.Optimizer, log_verbose=true)
In the example above ProxSDP.jl v1.0.0
does not converge after 100 thousand iterations. In contrast, running MOSEK.jl 0.9.11
via:
using MosekTools
solution, solution_objective = solve_sdpa_jump(c, F, solver=Mosek.Optimizer)
converges to a solution with objective value very close to the one provided by sdplib
's readme (-7.3431).
Let me know if you think that there is some inefficiency in the above code. Are there any plans for publishing the code used in the experiments?
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.