SDP relaxations: max-cut

Solves a semidefinite programming relaxation of the MAXCUT graph problem:

max   0.25 * L•X
s.t.  diag(X) == e
      X ≽ 0

Where L is the weighted graph Laplacian. Uses this relaxation to generate a solution to the original MAXCUT problem using the method from the paper:

Goemans, M. X., & Williamson, D. P. (1995). Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming. Journal of the ACM (JACM), 42(6), 1115-1145.

using JuMP
import LinearAlgebra
import Random
import SCS
import Test

"""
    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)
    # We now have `X ≈ `F.U * D² * F.U'` where:
    D = LinearAlgebra.Diagonal(sqrt.(F.S))
    # So `X ≈ U' * U` where `U` is:
    return (F.U * D)'
end

function solve_max_cut_sdp(num_vertex, weights)
    # Calculate the (weighted) Laplacian of the graph: L = D - W.
    L = LinearAlgebra.diagm(0 => weights * ones(num_vertex)) - weights
    # Solve the SDP relaxation
    model = Model(SCS.Optimizer)
    set_silent(model)
    # Start with X as the identity matrix to avoid numerical issues.
    @variable(
        model,
        X[i = 1:num_vertex, j = 1:num_vertex],
        PSD,
        start = (i == j ? 1.0 : 0.0),
    )
    @objective(model, Max, 1 / 4 * LinearAlgebra.dot(L, X))
    @constraint(model, LinearAlgebra.diag(X) .== 1)
    optimize!(model)
    @assert termination_status(model) == MOI.OPTIMAL
    opt_X = value(X)
    V = svd_cholesky(opt_X)
    # Generate random vector on unit sphere.
    Random.seed!(num_vertex)
    r = rand(size(V, 1))
    r /= LinearAlgebra.norm(r)
    # Iterate over vertices, and assign each vertex to a side of cut.
    cut = ones(num_vertex)
    for i in 1:num_vertex
        if LinearAlgebra.dot(r, V[:, i]) <= 0
            cut[i] = -1
        end
    end
    println("Solution:")
    print(" (S, S′) = ({")
    print(join(findall(cut .== -1), ", "))
    print("}, {")
    print(join(findall(cut .== 1), ", "))
    println("})")
    #  (S, S′)  = ({1}, {2, 3, 4})
    return cut, 0.25 * sum(L .* (cut * cut'))
end

function example_max_cut_sdp()
    println()
    println("Example 1:")
    #   [1] --- 5 --- [2]
    #
    # Solution:
    #  (S, S′)  = ({1}, {2})
    cut, cut_val = solve_max_cut_sdp(2, [0.0 5.0; 5.0 0.0])
    Test.@test cut[1] != cut[2]

    println()
    println("Example 2:")
    #   [1] --- 5 --- [2]
    #    |  \          |
    #    |    \        |
    #    7      6      1
    #    |        \    |
    #    |          \  |
    #   [3] --- 1 --- [4]
    #
    # Solution:
    #  (S, S′)  = ({1}, {2, 3, 4})
    W = [
        0.0 5.0 7.0 6.0
        5.0 0.0 0.0 1.0
        7.0 0.0 0.0 1.0
        6.0 1.0 1.0 0.0
    ]
    cut, cut_val = solve_max_cut_sdp(4, W)
    Test.@test cut[1] != cut[2]
    Test.@test cut[2] == cut[3] == cut[4]

    println()
    println("Example 3:")
    #   [1] --- 1 --- [2]
    #    |             |
    #    |             |
    #    5             9
    #    |             |
    #    |             |
    #   [3] --- 2 --- [4]
    #
    # Solution:
    #  (S, S′)  = ({1, 4}, {2, 3})
    W = [
        0.0 1.0 5.0 0.0
        1.0 0.0 0.0 9.0
        5.0 0.0 0.0 2.0
        0.0 9.0 2.0 0.0
    ]
    cut, cut_val = solve_max_cut_sdp(4, W)
    Test.@test cut[1] == cut[4]
    Test.@test cut[2] == cut[3]
    Test.@test cut[1] != cut[2]
    return
end

example_max_cut_sdp()

Example 1:
Solution:
 (S, S′) = ({1}, {2})

Example 2:
Solution:
 (S, S′) = ({2, 3, 4}, {1})

Example 3:
Solution:
 (S, S′) = ({2, 3}, {1, 4})