The source files for all examples can be found in /examples.
Minimizing the sum of the k-largest λ
We show how to find the sum of absolute value of the k largest eigenvalues of a symmetric matrix $A \in \mathbb{S}^n$. This problem can be solved as a semidefinite program. The primal and dual forms are stated in Alizadeh [1]:
\[\begin{array}{llll} \text{maximize} & \text{Tr}(AY) - \text{Tr}(AW) & \text{minimize} & kz + Tr(U) + Tr(V) \\ \text{subject to} & \text{Tr}(Y + W) = k & \text{subject to} & zI + V - A \succeq 0 \\ & 0 \preceq Y \preceq I & & zI + U + A \succeq 0 \\ & 0 \preceq W \preceq I & & U, V \succeq 0, \end{array}\]
where $Y, W$ are the variables of the primal and $U, V$ are the variables of the dual problem.
using LinearAlgebra, JuMP, COSMO, Random
rng = Random.MersenneTwister(212)
n = 10
A = 5 .* randn(rng, 10, 10)
A = Symmetric(A, :U)10×10 LinearAlgebra.Symmetric{Float64, Matrix{Float64}}:
11.7906 3.8406 -2.27707 … -5.49545 5.94112 -15.7009
3.8406 0.354475 -4.15618 -0.521508 9.31386 5.21261
-2.27707 -4.15618 0.0940684 3.97709 4.79056 2.02411
3.47812 8.49617 3.51557 -10.0605 1.27571 -5.10486
-4.46697 -1.10636 3.91133 0.921127 -2.03889 9.29371
0.814444 -1.12582 8.97819 … 3.19451 -2.19222 -2.27626
-6.08621 8.48159 -1.36896 1.34568 2.53833 -1.43798
-5.49545 -0.521508 3.97709 -3.00852 4.92574 7.4608
5.94112 9.31386 4.79056 4.92574 1.91581 -3.88751
-15.7009 5.21261 2.02411 7.4608 -3.88751 4.08394We are interested in minimizing the sum of absolute values of the k=3 largest eigenvalues. Let's formulate the problem in JuMP with COSMO as the backend solver:
k = 3
model = JuMP.Model(optimizer_with_attributes(COSMO.Optimizer, "verbose" => true));
@variable(model, Y[1:n, 1:n], PSD);
@variable(model, W[1:n, 1:n], PSD);
@objective(model, Max, tr(A * Y) - tr(A * W));
@constraint(model, tr(Y + W) == k);
@constraint(model, Symmetric(I - Y) in PSDCone());
@constraint(model, Symmetric(I - W) in PSDCone());
status = JuMP.optimize!(model)------------------------------------------------------------------
COSMO v0.8.9 - A Quadratic Objective Conic Solver
Michael Garstka
University of Oxford, 2017 - 2022
------------------------------------------------------------------
Problem: x ∈ R^{110},
constraints: A ∈ R^{221x110} (240 nnz),
matrix size to factor: 331x331,
Floating-point precision: Float64
Sets: DensePsdConeTriangle of dim: 55 (10x10)
DensePsdConeTriangle of dim: 55 (10x10)
DensePsdConeTriangle of dim: 55 (10x10)
DensePsdConeTriangle of dim: 55 (10x10)
ZeroSet of dim: 1
... and 0 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.14ms
Iter: Objective: Primal Res: Dual Res: Rho:
1 -6.3802e+04 1.7856e+02 6.3659e+00 1.0000e-01
25 -2.5676e+00 7.4811e-01 2.1078e-02 1.0000e-01
50 -9.5448e+01 1.3991e-01 1.9139e+00 1.3044e+01
75 -7.9281e+01 7.4365e-14 2.2302e-12 1.3044e+01
------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 89 (incl. 14 safeguarding iter)
Optimal objective: -79.28
Runtime: 0.01s (10.49ms)opt_objective = JuMP.objective_value(model)79.28110570563726Now, we can check the solution by computing the sum of the absolute value of the 3-largest eigenvalues:
k_λ_abs = sum(sort(abs.(eigen(A).values), rev = true)[1:k])79.28110570563712Solve the dual
Alternatively, we can solve the dual problem:
model = JuMP.Model(with_optimizer(COSMO.Optimizer, verbose=true));
@variable(model, V[1:n, 1:n], PSD);
@variable(model, U[1:n, 1:n], PSD);
@variable(model, z);
@objective(model, Min, k * z + tr(V) + tr(U));
@constraint(model, Symmetric(z .* diagm(0 => ones(n)) + V - A) in PSDCone());
@constraint(model, Symmetric(z .* diagm(0 => ones(n)) + U + A) in PSDCone());
status = JuMP.optimize!(model)opt_objective = JuMP.objective_value(model)79.28110570563726This gives the same result.
Problem with A as variable
Above problems are mostly helpful for illustrative purpose. It is obviously easier to find the sum of the k-largest eigenvalues by simply computing the eigenvalues of $A$. However, above results become useful if finding $A$ itself is part of the problem. For example, assume we want to find a valid matrix $A$ under the constraints: $C\, \text{vec}(A) = b$ with the minimum sum of absolute values of the k-largest eigenvalues. We can then solve the equivalent problem:
\[\begin{array}{ll} \text{minimize} & kz + Tr(U) + Tr(V) \\ \text{subject to} & C \text{vec}(A) = b \\ & zI + V - A \succeq 0 \\ & zI + U + A \succeq 0 \\ & U, V \succeq 0. \end{array}\]
References
[1] Alizadeh - Interior point methods in semidefinite programming with applications to combinatorial optimization (1995)
This page was generated using Literate.jl.