Sparse constraints with COSMO.jl

Specific Domains Optimization (Mathematical)
1

I am looking to solve a very large but sparse optimization problem using COSMO.jl. When creating my constraint set

COSMO.Constraint(sparse(stack(A_i_vecs)), sparsevec(p.B), COSMO.PsdCone)

where A = sparse(stack(A_i_vecs)) is the matrix, and b = sparsevec(p.B) is the matrix constraint of the form Ax + b \in \mathcal{K} as in Getting Started · COSMO.jl

However, when I attempt to apply their decomposition techniques, it fails to create a sparse PsdCone and instead uses densePsdCone. I can supply more code as needed, but currently my questions are: what conditions are needed on A and b to allow the sparse PsdCone to be used (my objective and constraint above are sparse), and how can I find out what sparsity cliques were made from COSMO? The printout (even with Verbose on) does not give me any information.

2

It’s hard to say anything without a reproducible example. Is your matrix symmetric etc?

You could also consider using the JuMP interface. The COSMO documentation has a relevant section:

https://oxfordcontrol.github.io/COSMO.jl/stable/decomposition/

3

When using JuMP, the solver still does not create the sparse structure that I would expect.

Creating a minimal example might take a little bit, so I will try to make something tomorrow. I think I need to do a little more to isolate the problem.

For now, I do know that the matrix A is not symmetric, and that it is less than 1% nonzero. So, COSMO should be correctly using the sparse data structures available at a minimum. From the COSMO documentation, it seems like the decomposition techniques require symmetry… Chordal Decomposition · COSMO.jl . Perhaps I should dig a little deeper before making a minimal example (there’s a lot of code going on).

If I still think this is a cosmo / jump issue and not a problem formulation issue, I will make a comment here (should I do an @ in case it takes a little while for me to create it?).

When uploading a minimal example, I think uploading the matrices I am using is the easiest way to make the example work. Besides google drive (first thing that comes to mind), is there a good way to do this?

Thanks for the help!

4

No need to @. I tend to read every post :smile:

This example (from Chordal Decomposition · COSMO.jl) may help so you the comparison between JuMP and COSMO:

julia> using COSMO, JuMP, LinearAlgebra, SparseArrays

julia> begin
           A1 = [-4.0 0.0 -2.0 0.0 0.0 -1.0 0.0 0.0 0.0; 0.0 -3.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0; -2.0 -1.0 -2.0 0.0 0.0 5.0 4.0 -4.0 0.0; 0.0 0.0 0.0 -4.0 -5.0 0.0 0.0 3.0 0.0; 0.0 0.0 0.0 -5.0 4.0 0.0 0.0 2.0 0.0; -1.0 0.0 5.0 0.0 0.0 5.0 -4.0 -4.0 -5.0; 0.0 0.0 4.0 0.0 0.0 -4.0 -1.0 -1.0 -3.0; 0.0 0.0 -4.0 3.0 2.0 -4.0 -1.0 2.0 -2.0; 0.0 0.0 0.0 0.0 0.0 -5.0 -3.0 -2.0 -3.0];
           A2 = [-5.0 0.0 3.0 0.0 0.0 -2.0 0.0 0.0 0.0; 0.0 -3.0 -5.0 0.0 0.0 0.0 0.0 0.0 0.0; 3.0 -5.0 3.0 0.0 0.0 5.0 -4.0 -5.0 0.0; 0.0 0.0 0.0 3.0 2.0 0.0 0.0 -2.0 0.0; 0.0 0.0 0.0 2.0 4.0 0.0 0.0 -3.0 0.0; -2.0 0.0 5.0 0.0 0.0 1.0 -5.0 -2.0 -4.0; 0.0 0.0 -4.0 0.0 0.0 -5.0 -2.0 -3.0 3.0; 0.0 0.0 -5.0 -2.0 -3.0 -2.0 -3.0 5.0 3.0; 0.0 0.0 0.0 0.0 0.0 -4.0 3.0 3.0 -4.0];
           B = [-0.11477375644968069 0.0 6.739182490600791 0.0 0.0 -1.2185593245043502 0.0 0.0 0.0; 0.0 1.2827680528587497 -5.136452036888789 0.0 0.0 0.0 0.0 0.0 0.0; 6.739182490600791 -5.136452036888789 7.344770673489607 0.0 0.0 -0.2224400187044442 -10.505300166831221 -1.2627361794562273 0.0; 0.0 0.0 0.0 10.327710040060499 8.91534585379813 0.0 0.0 -6.525873789637007 0.0; 0.0 0.0 0.0 8.91534585379813 0.8370459338528677 0.0 0.0 -6.210900615408826 0.0; -1.2185593245043502 0.0 -0.2224400187044442 0.0 0.0 -3.8185953011245024 -0.994033914192722 2.8156077981712997 1.4524716674219218; 0.0 0.0 -10.505300166831221 0.0 0.0 -0.994033914192722 0.029162208619863517 -2.8123790276830745 7.663416446183705; 0.0 0.0 -1.2627361794562273 -6.525873789637007 -6.210900615408826 2.8156077981712997 -2.8123790276830745 4.71893305728242 6.322431630550857; 0.0 0.0 0.0 0.0 0.0 1.4524716674219218 7.663416446183705 6.322431630550857 0.5026094532322212];
           c = [-0.21052661285686525, -1.263324575834677];
       end;

julia> begin
           model = JuMP.Model(COSMO.Optimizer)
           @variable(model, x[1:2]);
           @objective(model, Min, c' * x )
           @constraint(model, Symmetric(B - A1  .* x[1] - A2 .* x[2])  in JuMP.PSDCone());
           JuMP.optimize!(model)
       end
------------------------------------------------------------------
          COSMO v0.8.8 - A Quadratic Objective Conic Solver
                         Michael Garstka
                University of Oxford, 2017 - 2022
------------------------------------------------------------------

Problem:  x ∈ R^{7},
          constraints: A ∈ R^{30x7} (58 nnz),
          matrix size to factor: 37x37,
          Floating-point precision: Float64
Sets:     PsdConeTriangle of dim: 15 (5x5)
          PsdConeTriangle of dim: 6 (3x3)
          PsdConeTriangle of dim: 6 (3x3)
          PsdConeTriangle of dim: 3 (2x2)
Decomp:   Num of original PSD cones: 1
          Num decomposable PSD cones: 1
          Num PSD cones after decomposition: 4
          Merge Strategy: CliqueGraphMerge
Settings: ϵ_abs = 1.0e-05, ϵ_rel = 1.0e-05,
          ϵ_prim_inf = 1.0e-04, ϵ_dual_inf = 1.0e-04,
          ρ = 0.1, σ = 1e-06, α = 1.6,
          max_iter = 5000,
          scaling iter = 10 (on),
          check termination every 25 iter,
          check infeasibility every 40 iter,
          KKT system solver: QDLDL
Acc:      Anderson Type2{QRDecomp},
          Memory size = 15, RestartedMemory,	
          Safeguarded: true, tol: 2.0
Setup Time: 0.1ms

Iter:	Objective:	Primal Res:	Dual Res:	Rho:
1	-2.2362e+00	8.9140e+00	1.2209e+00	1.0000e-01
25	-1.4134e+00	1.2134e-06	1.2496e-06	1.0000e-01

------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 25
Optimal objective: -1.413
Runtime: 0.003s (3.4ms)


julia> begin
           A = -sparse(hcat(vec(A1), vec(A2)))
           b = sparse(vec(B))
           c_model = COSMO.Model()
           constraints = [COSMO.Constraint(A, b, COSMO.PsdCone)]
           # Seems like `decompose = false` needed here. Not sure if intended or a bug
           settings = COSMO.Settings(verbose = true, decompose = false)
           COSMO.assemble!(c_model, zeros(2, 2), c, constraints; settings)
           results = COSMO.optimize!(c_model)
       end
------------------------------------------------------------------
          COSMO v0.8.8 - A Quadratic Objective Conic Solver
                         Michael Garstka
                University of Oxford, 2017 - 2022
------------------------------------------------------------------

Problem:  x ∈ R^{2},
          constraints: A ∈ R^{81x2} (78 nnz),
          matrix size to factor: 83x83,
          Floating-point precision: Float64
Sets:     PsdCone of dim: 81 (9x9)
Settings: ϵ_abs = 1.0e-05, ϵ_rel = 1.0e-05,
          ϵ_prim_inf = 1.0e-04, ϵ_dual_inf = 1.0e-04,
          ρ = 0.1, σ = 1e-06, α = 1.6,
          max_iter = 5000,
          scaling iter = 10 (on),
          check termination every 25 iter,
          check infeasibility every 40 iter,
          KKT system solver: QDLDL
Acc:      Anderson Type2{QRDecomp},
          Memory size = 15, RestartedMemory,	
          Safeguarded: true, tol: 2.0
Setup Time: 0.1ms

Iter:	Objective:	Primal Res:	Dual Res:	Rho:
1	-2.2176e+00	6.3032e+00	1.2633e+00	1.0000e-01
25	-1.4134e+00	1.6825e-11	1.3863e-11	1.0000e-01

------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 25
Optimal objective: -1.413
Runtime: 0.002s (2.16ms)

>>> COSMO - Results
Status: Solved
Iterations: 25 (incl. 0 safeguarding iterations)
Optimal Objective: -1.413
Runtime: 2.16ms
Setup Time: 0.1ms

Avg Iter Time: 0.04ms

julia> begin
           # If using Triangle, you need to scale the off-diagonal by √2
           scale(i, j) = i == j ? 1.0 : sqrt(2)
           upper_tri(A) = [scale(i, j) * A[i, j] for j in 1:size(A, 2) for i in 1:j]
           A = -sparse(hcat(upper_tri(A1), upper_tri(A2)))
           b = sparse(upper_tri(B))
           constraints = [COSMO.Constraint(A, b, COSMO.PsdConeTriangle)]
           c_model = COSMO.Model()
           settings = COSMO.Settings(verbose = true)
           COSMO.assemble!(c_model, zeros(2, 2), c, constraints; settings)
           results = COSMO.optimize!(c_model)
       end
------------------------------------------------------------------
          COSMO v0.8.8 - A Quadratic Objective Conic Solver
                         Michael Garstka
                University of Oxford, 2017 - 2022
------------------------------------------------------------------

Problem:  x ∈ R^{7},
          constraints: A ∈ R^{30x7} (58 nnz),
          matrix size to factor: 37x37,
          Floating-point precision: Float64
Sets:     PsdConeTriangle of dim: 15 (5x5)
          PsdConeTriangle of dim: 6 (3x3)
          PsdConeTriangle of dim: 6 (3x3)
          PsdConeTriangle of dim: 3 (2x2)
Decomp:   Num of original PSD cones: 1
          Num decomposable PSD cones: 1
          Num PSD cones after decomposition: 4
          Merge Strategy: CliqueGraphMerge
Settings: ϵ_abs = 1.0e-05, ϵ_rel = 1.0e-05,
          ϵ_prim_inf = 1.0e-04, ϵ_dual_inf = 1.0e-04,
          ρ = 0.1, σ = 1e-06, α = 1.6,
          max_iter = 5000,
          scaling iter = 10 (on),
          check termination every 25 iter,
          check infeasibility every 40 iter,
          KKT system solver: QDLDL
Acc:      Anderson Type2{QRDecomp},
          Memory size = 15, RestartedMemory,	
          Safeguarded: true, tol: 2.0
Setup Time: 0.2ms

Iter:	Objective:	Primal Res:	Dual Res:	Rho:
1	-2.2362e+00	8.9140e+00	1.2209e+00	1.0000e-01
25	-1.4134e+00	1.2134e-06	1.2496e-06	1.0000e-01

------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 25
Optimal objective: -1.413
Runtime: 0.004s (3.67ms)

>>> COSMO - Results
Status: Solved
Iterations: 25 (incl. 0 safeguarding iterations)
Optimal Objective: -1.413
Runtime: 3.67ms
Setup Time: 0.2ms

Avg Iter Time: 0.03ms