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}}:
3.75235 -10.1468 11.8536 … 8.2426 -2.93578 -1.7936
-10.1468 3.47152 -1.2742 -4.50988 -5.62554 4.1421
11.8536 -1.2742 -5.39872 -3.94134 -5.66674 3.24501
5.08674 -2.33596 -0.452489 0.619368 -2.55535 6.84733
8.83628 1.75814 4.58552 -2.76005 -4.13215 3.39184
4.19718 3.64755 8.93206 … -3.95832 0.445552 -5.50658
0.308163 -5.51088 0.424802 -11.2947 -1.35679 5.7731
8.2426 -4.50988 -3.94134 7.45319 5.59746 0.668617
-2.93578 -5.62554 -5.66674 5.59746 -3.41909 4.27924
-1.7936 4.1421 3.24501 0.668617 4.27924 -7.90178We 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.11 - 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.63ms
Iter: Objective: Primal Res: Dual Res: Rho:
1 -5.6886e+04 1.6340e+02 4.8007e+00 1.0000e-01
25 3.7049e-04 8.4997e-01 1.7500e-08 1.0000e-01
50 -4.7558e+00 1.6838e-05 2.3707e+01 4.9553e+03
75 -3.3266e+01 1.2198e-03 2.2344e+01 4.9553e+03
100 -7.3809e+01 1.2839e-03 1.3274e+00 1.0090e+02
125 -7.4070e+01 9.6939e-05 3.1665e-02 1.0090e+02
150 -7.4069e+01 2.0422e-15 5.6843e-13 1.0090e+02
------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 160 (incl. 10 safeguarding iter)
Optimal objective: -74.07
Runtime: 0.165s (165.19ms)opt_objective = JuMP.objective_value(model)74.06893065551017Now, 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])74.06893065551021Solve the dual
Alternatively, we can solve the dual problem:
model = JuMP.Model(optimizer_with_attributes(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)------------------------------------------------------------------
COSMO v0.8.11 - A Quadratic Objective Conic Solver
Michael Garstka
University of Oxford, 2017 - 2022
------------------------------------------------------------------
Problem: x ∈ R^{111},
constraints: A ∈ R^{220x111} (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)
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.6ms
Iter: Objective: Primal Res: Dual Res: Rho:
1 -3.1509e+02 2.7621e+01 7.8548e+00 1.0000e-01
25 7.4069e+01 1.3454e-07 2.0250e-08 1.0000e-01
------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 25
Optimal objective: 74.07
Runtime: 0.003s (2.84ms)opt_objective = JuMP.objective_value(model)74.068931629612This 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.