High Memory Usage with Convex.jl When Solving Optimization Problem with Sparse Matrices

[CaG21]

Hi everyone,

I am using Convex.jl to solve an optimization problem. Even though it is relatively large, its constraints can be modeled using sparse matrices, which occupy a small amount of memory (about 200 MB). However, as soon as I run the solve! command, the memory usage increases significantly, quickly saturating the RAM. It appears as though the Convex.jl package is converting the sparse matrices into dense matrices.

I have prepared a MWE to illustrate the issue. In this example, the memory usage increases from few tens of MB to approximately 5 GB.

Here’s the code for the MWE:

using Convex
const MOI = Convex.MOI
using ECOS
using SparseArrays

x = Variable(300, 101);

A = sprandn(30000, 30300, 0.01);
B = sprandn(50000, 30300, 0.01);

model = minimize(-x[end, end]);

model.constraints += A * vec(x) == 0;
model.constraints += B * vec(x) <= 0;

opt = MOI.OptimizerWithAttributes(ECOS.Optimizer, MOI.Silent() => true);

solve!(model, opt)

Could someone help me understand why this is happening and how to avoid this problem?

Thanks in advance for your help!

[odow]

I assume somewhere Convex is converting your matrices to dense.

I’ll take a look (Sparse matrices in A * x · Issue #697 · jump-dev/Convex.jl · GitHub), but before then, you should use JuMP:

using ECOS
using JuMP
using SparseArrays
A = sprandn(30_000, 30_300, 0.01);
B = sprandn(50_000, 30_300, 0.01);
model = Model(ECOS.Optimizer)
set_silent(model)
@variable(model, x[1:300, 1:101])
@objective(model, Min, -x[end, end])
@constraint(model, A * vec(x) == 0)
@constraint(model, B * vec(x) <= 0)
optimize!(model)

[odow]

Are you using the latest version of Convex.jl?

I get this:

julia> using Convex

julia> import SCS

julia> import SparseArrays

julia> function main(m, n)
           A = SparseArrays.sprandn(m * (n - 1), m * n, 0.01);
           B = SparseArrays.sprandn(m * (n - 1), m * n, 0.01);
           @show SparseArrays.nnz(A) + SparseArrays.nnz(B)
           x = Variable(m, n)
           model = minimize(-x[end, end], [A * vec(x) == 0, B * vec(x) <= 0])
           solve!(model, SCS.Optimizer)
           return model
       end
main (generic function with 1 method)

julia> @time main(30, 101)
SparseArrays.nnz(A) + SparseArrays.nnz(B) = 181890
[ Info: [Convex.jl] Compilation finished: 0.11 seconds, 21.368 MiB of memory allocated
------------------------------------------------------------------
	       SCS v3.2.4 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 3030, constraints m: 6000
cones: 	  z: primal zero / dual free vars: 3000
	  l: linear vars: 3000
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: 100000, normalize: 1, rho_x: 1.00e-06
	  acceleration_lookback: 10, acceleration_interval: 10
	  compiled with openmp parallelization enabled
lin-sys:  sparse-direct-amd-qdldl
	  nnz(A): 181890, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 3.55e-01  9.37e-01  1.71e-01 -8.56e-02  1.00e-01  5.11e+00 
    25| 1.90e-07  2.48e-08  3.77e-08 -1.88e-08  1.00e-01  5.40e+00 
------------------------------------------------------------------
status:  solved
timings: total: 5.40e+00s = setup: 5.08e+00s + solve: 3.21e-01s
	 lin-sys: 3.06e-01s, cones: 4.66e-04s, accel: 4.63e-04s
------------------------------------------------------------------
objective = -0.000000
------------------------------------------------------------------
  5.571585 seconds (4.61 k allocations: 39.316 MiB)
Problem statistics
  problem is DCP         : true
  number of variables    : 1 (3_030 scalar elements)
  number of constraints  : 2 (6_000 scalar elements)
  number of coefficients : 187_892
  number of atoms        : 12

Solution summary
  termination status : OPTIMAL
  primal status      : FEASIBLE_POINT
  dual status        : FEASIBLE_POINT
  objective value    : -0.0

Expression graph
  minimize
   └─ Convex.NegateAtom (affine; real)
      └─ index (affine; real)
         └─ 30×101 real variable (id: 563…038)
  subject to
   ├─ == constraint (affine)
   │  └─ + (affine; real)
   │     ├─ * (affine; real)
   │     │  ├─ …
   │     │  └─ …
   │     └─ Convex.NegateAtom (constant; negative)
   │        └─ …
   └─ ≤ constraint (affine)
      └─ + (affine; real)
         ├─ * (affine; real)
         │  ├─ …
         │  └─ …
         └─ Convex.NegateAtom (constant; negative)
            └─ …


julia> @time main(300, 101)
SparseArrays.nnz(A) + SparseArrays.nnz(B) = 18172888
[ Info: [Convex.jl] Compilation finished: 3.96 seconds, 1.998 GiB of memory allocated
------------------------------------------------------------------
	       SCS v3.2.4 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 30300, constraints m: 60000
cones: 	  z: primal zero / dual free vars: 30000
	  l: linear vars: 30000
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: 100000, normalize: 1, rho_x: 1.00e-06
	  acceleration_lookback: 10, acceleration_interval: 10
	  compiled with openmp parallelization enabled
lin-sys:  sparse-direct-amd-qdldl
	  nnz(A): 18172888, nnz(P): 0
^C

where Convex compiles the problem in 4 seconds, but then SCS gets stuck trying to factor the matrix.

Is this your real problem? The problem is difficult to solve, even by a specialized linear optimizer like Gurobi. 0.01 seems small, but the problem is still quite “dense”.

[CaG21]

I tried to use JuMP. The issue seems to be mitigated for the MWE, with 3 GB of RAM used, instead of 5.

[odow]

What is the output of import Pkg; Pkg.status()? Did you check if you are using the latest version of Convex? @ericphanson fixed some stuff with sparse matrices recently.

[CaG21]

I am actually using the version 0.15.4. I’ll try to update to the latest version.

The MWE is just an example to reproduce the issue I thought I’ve encountered, i.e., a big sparse problem seems to use more memory than expected, but it is not fully represented of my original optimization problem. In my orginal problem, I have also some second-order cone contraints (that I did not report in the MWE, since, in the original problem, they seem to have a minimal impact from the memory point of view), and matrices are bigger but sparser (<0.001).

[odow]

I am actually using the version 0.15.4

Yes, v0.16 should be much faster. See:

ANN: Convex.jl v0.16 is released

Dear all, We are pleased to announce the release of Convex.jl v0.16.0. This release of Convex involved a substantial internal rewrite of Convex.jl to integrate better with MathOptInterface. The result is a faster, simpler, and more maintainable package that is nearly ready for a 1.0 release. Detailed release notes are available at: Release notes · Convex.jl Most of the very large number of changes (131 commits, 270 files changed) affect only the internal representation of Convex.jl models, a…

[CaG21]

I updated to v0.16. It seems faster and the issue is a lot mitigated. It has a spike in RAM usage at the beginning, but later on it sets to the expected value of memory usage.

OT: Even though the update was mostly painfree, I was forced to change part of my code since I got the error distance_to_set using the distance metric MathOptInterface.Utilities.ProjectionUpperBoundDistance() for set type MathOptInterface.Nonpositives has not been implemented yet.
This error is generated by the constraint
push!(model.constraints, sum(γ, dims=1) <= 1)
only when I fix γ (with γ = Variable(n, m, BinVar)).
I solved this issue by rewriting the constraint with a for cycle

for i in 1:m
    push!(model.constraints, sum(γ[:, i]) <= 1)
end

I do not know if it is worth opening a post on this.

[odow]

This error should only occur if γ is a constant, and if it is, why are you adding a constraint to it?

I have a reproducible example, so I’ll open an issue:

julia> using Convex, HiGHS

julia> x = Variable(2, 2)
Variable
size: (2, 2)
sign: real
vexity: affine
id: 909…059

julia> fix!(x, [1 0; 0 1])
Variable
size: (2, 2)
sign: real
vexity: constant
id: 909…059
value: [1 0; 0 1]

julia> problem = satisfy([sum(x; dims = 1) <= 1])
Problem statistics
  problem is DCP         : true
  number of variables    : 1 (4 scalar elements)
  number of constraints  : 1 (2 scalar elements)
  number of coefficients : 5
  number of atoms        : 4

Solution summary
  termination status : OPTIMIZE_NOT_CALLED
  primal status      : NO_SOLUTION
  dual status        : NO_SOLUTION

Expression graph
  satisfy
   └─ nothing
  subject to
   └─ ≤ constraint (constant)
      └─ + (constant; real)
         ├─ * (constant; real)
         │  ├─ …
         │  └─ …
         └─ Convex.NegateAtom (constant; negative)
            └─ …


julia> solve!(problem, HiGHS.Optimizer)
ERROR: distance_to_set using the distance metric MathOptInterface.Utilities.ProjectionUpperBoundDistance() for set type MathOptInterface.Nonpositives has not been implemented yet.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] distance_to_set(d::MathOptInterface.Utilities.ProjectionUpperBoundDistance, ::Matrix{Float64}, set::MathOptInterface.Nonpositives)
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/gBojA/src/Utilities/distance_to_set.jl:75

[odow]

This was a bug in Convex so thanks for reporting it

I have a fix: Fix distance_to_set for sets without a definition in MOI by odow · Pull Request #699 · jump-dev/Convex.jl · GitHub