Symmetry reducting a strong relaxation of the quadratic assigment problem

Here, we are going to show how to load a QAP from QABLib, formulate a strong semidefinite relaxation of it, symmetry reduce it, and finally solve it.

Quadratic assigment problems

QAPs are given by two quadratic matrices $A$ and $B$. The objective is to permute the rows and columns of $B$, such that the inner product between the matrices is minimized.

$\mathrm{QAP}(A,B) = \min_{\phi\in S_n} \sum_{i,j=1}^n a_{ij}b_{\phi(i)\phi(j)}$

QAPs are notoriously hard to solve exactly, but there exist strong polynomial time relaxations, such as the following semidefinite programming relaxation:

\[\begin{aligned} \min\enspace & \langle B\otimes A ,Y\rangle\\ \mathrm{s.t.}\enspace & \langle I_n\otimes E_{jj},Y\rangle=1 \text{ for }j\in [n],\\ & \langle E_{jj}\otimes I_n,Y\rangle=1 \text{ for }j\in [n],\\ & \langle I_n\otimes (J_n-I_n)+(J_n-I_n)\otimes I_n,Y\rangle =0, \\ & \langle J_{n^2},Y\rangle = n^2,\\ & Y\in D^{n^2}, \end{aligned}\]

But in practice this relaxation is often too big to be solved directly.

Loading the data of a QAP

using SparseArrays, LinearAlgebra

file = joinpath(@__DIR__, "esc16j.dat")
data = open(file) do file
    read(file, String)
end
data = split(data, [' ', '\n', '\r'], keepempty = false)

n = parse(Int64, data[1])
A = zeros(Int64, n, n)
B = zeros(Int64, n, n)

pos = 2
for x = 1:n
    for y = 1:n
        A[x, y] = parse(Int64, data[pos])
        global pos += 1
    end
end
for x = 1:n
    for y = 1:n
        B[x, y] = parse(Int64, data[pos])
        global pos += 1
    end
end

Building the SDP (in vectorized standard form)

n = size(A, 1)

# Objective
CPrg = sparse(kron(B, A))

In = sparse(I, n, n)
Jn = ones(n, n)

# Vectorised constraint matrices as rows of large matrix APrg
APrg = spzeros(2n + 1, n^4)
bPrg = zeros(2n + 1)
currentRow = 1

for j = 1:n
    Ejj = spzeros(n, n)
    Ejj[j, j] = 1.0
    APrg[currentRow, :] = vec(kron(In, Ejj))
    bPrg[currentRow] = 1
    global currentRow += 1
    # Last constraint is linearly dependent on others
    if (j < n)
        APrg[currentRow, :] = vec(kron(Ejj, In))
        bPrg[currentRow] = 1
        global currentRow += 1
    end
end

APrg[currentRow, :] = vec(kron(In, Jn - In) + kron(Jn - In, In))
bPrg[currentRow] = 0
currentRow += 1
APrg[currentRow, :] = vec(ones(n^2, n^2))
bPrg[currentRow] = n^2

CPrg = sparse(vec(0.5 * (CPrg + CPrg')));

Symmetry reducing the SDP

We first determine an optimal admissible partition subspace

using SDPSymmetryReduction
P = admissible_subspace(CPrg, APrg, bPrg; verbose=true)
dim(P)

And then we block-diagonalize it

blkD = blockDiagonalize(P, true);

[ Info: Determining eigen-decomposition over Float64...
[ Info:   0.034041 seconds (76 allocations: 2.616 MiB)
[ Info: Determining the algebra-isomorphism...
[ Info:   0.002226 seconds (334 allocations: 550.938 KiB)
[ Info: Calculating image of the basis of the algebra...
[ Info:   0.011301 seconds (5.94 k allocations: 1.516 MiB)

Determining the coefficients of the reduced SDP

PMat = spzeros(Bool, n^4, dim(P))
for c in eachindex(P.matrix)
    PMat[c, P.matrix[c]] = 1
end

newA = APrg * PMat
newB = bPrg
newC = CPrg' * PMat;

Solving the reduced SDP with JuMP and CSDP

using JuMP, CSDP, MutableArithmetics
m = Model(CSDP.Optimizer)

# Initialize variables corresponding parts of the partition P
# >= 0 because the original SDP-matrices are entry-wise nonnegative
x = @variable(m, x[1:dim(P)] >= 0)

@constraint(m, newA * x .== newB)
@objective(m, Min, newC * x)

psdBlocks = @rewrite(sum(x[i] * blkD.blks[i] for i = 1:dim(P)));

for blk in psdBlocks
    if size(blk, 1) > 1
        @constraint(m, blk in PSDCone())
    else
        @constraint(m, blk .>= 0)
    end
end

optimize!(m)

Iter: 15 Ap: 9.55e-01 Pobj:  1.5743403e+01 Ad: 9.56e-01 Dobj:  1.5743403e+01
Success: SDP solved
Primal objective value: 1.5743403e+01
Dual objective value: 1.5743403e+01
Relative primal infeasibility: 7.80e-11
Relative dual infeasibility: 1.98e-09
Real Relative Gap: -4.40e-09
XZ Relative Gap: 3.48e-09
DIMACS error measures: 7.80e-11 0.00e+00 2.55e-08 0.00e+00 -4.40e-09 3.48e-09
CSDP 6.2.0
Iter:  0 Ap: 0.00e+00 Pobj: -7.0720000e+04 Ad: 0.00e+00 Dobj:  0.0000000e+00
Iter:  1 Ap: 3.30e-01 Pobj: -7.6764746e+04 Ad: 1.00e+00 Dobj:  5.2143222e+02
Iter:  2 Ap: 9.31e-01 Pobj: -7.9106606e+03 Ad: 1.00e+00 Dobj:  5.4816466e+02
Iter:  3 Ap: 9.28e-01 Pobj: -7.8457045e+02 Ad: 1.00e+00 Dobj:  5.5461228e+02
Iter:  4 Ap: 8.82e-01 Pobj: -1.6043789e+02 Ad: 1.00e+00 Dobj:  4.4338310e+02
Iter:  5 Ap: 6.36e-01 Pobj: -9.4557913e+01 Ad: 1.00e+00 Dobj:  2.0788282e+02
Iter:  6 Ap: 5.84e-01 Pobj: -5.2224326e+01 Ad: 9.82e-01 Dobj:  1.2940404e+02
Iter:  7 Ap: 7.28e-01 Pobj: -2.9510313e+01 Ad: 8.77e-01 Dobj:  1.2176895e+02
Iter:  8 Ap: 4.98e-01 Pobj: -2.1989634e+01 Ad: 1.00e+00 Dobj:  4.6134467e+01
Iter:  9 Ap: 6.76e-01 Pobj: -1.6046716e+01 Ad: 1.00e+00 Dobj:  1.7936308e+01
Iter: 10 Ap: 8.57e-01 Pobj: -1.0869645e+01 Ad: 1.00e+00 Dobj:  1.5229224e+00
Iter: 11 Ap: 7.61e-01 Pobj: -9.1547052e+00 Ad: 1.00e+00 Dobj: -4.9762779e+00
Iter: 12 Ap: 8.31e-01 Pobj: -8.2435512e+00 Ad: 1.00e+00 Dobj: -6.8076747e+00
Iter: 13 Ap: 7.28e-01 Pobj: -8.0063078e+00 Ad: 1.00e+00 Dobj: -7.4390623e+00
Iter: 14 Ap: 8.12e-01 Pobj: -7.8627214e+00 Ad: 1.00e+00 Dobj: -7.6741118e+00
Iter: 15 Ap: 8.12e-01 Pobj: -7.8129467e+00 Ad: 1.00e+00 Dobj: -7.7702774e+00
Iter: 16 Ap: 8.51e-01 Pobj: -7.7984854e+00 Ad: 1.00e+00 Dobj: -7.7881841e+00
Iter: 17 Ap: 9.56e-01 Pobj: -7.7947172e+00 Ad: 9.94e-01 Dobj: -7.7935255e+00
Iter: 18 Ap: 9.90e-01 Pobj: -7.7942416e+00 Ad: 9.95e-01 Dobj: -7.7945783e+00
Iter: 19 Ap: 1.00e+00 Pobj: -7.7942211e+00 Ad: 9.98e-01 Dobj: -7.7945005e+00
Iter: 20 Ap: 1.00e+00 Pobj: -7.7942188e+00 Ad: 1.00e+00 Dobj: -7.7942311e+00

termination_status(m)

OPTIMAL::TerminationStatusCode = 1

objective_value(m)

7.794218656165488

This page was generated using Literate.jl.