Простые примеры полуопределенного программирования

В этом руководстве собраны примеры небольших конических программ из области полуопределенного программирования (SDP).

В этом руководстве используются следующие пакеты.

using JuMP
import LinearAlgebra
import Plots
import Random
import SCS
import Test

Максимальный разрез с использованием SDP

Задача о максимальном разрезе — классический пример из теории графов. В ней нужно разделить граф на два комплементарных множества так, чтобы вес ребер между ними был максимальным. Эта задача является трудной в нелинейном программировании, но можно получить приближенное решение с помощью ослабления полуопределенного программирования:

где — это взвешенный лапласиан графа, а — вектор единиц. Дополнительные сведения см. в работе Goemans1995.

"""
    svd_cholesky(X::AbstractMatrix, rtol)


Return the matrix `U` of the Cholesky decomposition of `X` as `U' * U`.
Note that we do not use the `LinearAlgebra.cholesky` function because it
requires the matrix to be positive definite while `X` may be only
positive *semi*definite.

We use the convention `U' * U` instead of `U * U'` to be consistent with
`LinearAlgebra.cholesky`.
"""

function svd_cholesky(X::AbstractMatrix)
    F = LinearAlgebra.svd(X)
    # Теперь мы имеем `X ≈ `F.U * D² * F.U'`, где:
    D = LinearAlgebra.Diagonal(sqrt.(F.S))
    # Следовательно, `X ≈ U' * U`, где `U`:
    return (F.U * D)'
end

function solve_max_cut_sdp(weights)
    N = size(weights, 1)
    # Вычисляем (взвешенный) лапласиан графа: L = D - W.
    L = LinearAlgebra.diagm(0 => weights * ones(N)) - weights
    model = Model(SCS.Optimizer)
    set_silent(model)
    @variable(model, X[1:N, 1:N], PSD)
    for i in 1:N
        set_start_value(X[i, i], 1.0)
    end
    @objective(model, Max, 0.25 * LinearAlgebra.dot(L, X))
    @constraint(model, LinearAlgebra.diag(X) .== 1)
    optimize!(model)
    @assert is_solved_and_feasible(model)
    V = svd_cholesky(value(X))
    Random.seed!(N)
    r = rand(N)
    r /= LinearAlgebra.norm(r)
    cut = [LinearAlgebra.dot(r, V[:, i]) > 0 for i in 1:N]
    S = findall(cut)
    T = findall(.!cut)
    println("Solution:")
    println(" (S, T) = ({", join(S, ", "), "}, {", join(T, ", "), "})")
    return S, T
end

solve_max_cut_sdp (generic function with 1 method)

Для следующего графа:

[1] --- 5 --- [2]

Решением будет (S, T) = ({1}, {2}):

S, T = solve_max_cut_sdp([0 5; 5 0])

([2], [1])

Для следующего графа:

[1] --- 5 --- [2]
 |  \          |
 |    \        |
 7      6      1
 |        \    |
 |          \  |
[3] --- 1 --- [4]

Решением будет (S, T) = ({1}, {2, 3, 4}):

S, T = solve_max_cut_sdp([0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0])

([1], [2, 3, 4])

Для следующего графа:

[1] --- 1 --- [2]
 |             |
 |             |
 5             9
 |             |
 |             |
[3] --- 2 --- [4]

Решением будет (S, T) = ({1, 4}, {2, 3}):

S, T = solve_max_cut_sdp([0 1 5 0; 1 0 0 9; 5 0 0 2; 0 9 2 0])

([1, 4], [2, 3])

Кластеризация методом k-средних с использованием SDP

Данное множество точек в необходимо распределить по кластерам.

Дополнительные сведения см. в работе Peng2007.

function example_k_means_clustering()
    a = [[2.0, 2.0], [2.5, 2.1], [7.0, 7.0], [2.2, 2.3], [6.8, 7.0], [7.2, 7.5]]
    m = length(a)
    num_clusters = 2
    W = zeros(m, m)
    for i in 1:m, j in i+1:m
        W[i, j] = W[j, i] = exp(-LinearAlgebra.norm(a[i] - a[j]) / 1.0)
    end
    model = Model(SCS.Optimizer)
    set_silent(model)
    @variable(model, Z[1:m, 1:m] >= 0, PSD)
    @objective(model, Min, LinearAlgebra.tr(W * (LinearAlgebra.I - Z)))
    @constraint(model, [i = 1:m], sum(Z[i, :]) .== 1)
    @constraint(model, LinearAlgebra.tr(Z) == num_clusters)
    optimize!(model)
    @assert is_solved_and_feasible(model)
    Z_val = value.(Z)
    current_cluster, visited = 0, Set{Int}()
    for i in 1:m
        if !(i in visited)
            current_cluster += 1
            println("Cluster $current_cluster")
            for j in i:m
                if isapprox(Z_val[i, i], Z_val[i, j]; atol = 1e-3)
                    println(a[j])
                    push!(visited, j)
                end
            end
        end
    end
    return
end

example_k_means_clustering()

Cluster 1
[2.0, 2.0]
[2.5, 2.1]
[2.2, 2.3]
Cluster 2
[7.0, 7.0]
[6.8, 7.0]
[7.2, 7.5]

Задача о корреляции

Даны три случайные величины , и , а также границы двух из трех коэффициентов корреляции:

Задача заключается в определении верхней и нижней границ третьего коэффициента корреляции .

Мы решаем SDP, используя следующее положительно полуопределенное свойство матрицы корреляции:

function example_correlation_problem()
    model = Model(SCS.Optimizer)
    set_silent(model)
    @variable(model, X[1:3, 1:3], PSD)
    S = ["A", "B", "C"]
    ρ = Containers.DenseAxisArray(X, S, S)
    @constraint(model, [i in S], ρ[i, i] == 1)
    @constraint(model, -0.2 <= ρ["A", "B"] <= -0.1)
    @constraint(model, 0.4 <= ρ["B", "C"] <= 0.5)
    @objective(model, Max, ρ["A", "C"])
    optimize!(model)
    @assert is_solved_and_feasible(model)
    println("An upper bound for ρ_AC is $(value(ρ["A", "C"]))")
    @objective(model, Min, ρ["A", "C"])
    optimize!(model)
    @assert is_solved_and_feasible(model)
    println("A lower bound for ρ_AC is $(value(ρ["A", "C"]))")
    return
end

example_correlation_problem()

An upper bound for ρ_AC is 0.8719220302987134
A lower bound for ρ_AC is -0.97799895940284

Задача о минимальном искажении

Этот пример взят из вычислительной геометрии, в частности из задачи вложения общего конечного метрического пространства в евклидово пространство.

Известно, что 4-точечное метрическое пространство, определяемое звездным графом:

  [1]
    \
     1
      \
      [0] —- 1 -- [2]
      /
     1
    /
  [3]

не может быть точно вложено в евклидово пространство любой размерности, где расстояния вычисляются по длине кратчайшего пути между вершинами. Вложение с сохранением расстояний потребовало бы, чтобы три конечных узла образовывали равносторонний треугольник с длиной стороны 2, а центральный узел (0) отображался в равноудаленную точку на расстоянии 1. Это невозможно, так как из неравенства треугольника в евклидовом пространстве следует, что все точки должны быть одновременно коллинеарными.

Здесь мы сформулируем и решим задачу SDP для вычисления наилучшего возможного вложения, то есть вложения , связывающего каждую вершину с вектором с минимальным искажением таким, что

для всех ребер графа, где ] — это расстояние в метрическом пространстве графа.

Любое вложение может описываться матрицей Грама , которая является положительно полуопределенной и такой, что

Элемент ] матрицы представляет внутреннее произведение на .

Поэтому мы налагаем ограничение

для всех ребер графа и минимизируем , в результате чего получается приведенная ниже формулировка SDP. Так как за исходную можно взять любую точку, мы фиксируем первую вершину в точке 0.

Дополнительные сведения см. в работах Matousek 2013, Linial 2002.

function example_minimum_distortion()
    model = Model(SCS.Optimizer)
    set_silent(model)
    D = [
        0.0 1.0 1.0 1.0
        1.0 0.0 2.0 2.0
        1.0 2.0 0.0 2.0
        1.0 2.0 2.0 0.0
    ]
    @variable(model, c² >= 1.0)
    @variable(model, Q[1:4, 1:4], PSD)
    for i in 1:4, j in (i+1):4
        @constraint(model, D[i, j]^2 <= Q[i, i] + Q[j, j] - 2 * Q[i, j])
        @constraint(model, Q[i, i] + Q[j, j] - 2 * Q[i, j] <= c² * D[i, j]^2)
    end
    fix(Q[1, 1], 0)
    @objective(model, Min, c²)
    optimize!(model)
    Test.@test is_solved_and_feasible(model)
    Test.@test objective_value(model) ≈ 4 / 3 atol = 1e-4
    # Восстановим вложение с минимальным искажением:
    X = [zeros(3) sqrt(value.(Q)[2:end, 2:end])]
    return Plots.plot(
        X[1, :],
        X[2, :],
        X[3, :];
        seriestype = :mesh3d,
        connections = ([0, 0, 0, 1], [1, 2, 3, 2], [2, 3, 1, 3]),
        legend = false,
        fillalpha = 0.1,
        lw = 3,
        ratio = :equal,
        xlim = (-1.1, 1.1),
        ylim = (-1.1, 1.1),
        zlim = (-1.5, 1.0),
        zticks = -1:1,
        camera = (60, 30),
    )
end

example_minimum_distortion()

Числа Ловаса

Число Ловаса графа, также известное как тета-функция Ловаса, — это число, лежащее между двумя важными и связанными числами, которые трудно определить вычислительными методами, а именно хроматическим и кликовым числами графа. Однако число Ловаса можно эффективно вычислить как оптимальное значение полуопределенной программы.

Рассмотрим пятиугольный граф:

     [5]
    /   \
   /     \
 [1]     [4]
  |       |
  |       |
 [2] --- [3]

с пятью вершинами и ребрами. Известно, что его число Ловаса точно равно и лежит между 2 (наибольший размер клики) и 3 (наименьшее число, необходимое для раскраски вершины).

Пусть являются целыми числами, такими, что . определяется как симметричная матрица с элементами и , равными 1, и всеми остальными элементами, равными 0. Пусть — это множество ребер графа; в этом примере содержит (1,2), (2,3), (3,4), (4,5), (5,1) и их транспонированные элементы. Число Ловаса можно вычислить с помощью следующей программы:

где — это матрица, заполненная единицами, а — матрица тождественности.

Дополнительные сведения см. в работах Barvinok 2002, Knuth 1994.

function example_theta_problem()
    model = Model(SCS.Optimizer)
    set_silent(model)
    E = [(1, 2), (2, 3), (3, 4), (4, 5), (5, 1)]
    @variable(model, X[1:5, 1:5], PSD)
    for i in 1:5
        for j in (i+1):5
            if !((i, j) in E || (j, i) in E)
                A = zeros(Int, 5, 5)
                A[i, j] = 1
                A[j, i] = 1
                @constraint(model, LinearAlgebra.dot(A, X) == 0)
            end
        end
    end
    @constraint(model, LinearAlgebra.tr(LinearAlgebra.I * X) == 1)
    J = ones(Int, 5, 5)
    @objective(model, Max, LinearAlgebra.dot(J, X))
    optimize!(model)
    Test.@test is_solved_and_feasible(model)
    Test.@test objective_value(model) ≈ sqrt(5) rtol = 1e-4
    println("The Lovász number is: $(objective_value(model))")
    return
end

example_theta_problem()

The Lovász number is: 2.2360678617948677

Устойчивые множества неопределенности

В этом примере вычисляется значение риска для множества неопределенности, основанного на данных. Доступны выражения в замкнутой форме для оптимального значения. Дополнительные сведения см. в работе Bertsimas 2018.

function example_robust_uncertainty_sets()
    R, d, 𝛿, ɛ = 1, 3, 0.05, 0.05
    N = ceil((2 + 2 * log(2 / 𝛿))^2) + 1
    c, μhat, M = randn(d), rand(d), rand(d, d)
    Σhat = 1 / (d - 1) * (M - ones(d) * μhat')' * (M - ones(d) * μhat')
    Γ1(𝛿, N) = R / sqrt(N) * (2 + sqrt(2 * log(1 / 𝛿)))
    Γ2(𝛿, N) = 2 * R^2 / sqrt(N) * (2 + sqrt(2 * log(2 / 𝛿)))
    model = Model(SCS.Optimizer)
    set_silent(model)
    @variable(model, Σ[1:d, 1:d], PSD)
    @variable(model, u[1:d])
    @variable(model, μ[1:d])
    @constraint(model, [Γ1(𝛿 / 2, N); μ - μhat] in SecondOrderCone())
    @constraint(model, [Γ2(𝛿 / 2, N); vec(Σ - Σhat)] in SecondOrderCone())
    @constraint(model, [((1-ɛ)/ɛ) (u - μ)'; (u-μ) Σ] >= 0, PSDCone())
    @objective(model, Max, c' * u)
    optimize!(model)
    @assert is_solved_and_feasible(model)
    exact =
        μhat' * c +
        Γ1(𝛿 / 2, N) * LinearAlgebra.norm(c) +
        sqrt((1 - ɛ) / ɛ) *
        sqrt(c' * (Σhat + Γ2(𝛿 / 2, N) * LinearAlgebra.I) * c)
    Test.@test objective_value(model) ≈ exact atol = 1e-2
    return
end

example_robust_uncertainty_sets()