Chordal decomposition
This tutorial was generated using Literate.jl. Download the source as a .jl file.
The purpose of this tutorial is to show how to use MathOptChordalDecomposition.jl to improve the performance of models with PSD constraints.
Required packages
This tutorial uses the following packages:
using JuMP
import MathOptChordalDecomposition
import SCS
import SparseArraysBackground
Chordal decomposition is a technique for decomposing a large PSD constraint into a set of smaller PSD constraints and some linear equality constraints.
If the original PSD constraint is sparse, the decomposed problem can be faster to solve than the original.
For more information on chordal decomposition, watch Michael Garstka's talk at JuMP-dev 2019.
Some solvers, such as Clarabel.jl and COSMO.jl implement chordal decomposition internally. Others, such as SCS.jl do not implement chordal decomposition.
The Julia package MathOptChordalDecomposition.jl is a MathOptInterface layer that implements chordal decomposition of sparse semidefinite constraints. It can be used to wrap any solver which supports PSD constraints and does not implement chordal decomposition internally.
JuMP Model
To demonstrate the benefits of chordal decomposition, we use the mcp124-1 model from SDPLIB.
model = read_from_file(joinpath(@__DIR__, "mcp124-1.dat-s"))A JuMP Model
├ solver: none
├ objective_sense: MIN_SENSE
│ └ objective_function_type: AffExpr
├ num_variables: 124
├ num_constraints: 1
│ └ Vector{AffExpr} in MOI.PositiveSemidefiniteConeTriangle: 1
└ Names registered in the model: noneThis model has 124 decision variables and one PSD constraint. This PSD constraint is sparse, which means that many elements of the matrix are zero.
To view the matrix, use all_constraints to get a list of the constraints, then use constraint_object to get the function and set form of the constraint:
S = MOI.PositiveSemidefiniteConeTriangle
constraints = all_constraints(model, Vector{AffExpr}, S)
con = constraint_object(constraints[1]);
con.setMathOptInterface.PositiveSemidefiniteConeTriangle(124)con.func7750-element Vector{AffExpr}:
_[1] - 0.25
0
_[2] - 0.5
0
0
_[3] - 0.25
0
0
0
_[4]
⋮
0
0
0
0
0
0
0
0
_[124] - 1The constraint function is given in vectorized form. Use reshape_vector to convert it into a matrix:
F = reshape_vector(con.func, SymmetricMatrixShape(con.set.side_dimension))124×124 LinearAlgebra.Symmetric{AffExpr, Matrix{AffExpr}}:
_[1] - 0.25 0 0 0 … 0 0
0 _[2] - 0.5 0 0 0 0
0 0 _[3] - 0.25 0 0 0
0 0 0 _[4] 0 0
0 0 0 0 0 0
0 0 0 0 … 0 0
0 0 0 0 0 0
0 0 0.25 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
⋮ ⋱
0 0 0 0 … 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 … 0 0
0 0 0 0 0 0
0 0 0 0 _[123] - 1 0
0 0 0 0 0 _[124] - 1The F matrix is dense, but many elements are zero. Use SparseArrays.sparse to turn it into a sparse matrix:
A = SparseArrays.sparse(F)124×124 SparseArrays.SparseMatrixCSC{AffExpr, Int64} with 422 stored entries:
⎡⠑⢄⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠀⡀⠀⠀⠀⠂⡀⢀⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⢀⠀⢀⠀⠀⠀⎤
⎢⠈⠀⡑⢌⠀⠠⠀⠀⠈⢀⠀⠠⠃⠀⡀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠐⠠⠀⠄⠀⠀⠀⠀⠀⠠⠀⠀⠀⠀⠄⠀⎥
⎢⠀⠀⠀⡀⠑⢄⠀⠀⠀⠠⠀⠀⠀⠄⠀⠠⢈⠀⠀⠀⠠⠀⠁⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠀⠠⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠐⡀⠁⠀⡂⠀⠀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⣀⠀⠀⠠⠠⠀⠑⎥
⎢⠀⠀⠂⢀⠀⡀⠀⠀⠑⢄⠀⠀⠡⠀⠂⠀⢀⠀⠀⠀⠀⠀⠂⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⢀⎥
⎢⠀⠀⠀⡀⠀⠀⠀⠀⠀⠀⠑⢄⠐⠔⠀⠀⠂⠁⠀⠀⠀⠂⠀⠀⠀⠀⣀⠀⠀⠀⠀⠆⠀⠂⠠⠤⠀⠆⠂⠄⎥
⎢⠀⠀⠉⠀⠀⠄⢀⠀⠁⠂⢐⠄⠑⢄⡈⠀⠄⠀⠀⠀⡄⠀⠀⠀⠈⠀⠀⠈⠀⡁⠀⠀⠀⠀⠀⠀⠙⠀⠀⠀⎥
⎢⠈⠀⠀⠈⠀⡀⠄⠈⠈⠀⠀⠀⠂⠈⠑⢄⠂⠀⠠⠀⠀⠀⠀⢄⠀⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠈⠈⠀⠀⡀⎥
⎢⠀⠠⠀⠀⠂⠐⠠⠠⠀⠐⠌⠀⠀⠁⠈⠀⠑⢄⠀⠀⠀⠀⡀⠀⡀⠂⢀⠀⠀⠀⠀⠂⠄⠀⠀⠢⠀⠀⡀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠀⠑⢄⠀⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⠀⠂⠀⢀⢀⠀⠀⠈⠈⎥
⎢⠠⠀⠀⠁⠀⠂⠁⠀⠀⠀⠠⠀⠀⠉⠀⠀⠀⠀⠀⠀⠑⢄⠀⠄⠀⠠⠀⠀⠀⠁⠀⡀⠀⠈⠀⠠⠀⢀⡀⠀⎥
⎢⠀⢈⠀⠀⠁⠀⠀⠀⢈⠀⠀⠀⠀⠀⠀⢄⠀⠈⠀⠀⠀⠄⠑⢄⢀⠀⡀⠀⠀⠀⠀⠀⠈⠀⠀⠐⠀⠀⠈⠀⎥
⎢⠀⠀⠐⡀⠀⠁⠀⠀⠀⠀⠀⠀⠂⠀⠀⠀⠠⠈⠂⠀⠀⡀⠀⠐⠑⢄⠁⠐⠀⠀⠀⠀⠀⠀⠀⠈⠀⠐⠀⠀⎥
⎢⡀⠀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠘⡀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠈⢁⠀⠑⢄⠀⠠⠐⠀⣀⠀⢠⠁⠀⢀⠂⢀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠄⠠⠀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠀⠀⠀⡀⢑⢔⠀⠀⠀⠀⠀⠠⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠄⠀⠀⠀⠐⠠⠀⠀⠀⠀⠠⠀⠀⠀⠀⠐⠀⠀⠀⠑⢄⠀⠀⠈⠁⠀⠀⠁⡀⎥
⎢⠀⠀⠀⡀⠀⠄⠈⢠⠀⠀⠠⠀⠀⠀⠀⠀⠀⠁⠈⠀⡀⠀⠂⠀⠀⠀⠀⠘⠀⠀⠀⠀⠑⢄⠀⠀⢀⠄⠆⠀⎥
⎢⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⡆⠀⠀⡀⠀⠠⡀⠀⢐⠀⡀⢀⠀⡀⠀⠄⠒⠀⡀⠆⠀⠀⠀⢑⣴⠀⠀⠀⠄⎥
⎢⠀⠐⠀⠀⠀⠀⠀⡂⠀⠐⠠⠄⠓⠀⠂⠀⠀⠀⠀⠀⠀⢀⠀⠀⢀⠀⠀⢀⠀⠀⠀⠀⠀⠔⠀⠀⠑⢄⠀⠀⎥
⎣⠀⠀⠀⠁⠀⠂⢄⠀⠀⢀⠈⠄⠀⠀⠀⠠⠀⠈⡂⠀⠀⠈⠂⠀⠀⠀⠈⢀⠀⠀⠁⠠⠈⠁⠀⠄⠀⠀⠑⢄⎦The sparse matrix has 422 non-zeros, which is a density of 2.7%:
SparseArrays.nnz(A) / size(A, 1)^20.027445369406867846Solution speed
SCS.jl is a first-order solver that does not exploit the sparsity of PSD constraints. Let's solve it and see how long it took:
set_optimizer(model, SCS.Optimizer)
@time optimize!(model)------------------------------------------------------------------
SCS v3.2.7 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 124, constraints m: 7750
cones: s: psd vars: 7750, 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: 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): 124, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 9.71e-01 3.87e+00 3.12e+02 2.65e+02 1.00e-01 4.38e-03
250| 1.20e-03 3.27e-04 2.64e-01 1.42e+02 2.92e-02 8.12e-01
500| 2.97e-04 7.52e-05 4.58e-02 1.42e+02 2.92e-02 1.62e+00
750| 2.40e-04 6.14e-05 3.89e-02 1.42e+02 2.92e-02 2.42e+00
1000| 1.91e-04 5.10e-05 3.30e-02 1.42e+02 2.92e-02 3.23e+00
1250| 1.55e-04 4.34e-05 2.78e-02 1.42e+02 2.92e-02 4.04e+00
1500| 1.34e-04 3.65e-05 2.34e-02 1.42e+02 2.92e-02 4.83e+00
1750| 2.19e-01 2.84e-01 3.74e+00 1.44e+02 2.92e-02 5.63e+00
2000| 1.01e-04 2.57e-05 1.67e-02 1.42e+02 2.92e-02 6.42e+00
2250| 9.41e-05 2.15e-05 1.42e-02 1.42e+02 2.92e-02 7.24e+00
------------------------------------------------------------------
status: solved
timings: total: 7.24e+00s = setup: 1.95e-03s + solve: 7.24e+00s
lin-sys: 1.64e-01s, cones: 7.00e+00s, accel: 1.36e-02s
------------------------------------------------------------------
objective = 141.972713
------------------------------------------------------------------
7.450631 seconds (152.84 k allocations: 8.637 MiB, 2.79% compilation time: 56% of which was recompilation)In comparison, if we wrap SCS.Optimizer in MathOptChordalDecomposition.Optimizer, then the problem takes less than 1 second to solve:
set_optimizer(model, () -> MathOptChordalDecomposition.Optimizer(SCS.Optimizer))
@time optimize!(model)------------------------------------------------------------------
SCS v3.2.7 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 1155, constraints m: 8781
cones: z: primal zero / dual free vars: 7750
s: psd vars: 1031, ssize: 115
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): 2186, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 2.75e+01 1.00e+00 9.93e+03 -4.81e+03 1.00e-01 3.02e-03
250| 3.65e-03 9.48e-04 1.32e-03 1.42e+02 8.07e-01 1.07e-01
500| 2.63e-04 7.03e-05 1.77e-04 1.42e+02 2.57e+00 2.13e-01
------------------------------------------------------------------
status: solved
timings: total: 2.13e-01s = setup: 2.55e-03s + solve: 2.10e-01s
lin-sys: 6.14e-02s, cones: 1.30e-01s, accel: 2.76e-03s
------------------------------------------------------------------
objective = 141.988826
------------------------------------------------------------------
0.333878 seconds (84.60 k allocations: 6.947 MiB, 34.54% compilation time: 90% of which was recompilation)The difference in performance is because of the chordal decomposition. The decomposed problem introduced new variables (there are now 1,155 variables instead of 124) and constraints (there are now 115 PSD constraints instead of one), but each PSD constraint is smaller than the original.
decom = unsafe_backend(model)MathOptChordalDecomposition.Optimizer{CliqueTrees.MF}
├ ObjectiveSense: MIN_SENSE
├ ObjectiveFunctionType: MOI.ScalarAffineFunction{Float64}
├ NumberOfVariables: 1155
└ NumberOfConstraints: 116
├ MOI.VectorAffineFunction{Float64} in MOI.Zeros: 1
└ MOI.VectorOfVariables in MOI.PositiveSemidefiniteConeTriangle: 115With a bit of effort, we can compute the number of PSD constraints of each size:
count_by_size = Dict{Int,Int}()
for ci in MOI.get(decom, MOI.ListOfConstraintIndices{MOI.VectorOfVariables,S}())
set = MOI.get(decom, MOI.ConstraintSet(), ci)
n = set.side_dimension
count_by_size[n] = get(count_by_size, n, 0) + 1
end
count_by_sizeDict{Int64, Int64} with 10 entries:
5 => 7
4 => 15
6 => 3
7 => 3
2 => 33
10 => 2
9 => 2
8 => 3
3 => 35
1 => 12The largest PSD constraint is now of size 10, which is much smaller than the original 124-by-124 matrix.