Improve speed SDP variable entries continuous 0-1 in Convex.jl

Kieran_Marray · February 6, 2024, 12:05pm

Hi

I am trying to solve a large SDP with two matrix constraints using SCS in Convex.jl. I want to constrain all entries of my variable (a matrix) to all be continuous between 0 and 1, and was wondering what the most efficient way to set up such a problem is. Currently, setup is very slow, and the problem only works for N between 6000 and 7000 before giving memory errors on a large server.

Here is a minimum working example with small N

using Convex, SCS
const N = 100
const c = 0.3
B = ones(Float64, (N,N))
A = ones(Float64, (N,N)) # example of more complicated matrix that I want to use
X = Variable((N,N), Positive())
problem=minimize(sum(X), [eigmax(A .* X) < c, X <= B])
solve!(problem, SCS.Optimizer)```

odow · February 6, 2024, 7:43pm

When N = 6000, the X matrix is dense with 36,000,000 elements. This is very large. Even if you can build such a problem, SCS will likely have a hard time solving it.

How large N do you hope to solve?

using Convex, SCS

function solve_convex(N, c)
    A = ones(Float64, (N, N))
    X = Variable((N, N), Positive())
    problem = minimize(sum(X), [eigmax(A .* X) <= c, X <= 1])
    # Just so we can measure the time to construct in Convex
    solver = Convex.MOI.OptimizerWithAttributes(SCS.Optimizer, "max_iters" => 1)
    solve!(problem, solver)
    return X.value
end

Here are some examples:

julia> @time solve_convex(800, 0.3);
------------------------------------------------------------------
	       SCS v3.2.4 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 640002, constraints m: 1920002
cones: 	  z: primal zero / dual free vars: 319601
	  l: linear vars: 1280001
	  s: psd vars: 320400, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
	  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
	  max_iters: 1, normalize: 1, rho_x: 1.00e-06
	  acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct-amd-qdldl
	  nnz(A): 2880402, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 6.04e+00  1.04e+02  1.18e+06 -5.91e+05  1.00e-01  3.46e+00 
     1| 6.04e+00  1.04e+02  1.18e+06 -5.91e+05  1.00e-01  3.57e+00 
------------------------------------------------------------------
status:  solved (inaccurate - reached max_iters)
timings: total: 3.59e+00s = setup: 2.56e+00s + solve: 1.04e+00s
	 lin-sys: 5.67e-02s, cones: 4.48e-01s, accel: 3.70e-08s
------------------------------------------------------------------
objective = -591280.501761 (inaccurate)
------------------------------------------------------------------
┌ Warning: Problem status ITERATION_LIMIT; solution may be inaccurate.
└ @ Convex ~/.julia/packages/Convex/b2S4H/src/solution.jl:342
  5.908103 seconds (13.44 M allocations: 1.650 GiB, 6.93% gc time)

julia> @time solve_convex(1000, 0.3);
------------------------------------------------------------------
	       SCS v3.2.4 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 1000002, constraints m: 3000002
cones: 	  z: primal zero / dual free vars: 499501
	  l: linear vars: 2000001
	  s: psd vars: 500500, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
	  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
	  max_iters: 1, normalize: 1, rho_x: 1.00e-06
	  acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct-amd-qdldl
	  nnz(A): 4500502, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 6.09e+00  9.42e+01  1.34e+06 -6.70e+05  1.00e-01  5.43e+00 
     1| 6.09e+00  9.42e+01  1.34e+06 -6.70e+05  1.00e-01  5.61e+00 
------------------------------------------------------------------
status:  solved (inaccurate - reached max_iters)
timings: total: 5.64e+00s = setup: 4.03e+00s + solve: 1.61e+00s
	 lin-sys: 7.61e-02s, cones: 8.40e-01s, accel: 3.20e-08s
------------------------------------------------------------------
objective = -670302.744957 (inaccurate)
------------------------------------------------------------------
┌ Warning: Problem status ITERATION_LIMIT; solution may be inaccurate.
└ @ Convex ~/.julia/packages/Convex/b2S4H/src/solution.jl:342
 10.896022 seconds (21.00 M allocations: 2.439 GiB, 21.78% gc time)

julia> @time solve_convex(2000, 0.3);
------------------------------------------------------------------
	       SCS v3.2.4 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 4000002, constraints m: 12000002
cones: 	  z: primal zero / dual free vars: 1999001
	  l: linear vars: 8000001
	  s: psd vars: 2001000, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
	  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
	  max_iters: 1, normalize: 1, rho_x: 1.00e-06
	  acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct-amd-qdldl
	  nnz(A): 18001002, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 6.81e+00  5.74e+01  1.63e+06 -8.16e+05  1.00e-01  3.10e+01 
     1| 6.81e+00  5.74e+01  1.63e+06 -8.16e+05  1.00e-01  3.18e+01 
------------------------------------------------------------------
status:  solved (inaccurate - reached max_iters)
timings: total: 3.22e+01s = setup: 2.09e+01s + solve: 1.13e+01s
	 lin-sys: 3.86e-01s, cones: 6.47e+00s, accel: 3.40e-08s
------------------------------------------------------------------
objective = -815705.625274 (inaccurate)
------------------------------------------------------------------
┌ Warning: Problem status ITERATION_LIMIT; solution may be inaccurate.
└ @ Convex ~/.julia/packages/Convex/b2S4H/src/solution.jl:342
 51.262972 seconds (84.00 M allocations: 9.390 GiB, 14.33% gc time)

The problem gets very large very quickly.

ericphanson · February 6, 2024, 9:31pm

By the way, Positive() doesn’t create a PSD constraint, just an element-wise one. So X isn’t necessarily positive semi-definite here.

Kieran_Marray · February 6, 2024, 9:40pm

Point taken. Ignore the large N thing. What is the most efficient way to implement the constraint that each element of X must be less than 1?

odow · February 6, 2024, 9:44pm

What is the most efficient way to implement the constraint that each element of X must be less than 1?

X <= 1