Error in sparsity pattern analysis

[sebastiendesignolle]

When importing a file written in the SDPA sparse format, I get the following error with Julia 1.9.3, MOI 1.20.0, JuMP1.14.1, and COSMO 0.8.8. You can file a minimal example here; make sure to change the extension to ensure proper format recognition (GitHub does not support .dat-s files).

julia> using JuMP, COSMO

julia> model = read_from_file("minimal.dat-s")
A JuMP Model
Minimization problem with:
Variable: 1
Objective function type: AffExpr
`Vector{AffExpr}`-in-`MathOptInterface.PositiveSemidefiniteConeTriangle`: 17 constraints
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.

julia> set_optimizer(model, COSMO.Optimizer)

julia> optimize!(model)
ERROR: Input matrix is not upper triangular or has an empty column
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] QDLDL.QDLDLWorkspace(triuA::SparseArrays.SparseMatrixCSC{Float64, Int64}, AtoPAPt::Vector{Int64}, Dsigns::Nothing, regularize_eps::Float64, regularize_delta::Float64)
    @ QDLDL ~/.julia/packages/QDLDL/i77pq/src/QDLDL.jl:75
  [3] qdldl(A::SparseArrays.SparseMatrixCSC{Float64, Int64}; amd_dense_scale::Float64, perm::Vector{Int64}, logical::Bool, Dsigns::Nothing, regularize_eps::Float64, regularize_delta::Float64)
    @ QDLDL ~/.julia/packages/QDLDL/i77pq/src/QDLDL.jl:179
  [4] find_graph!(ci::COSMO.ChordalInfo{Float64}, rows::Vector{Int64}, N::Int64, C::COSMO.PsdConeTriangle{Float64})
    @ COSMO ~/.julia/packages/COSMO/WSWfU/src/chordal_decomposition/trees.jl:636
  [5] _analyse_sparsity_pattern(ci::COSMO.ChordalInfo{Float64}, csp::Vector{Int64}, sets::Vector{COSMO.AbstractConvexSet}, C::COSMO.PsdConeTriangle{Float64}, k::Int64, psd_row_range::UnitRange{Int64}, sp_ind::Int64, merge_strategy::COSMO.OptionsFactory{COSMO.CliqueGraphMerge})
    @ COSMO ~/.julia/packages/COSMO/WSWfU/src/chordal_decomposition/chordal_decomposition.jl:63
  [6] analyse_sparsity_pattern!(ci::COSMO.ChordalInfo{Float64}, csp::Vector{Int64}, sets::Vector{COSMO.AbstractConvexSet}, C::COSMO.PsdConeTriangle{Float64}, k::Int64, psd_row_range::UnitRange{Int64}, sp_ind::Int64, merge_strategy::COSMO.OptionsFactory{COSMO.CliqueGraphMerge})
    @ COSMO ~/.julia/packages/COSMO/WSWfU/src/chordal_decomposition/chordal_decomposition.jl:55
  [7] find_sparsity_patterns!(ws::COSMO.Workspace{Float64})
    @ COSMO ~/.julia/packages/COSMO/WSWfU/src/chordal_decomposition/chordal_decomposition.jl:47
  [8] chordal_decomposition!(ws::COSMO.Workspace{Float64})
    @ COSMO ~/.julia/packages/COSMO/WSWfU/src/chordal_decomposition/chordal_decomposition.jl:18
  [9] macro expansion
    @ ./timing.jl:393 [inlined]
 [10] optimize!(ws::COSMO.Workspace{Float64})
    @ COSMO ~/.julia/packages/COSMO/WSWfU/src/solver.jl:90
 [11] optimize!
    @ ~/.julia/packages/COSMO/WSWfU/src/MOI_wrapper.jl:151 [inlined]
 [12] optimize!
    @ ~/.julia/packages/MathOptInterface/LQvlf/src/MathOptInterface.jl:85 [inlined]
 [13] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{COSMO.Optimizer{Float64}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/LQvlf/src/Utilities/cachingoptimizer.jl:316
 [14] optimize!
    @ ~/.julia/packages/MathOptInterface/LQvlf/src/Bridges/bridge_optimizer.jl:376 [inlined]
 [15] optimize!
    @ ~/.julia/packages/MathOptInterface/LQvlf/src/MathOptInterface.jl:85 [inlined]
 [16] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{COSMO.Optimizer{Float64}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/LQvlf/src/Utilities/cachingoptimizer.jl:316
 [17] optimize!(model::Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ JuMP ~/.julia/packages/JuMP/mvUVO/src/optimizer_interface.jl:439
 [18] optimize!(model::Model)
    @ JuMP ~/.julia/packages/JuMP/mvUVO/src/optimizer_interface.jl:409
 [19] top-level scope
[minimal.txt](https://github.com/oxfordcontrol/COSMO.jl/files/12583601/minimal.txt)

    @ REPL[3]:1

Note that I do not get this error with other solvers (I tried Mosek and SCS), neither with the toy example provided in the page explaining the format, so it may have a different origin than the format itself.

[sebastiendesignolle]

Ok, the error message is actually quite clear: my problem has empty columns. This is detected by other solvers but seems to not be preprocessed by COSMO. I'll adapt my input files accordingly, but would it be a potential feature to add?

[goulart-paul]

I think it is still a bug, because the solver is written such that it is (or at least should be...) robust to empty columns in the user input. The solver adds explicit small shifts to the diagonal of the KKT matrix specifically to avoid this condition being reached.

I suggest to leave this issue open until we can test the example you provided.

[sebastiendesignolle]

Indeed, even when I remove the empty columns, I still obtain an error. Actually not with the example above, but here is a minimal example in the dat-s format.

1
1
3
1
0 1 1 1 1 
0 1 1 2 1 
0 1 2 2 1 
1 1 1 1 1 
1 1 1 2 1 
1 1 1 3 1

Running the commands from my first message gives the exact same error, although no column is empty for the matrix F1 (see the details of the format here). Maybe the current sparsity pattern recognition only uses F0, which indeed has an empty column here?

[goulart-paul]

It appear to be a bug in the preprocessing stage that does chordal decomposition. This works:

julia> model = read_from_file("minimal.dat-s");

julia> set_optimizer(model, COSMO.Optimizer)

julia> set_optimizer_attribute(model,"decompose",false)

julia> optimize!(model)
------------------------------------------------------------------
          COSMO v0.8.7 - A Quadratic Objective Conic Solver
                         Michael Garstka
                University of Oxford, 2017 - 2022
------------------------------------------------------------------

Problem:  x ∈ R^{1},
          constraints: A ∈ R^{282x1} (7 nnz),
          matrix size to factor: 283x283,
          Floating-point precision: Float64
Sets:     PsdConeTriangle of dim: 78 (12x12)
          PsdConeTriangle of dim: 45 (9x9)
          PsdConeTriangle of dim: 45 (9x9)
          PsdConeTriangle of dim: 21 (6x6)
          PsdConeTriangle of dim: 21 (6x6)
          ... and 12 more
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.63ms

Iter:	Objective:	Primal Res:	Dual Res:	Rho:
1	-8.1491e-02	6.9120e+02	6.4174e-01	1.0000e-01
25	 4.5534e-20	2.8422e-13	3.8459e-16	1.0000e-01

------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 26 (incl. 1 safeguarding iter)
Optimal objective: 4.553e-20
Runtime: 0.004s (4.37ms)

I think @migarstka is probably much more likely to understand how to fix that though.

[migarstka]

Thanks for raising this. This is known behaviour.

COSMO expects PSD cones to have connected sparsity patterns (no zero rows/columns). If you have zero rows, this means you can split the PSD constraint into two decomposed constraints.

I agree it would be preferable to handle this automatically, so this could be a potential feature request. It would be better though if that happened in JuMP/MOI. I am surprised that this does not happen automatically, but I am assuming read_from_file just follows the strict SDPA format definition when the file content is translated into a MOI model.

[sebastiendesignolle]

Hello, thanks for your replies.

Deactivating the chordal decomposition indeed removes the bug, but I suspect some users (not me in my case actually) may want to take advantage of this possible acceleration.

I do agree that removing empty columns belong either to the user (writing the SDPA file) or to MOI. However, as far as I can see, the minimal example from my last message does not contain any empty column as far as the global sparsity pattern is concerned. The matrices look like:

F0
1 1 0
1 1 0
0 0 0

F1
1 1 1
1 0 0
1 0 0

The sum of absolute values of F0 and F1 has a sparsity pattern with no empty row/column, but the decomposition raises the same error.

[goulart-paul]

I think the problem there is that the union of the sparsity patterns does not contain the diagonal. Perhaps when analysing chordal sparsity internally it makes sense to just include an identity matrix into this union when computing the decomposition?