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)Maximize ScalarAffineFunction{Float64}:
0.0 - 1.0 ν[1] - 1.0 ν[2] - 1.0 ν[3] - 1.0 ν[4] - 1.0 ν[5] - 1.0 ν[6] - 1.0 ν[7] - 1.0 ν[8] - 1.0 ν[9] - 1.0 ν[10] - 1.0 ν[11] - 1.0 ν[12] - 1.0 ν[13] - 1.0 ν[14] - 1.0 ν[15] - 1.0 ν[16] - 1.0 ν[17] - 1.0 ν[18] - 1.0 ν[19] - 1.0 ν[20]
Subject to:
VectorAffineFunction{Float64}-in-Scaled{PositiveSemidefiniteConeTriangle}
┌ ┐
│-1.0528363306124202 + 1.0 ν[1] │
│0.43450460079681563 │
│0.02835428739065727 + 1.0 ν[2] │
│0.0 │
│1.1344555620998293 │
│-0.7602362063711745 + 1.0 ν[3] │
│0.0 │
│0.0 │
│0.8272574966701005 │
│-0.3209415773794654 + 1.0 ν[4] │
│0.0 │
│0.0 │
│0.0 │
│-0.6347750259939933 │
│-0.057547453373101486 + 1.0 ν[5]│
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│-0.21642154367708533 │
│1.9305238145486778 + 1.0 ν[6] │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│1.1299973966118715 │
│0.322146132811399 + 1.0 ν[7] │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│-0.6769163654348226 │
│-0.29178365714150484 + 1.0 ν[8] │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.31041697885051156 │
│-0.5610902671491463 + 1.0 ν[9] │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│-0.21626918917423013 │
│-0.9654005199504526 + 1.0 ν[10] │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.3282846997746994 │
│-0.8216864802613469 + 1.0 ν[11] │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│1.590643027786316 │
│0.8036336112584134 + 1.0 ν[12] │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│-0.10593008441102737 │
│-0.5521994921431328 + 1.0 ν[13] │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│0.0 │
│2.917733868611536 │
│0.03759227595991484 + 1.0 ν[14] │
│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.9482849268801556 │
│-1.3486136363403636 + 1.0 ν[15] │
│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.5069950085676354 │
│0.42053061837018363 + 1.0 ν[16] │
│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 │
│1.029204297794958 │
│0.1582417620725275 + 1.0 ν[17] │
│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.4501426417706077 │
│0.503905800867215 + 1.0 ν[18] │
│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.3830282797137492 │
│-1.2155534535643648 + 1.0 ν[19] │
│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 │
│-2.638863787749651 │
│-1.0998836346904342 + 1.0 ν[20] │
└ ┘ ∈ Scaled{PositiveSemidefiniteConeTriangle}(PositiveSemidefiniteConeTriangle(20))
------------------------------------------------------------------
COSMO v0.8.10 - 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.19ms
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.188s (188.47ms)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.