Home > Optimization, QAP > Quadratic assignment problems: bounds

This is part 3 of what will be a long series about how to solve QAPs using branch-and-bound. My goal to end up with a simple branch-and-bound algorithm in C#.  In my last post I wrote down the formulation for QAP, wrote a QAP class in C#, and showed how to create subproblems for branching.  Let’s focus on lower bounds for QAP.

We get a lower bound for QAP by relaxing some of the constraints of the problem, and then solving the relaxed version of the problem.  Different bounds  are obtained by relaxing the constraints in different ways.  Unfortunately in most cases the relaxations require a lot of math and are pretty darn hard to blog about.  So I will pick one of the oldest bounds for QAP, called the “Gilmore Lawler” bound, and try to sketch the motivation behind it.  I’m not going to try and prove anything about the bound – the point is that the relaxation that we end up with is fairly easy to code, which is exactly what I’m after.  If you’re after a nitty-gritty derivation of the bound, there’s the original paper by Lawler from 1963 (!), or check out Cela’s book on QAP.

In the case of QAP, the constraint is that “X” in our equation is a permutation matrix.  Permutation matrices are matrices where the rows and columns sum up to 1 and each entry is either a 0 or a 1.  The idea behind the bound is to not be so strict about these constraints.  Since QAP is a minimization problem, dropping constraints means that you’re going to get a solution whose value is lower than it would otherwise have been:  a lower bound.   Loosening these constraints (and here’s where it gets hand-wavy) results in a linear assignment problem – the matrix for the LAP is formed from the input matrices (A, B, C).  Solving the linear assignment problem gives us GLB.  A few months back I wrote some code to solve linear assignment problems, so now let’s use it to compute GLB.  The code is below.  A few notes:

• Just like the last post, LinearAlgebra.NewMatrix(n) creates a new matrix of size n.  In general, when I make reference to a method in LinearAlgebra, I assume you can write it yourself!
• LinearAssignment.ReducedCosts computes what is known as the “dual matrix” for the solution.  I won’t get into the details – the important fact is that each element of the dual matrix U_ij is a kind of prediction of how much the objective function will increase when i assign facility “i” to location “j”.  In other words, if U_ij is really big, then assigning facility i to location j is probably a bad idea.  For this reason, the dual matrix (I’ll just call it the “U matrix”) is very useful in a branch-and-bound algorithm: it helps us make branching decisions.
• UPDATED 7/21/2010: I have broken out AssignmentMatrix into a separate method. Other QAP bounds, such as the one referenced in the paper in this post, use this matrix.
```using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace NathanBrixius.AssignmentLib {
/// <summary>Computes the Gilmore-Lawler bound.</summary>
public class GilmoreLawler {
/// <summary>Computes the Gilmore-Lawler bound.</summary>
/// <param name="qap">The QAP.</param>
/// <param name="p">Primal solution.</param>
/// <param name="U">Reduced costs matrix (dual solution).</param>
/// <returns>The Gilmore-Lawler bound of the QAP.</returns>
public static double Bound(Qap qap, int[] p, double[][] U) {
int n = qap.Size;

double[][] L = AssignmentMatrix(qap);

// Solve LAP(L) to get bound - U is the reduced costs matrix.
int[] ip = new int[n];
double[] ys = new double[n];
double[] yt = new double[n];
double bound = LinearAssignment.Solve(L, p, ip, ys, yt) + qap.Shift;
LinearAssignment.ReducedCosts(L, U, ys, yt);
return bound;
}

public static double[][] AssignmentMatrix(Qap qap) {
int n = qap.Size;
// copy all elements of B_i into br_i except B_ii
double[][] br = LinearAlgebra.NewMatrix(n, n - 1);
for (int i = 0; i < n; i++) {
for (int j = 0; j < i; j++) {
br[i][j] = qap.B[i][j];
}
for (int j = i + 1; j < n; j++) {
br[i][j - 1] = qap.B[i][j];
}
Array.Sort(br[i]);
Array.Reverse(br[i]); // descending
}

// L_ij = A_ii B_jj + _ + c_ij
// (A'_i is computed inline below)
double[][] L = LinearAlgebra.NewMatrix(n);
double[] ar_i = new double[n - 1];
for (int i = 0; i < n; i++) {
for (int j = 0; j < i; j++) {
ar_i[j] = qap.A[i][j];
}
for (int j = i + 1; j < n; j++) {
ar_i[j - 1] = qap.A[i][j];
}
Array.Sort(ar_i); // ascending
for (int j = 0; j < n; j++) {
double l_ij = qap.A[i][i] * qap.B[j][j] + qap.C[i][j];
for (int k = 0; k < n - 1; k++) {
l_ij += ar_i[k] * br[j][k];
}
L[i][j] = l_ij;
}
}
return L;
}
}
}```