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.08394

We 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.8 - 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.21ms

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.7181e+01	1.3116e-01	2.2060e+00	1.2620e+01
75	-7.9281e+01	2.2097e-13	9.1979e-12	1.2620e+01

------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 89 (incl. 14 safeguarding iter)
Optimal objective: -79.28
Runtime: 0.014s (13.88ms)
opt_objective = JuMP.objective_value(model)
79.28110570563364

Now, 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.28110570563709

Solve 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.28110570563364

This 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.