The source files for all examples can be found in /examples.
Relaxed Two-Way Partitioning Problem
We consider the (nonconvex) two-way partitioning problem from Boyd and Vandenberghe (2004), p.219 [1]:
\[\begin{array}{ll} \text{minimize} & x^\top W x \\ \text{subject to} & x_i^2 = 1, \quad i=1,\dots, n, \end{array}\]
with $x \in \mathbf{R}^n$ and $W \in \mathbf{S}^n$. The problem can be interpreted as finding a partition of the points $i = 1,\dots, n$ into two sets, where the cost of two points in the same set is $W_{ij}$ and $-W_{ij}$ otherwise. Brute forcing the solution $x^*$ takes $2^n$ attempts and becomes quickly intractable, e.g. for $n \geq 30$. However, a lower-bound for this problem can be computed by solving the convex dual of the problem:
\[\begin{array}{ll} \text{maximize} & -\sum_i \nu \\ \text{subject to} & W + \text{diag}(\nu) \in \mathbf{S}_+^n. \end{array}\]
Solving this problem with optimal objective $d^*$ yields a lower bound as least as good as a lower bound based on the minimum eigenvalue $d^* \geq \lambda_{\text{min}}(W)$. Let's set up the problem for $n = 20$ and with a specially structured (banded) $W$:
using LinearAlgebra, Random, COSMO, JuMP, IterTools
rng = Random.MersenneTwister(212313);
n = 20;
W = diagm(0 => randn(rng, n));
W += diagm(1 => randn(rng, n-1));
W = Symmetric(W)20×20 LinearAlgebra.Symmetric{Float64, Matrix{Float64}}:
-1.05284 0.307241 0.0 … 0.0 0.0 0.0
0.307241 0.0283543 0.802181 0.0 0.0 0.0
0.0 0.802181 -0.760236 0.0 0.0 0.0
0.0 0.0 0.584959 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.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 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.318299 0.0 0.0
0.0 0.0 0.0 0.503906 0.270842 0.0
0.0 0.0 0.0 0.270842 -1.21555 -1.86596
0.0 0.0 0.0 0.0 -1.86596 -1.09988As you can see, the matrix $W$ imposes a structure on the LMI-constraint. Accordingly, COSMO will be able to use chordal decomposition to decompose the LMI constraint into multiple smaller constraints. This can make a significant difference in the performance of the algorithm for large $n$.
model = JuMP.Model(COSMO.Optimizer);
@variable(model, ν[1:n]);
@objective(model, Max, -ones(n)' * ν )
@constraint(model, Symmetric(W + diagm(ν) ) in JuMP.PSDCone());
JuMP.optimize!(model)------------------------------------------------------------------
COSMO v0.8.11 - A Quadratic Objective Conic Solver
Michael Garstka
University of Oxford, 2017 - 2022
------------------------------------------------------------------
Problem: x ∈ R^{38},
constraints: A ∈ R^{57x38} (56 nnz),
matrix size to factor: 95x95,
Floating-point precision: Float64
Sets: PsdConeTriangle of dim: 3 (2x2)
PsdConeTriangle of dim: 3 (2x2)
PsdConeTriangle of dim: 3 (2x2)
PsdConeTriangle of dim: 3 (2x2)
PsdConeTriangle of dim: 3 (2x2)
... and 14 more
Decomp: Num of original PSD cones: 1
Num decomposable PSD cones: 1
Num PSD cones after decomposition: 19
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 -6.0024e+02 1.7158e+01 5.9999e-01 1.0000e-01
25 2.3990e+01 4.7665e-01 2.3740e-03 1.0000e-01
50 2.6420e+01 3.2828e-01 3.0839e-03 1.0000e-01
75 2.7686e+01 1.0593e-01 1.0940e-02 1.0000e-01
100 2.8268e+01 7.7224e-04 1.0000e+00 1.8313e+03
125 2.7999e+01 1.0593e-01 9.7858e-06 6.4077e-01
150 2.8149e+01 1.3654e-10 3.6637e-15 6.4077e-01
------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 174 (incl. 24 safeguarding iter)
Optimal objective: 28.15
Runtime: 0.195s (195.3ms)Looking at the solver output you can see how the PSD constraint was decomposed into 19 PSD constraints. Let's look at the lower bound:
JuMP.objective_value(model)-28.149260383754584As $n$ is small, we can verify our result by finding the optimal solution by trying out all possible combinations:
function brute_force_optimisation(W, n)
opt_obj = Inf
opt_x = Inf * ones(n)
for xt in Iterators.product([[1.; -1.] for k = 1:n]...)
x = [i for i in xt]
obj_val = x' * W * x
if obj_val < opt_obj
opt_obj = obj_val
opt_x = x
end
end
return opt_obj, opt_x
end
opt_obj, opt_x = brute_force_optimisation(W, n)(-28.14926038252595, [-1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0])References
[1] Boyd and Vandenberghe - Convex Optimization, Cambridge University Press (2004)
This page was generated using Literate.jl.