mariohsouto / proxsdp.jl Goto Github PK

View Code? Open in Web Editor NEW

92.0 92.0 12.0 12.32 MB

Semidefinite programming optimization solver

License: MIT License

Julia 98.93% MATLAB 1.07%

conic-programs convex-optimization julia optimization sdp semidefinite-programming solver

proxsdp.jl's People

Contributors

Stargazers

Watchers

Forkers

migarstka florencefeng batmanabcdefg stjordanis blegat abhrodutta juliatagbot sasadat curioustauseef huliangbing xiaogaogaoxiao odow

proxsdp.jl's Issues

Arpack is non-deterministic

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!

Generalize eigsolver selection and interface

type dispatch on eigsolver
unify interface
allow external eigsolver

TagBot trigger issue

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!

ProxSDP does not return same solution as Mosek, is ProxSDP guaranteed to solve an SDP?

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)?

Allow variable bounds in the primal step

SolverName not implemented

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 Lanczos convergence

Check there is enough convergence lanczos:
https://github.com/JuliaLinearAlgebra/Arpack.jl/blob/a7cdb6d7f19f076f5fadd8b58a9c5a061c48322f/src/Arpack.jl#L191

Bridge is not working for this example

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

Build a wrapper for LA_SPEVX Lapack function

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.

Add polishing step

Run full eig decomps after convergence

Does ProxSDP assume decomposition of X...?

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

Is ProxSDP ready to be used for Julia 1.0?

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?

Exponential cone

Document solver options

OptimizeNotCalled()

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.

MOI geomean2 tests

ProxSDP is failing geomean2 tests that were added jump-dev/MathOptInterface.jl#962:
https://travis-ci.com/blegat/SolverTests/jobs/268551133

Cleanup manual

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

high precision?

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.

Slack and Unslack bridge

jump-dev/MathOptInterface.jl#528

Will enable F in S for SOC and SDCone.

Consequently we get RootDet and GeoMean for free

update for MOI 0.9

Follow the steps of: jump-dev/SCS.jl#144

Test against CBLIB

Update README JuMP sintax

Protect against invalid params

Just add an “else” in parse arg function

Test against SDPLIB

Change icon to the new one

Note that the icon of the tab is the old one

Power cone

randsdp test fails due Inf values of variables not captured prior to primal step

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

Time limit

SDP set

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.

Failure with MOI master

See https://github.com/blegat/SolverTests/runs/1934867168

Test failures from Convex.jl tests

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.

Add Warm-start

Exclude relentr

ProxSDP should exclude relentr test:
https://travis-ci.com/blegat/SolverTests/jobs/281969930

Run Arnoldi in case Lanczos fails

Lanczos is specialized for symmetric matrices, however, it might be possible that Arnoldi, although less efficient might be more robust.

Add to JuMP and SOS docs

In SOS: https://sums-of-squares.github.io/sos/index.html
In JuMP: http://www.juliaopt.org/JuMP.jl/dev/installation/

Chordal decomposition

Example fails

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

Add SparseArrays to deps in the toml's

Variance in solution time

The issue

~~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

Code used

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?