## Branch-and-bound algorithms for QAP

This is part 4 in a series of posts laying out a simple branch-and-bound solver for QAP in C#. Last time (several months ago!) I provided a simple bounding procedure for QAP.  I want to take a step back and give the high-level algorithm, then in subsequent posts I will fill in the implementation.

Branch-and-bound algorithms are a general framework for solving combinatorial optimization problems.  The idea is to repeatedly subdivide the search space into smaller subproblems of the same time.  For each subproblem we compute bounds, which may allow us to determine that the subproblem cannot lead to an improved solution to the original problem.

```  /// Solves QAP using branch-and-bound. ///
public class BranchAndBound {

private BranchAndBoundResults _results;

/// Create a new instance. ///
public BranchAndBound(Qap qap) {
_qap = qap;
}

/// Solve a QAP to optimality. ///
public void Solve() {

_results = new BranchAndBoundResults();
Stack stack = new Stack(); stack.Push(new BranchAndBoundNode(_qap)); while (stack.Count > 0) { BranchAndBoundNode node = stack.Pop(); if (node.LowerBound <= _results.Objective) { if (node_is_easy_to_solve) { EnumerateSolutions(node, _results); } else { node.LowerBound = GLB.Bound(node.Qap, U); if (node.LowerBound <= _results.Objective) { foreach (BranchAndBoundNode subNode in Branch(node, U).OrderBy(n => n.LowerBound)) { stack.Push(subNode); } } } } } } } /// A subproblem in a branch-and-bound algorithm. ///  public class BranchAndBoundNode { private readonly Qap _qap; /// The size of the subproblem. ///  public int Size { get { return _qap.Size; } } /// The lower bound of the subproblem. ///  public double LowerBound { get; set; } /// The QAP subproblem. ///  public Qap Qap { get { return _qap; } } } /// Branch-and-bound solution information. ///  public class BranchAndBoundResults { public int[] p { get; set; } public double Objective { get; set; } } ```

Let’s work through the pseudocode.  Our stack will contain the current set of subproblems that need to be solved.  We initialize the stack with our QAP instance, wrapping it inside a BranchAndBoundNode data  structure.  We’ll develop this data structure as we go, but from the code we can already see two uses:

• It stores the subproblem (a QAP),
• It stores the lower bound computed by the code from my last post.

The other thing we do is allocate a data structure which stores the results: the best assignment found so far, and the best objective value.  In the main loop of Solve() we repeatedly pop nodes from the stack.  We first compare the node’s lower bound to the best objective (also called the incumbent).  If the lower bound exceeds the incumbent value, there is no way this subproblem can lead to a better result, so we need not process it.  The incumbent value may have changed since the subproblem was created.  Then we check to see if the subproblem is “easy” – if so, we can simply enumerate over all possible solutions for the subproblem and update the results if needed.  Otherwise, we use our GLB code to compute a bound.  Again, we check the lower bound to see if we can eliminate (or fathom) the subproblem.  If not, we need to subdivide the node into subproblems, and push them onto the stack.

The top-level algorithm itself is not all that complicated!  It’s also quite general.  We have already described how to compute bounds, so all we need to do is specify:

• How to enumerate solutions for small problems.
• How to branch – i.e. how to subdivide BranchAndBoundNodes.

The motivation and code for enumerating and branching will be the subject of the next two posts.

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;
}
}
}```

## Quadratic Assignment Problems: branch and bound

In my last post I introduced quadratic assignment problems.  As I explained, the problem is to assign an equal number of facilities to locations, minimizing the transportation cost.  A compact way to write the problem is to store the “flows” in a matrix A, the “distances” in matrix B, and represent assignments as permutation matrices.  If you throw in a linear term C, then the problem can be written as: Branch-and-bound algorithms are the most common technique for finding optimal solutions for QAP.  The next few posts will focus on describing how branch-and-bound works and will result in a complete implementation of a QAP solver in C#.  Branch-and-bound is a “divide and conquer” solution technique for combinatorial problems. Branch-and-bound repeatedly divides a problem into simpler subproblems of the same type, which are then solved or further subdivided. The search space of a problem is partitioned into search spaces for its subproblems. The algorithm generates a search tree with the root corresponding to the original problem to be solved.  At each node in the search tree, a relaxed version of the subproblem is solved, typically by loosening the constraints of the subproblem. The solution to the relaxation is a lower bound to the original subproblem. If a lower bound at a given node exceeds the value of a previously known solution, then continued searching on this path cannot lead to an improved solution, and there’s no need to continue searching on that node’s subtree (this is called “fathoming”).

A little less formally, we take the search space and divide it into subregions.  Then for each piece, we solve a “relaxed” version of the problem to see if we can rule out the solution lying in that subregion.  If we can’t rule it out, we divide the subregion into pieces and repeat.  The “relaxation” part is bounding, the dividing is called branching.  The effectiveness of a branch-and-bound algorithm depends on 1) the strength of the bound, 2) how fast bounds can be computed, and 3) the branching decisions.  For example, if the bound is terrible we’ll never fathom subproblems and we’ll end up searching the entire tree.

So how do we apply branch-and-bound to QAP?  The solution to a QAP is an assignment of facilities to locations.  Therefore the subproblems are obtained by assigning some of the facilities to locations.  To branch, I can pick a facility and try assigning it to each of the available locations (or vice-versa).  As we will see, the logic for selecting which facility/location to branch on is extremely important.

Let’s start putting together some of the pieces of our branch-and-bound solver for QAP.  First, let’s define a QAP, which as I said before consists of the flow and distance matrices A and B, and a linear term C.

```  // A quadratic assignment problem.
public class Qap {
// A (flow matrix).
public double[][] A { get; set; }
// B (distance matrix).
public double[][] B { get; set; }
// C (linear costs).
public double[][] C { get; set; }
// Constant term.
public double Shift { get; set; }

// Size of (A, B, C) matrices.
public int Size {
get { return A == null ? 0 : A.Length; }
}

// Compute tr AXBX' + CX', where the permutation matrix X is determined by p.
public double Evaluate(int[] p) {
double obj = this.Shift;
for (int i = 0; i < A.Length; i++) {
obj += this.C[i][p[i]];
for (int j = 0; j < A[i].Length; j++) {
obj += this.A[i][j] * this.B[p[i]][p[j]];
}
}
return obj;
}
}```

To make the code shorter I am using C# 3.0 automatic properties. In addition to the matrices I have also added a constant term “Shift”, whose purpose will become clear in a moment.  I’ve also provided a method called Evaluate() that computes the objective given a permutation.  Let’s think about how branching works.  If we go back to the mathematical formulation of QAP at the start of this post, assigning a facility to a location means setting X_ij = 1.  If we plug X_ij into the formula above we see that we end up with a QAP of size n – 1 with:

A′ = A(ii)
B′ = B(jj)
C′ = C(ij) + 2 a_i b_j ,
Shift’ = A[i][i] B[j][j] + C[i][j] + Shift

Where  A(ij) denotes matrix A with row i and column j removed, a_i are the elements in row i of A, excluding A[i][i], and let b_j are the elements in column j of B, excluding b[j][j]. Here’s the code:

```    private static Qap Reduce(Qap qap, int i, int j) {
Qap r = new Qap();

r.A = Reduce(qap.A, i, i);
r.B = Reduce(qap.B, j, j);
r.C = Reduce(qap.C, i, j);

// Tack on: 2 a'_i b'_j (where a'_i is row i of A, with element i removed)
for (int ii = 0; ii < r.Size; ii++) {
for (int jj = 0; jj < r.Size; jj++) {
r.C[ii][jj] += 2 * qap.A[i][ii < i ? ii : ii + 1] * qap.B[jj < j ? jj : jj + 1][j];
}
}

r.Shift = qap.Shift + qap.A[i][i] * qap.B[j][j] + qap.C[i][j];

return r;
}

private static double[][] Reduce(double[][] M, int i, int j) {
double[][] R = LinearAlgebra.NewMatrix(M.Length - 1);
for (int ii = 0; ii < R.Length; ii++) {
for (int jj = 0; jj < R[ii].Length; jj++) {
R[ii][jj] = M[ii < i ? ii : ii + 1][jj < j ? jj : jj + 1];
}
}
return R;
}```

LinearAlgebra.NewMatrix() is a helper function that initializes the matrix. So now we have a data structure for QAP, and we now how to obtain subproblems.  Next time we’ll look at how to compute lower bounds.

Update 6/9/2010: Here is an additional method to read a problem in QAPLIB format.  It does not attempt to read the “C” matrix, and it is not totally robust. But it should work.

```    public static Qap Read(StreamReader reader) {
var tokens = content.Split(new char[] { ' ', '\n', '\r'}, StringSplitOptions.RemoveEmptyEntries);
int k = 0;
int size;
if (Int32.TryParse(tokens[k++], NumberStyles.Integer, NumberFormatInfo.InvariantInfo, out size)) {
Qap qap = new Qap();
qap.A = ReadMatrix(size, size, tokens, ref k);
qap.B = ReadMatrix(size, size, tokens, ref k);
return qap;
}
throw new InvalidDataException("Invalid QAP size.");
}

private static double[][] ReadMatrix(int m, int n, string[] tokens, ref int k) {
if (k + m * n > tokens.Length) {
throw new InvalidDataException("Not enough matrix data.");
}
double[][] matrix = LinearAlgebra.NewMatrix(m, n);
for (int i = 0; i < matrix.Length; i++) {
for (int j = 0; j < matrix[i].Length; j++) {
if (!Double.TryParse(tokens[k++], NumberStyles.Float, NumberFormatInfo.InvariantInfo, out matrix[i][j])) {
throw new InvalidDataException("Invalid matrix entry.");
}
}
}
return matrix;
}```

## Quadratic Assignment Problems, Part I

In a previous series I talked a bit about linear assignment problems (LAPs) and provided code for a simple LAP solver.  Now I want to spend a few posts on LAP’s big brother, the Quadratic Assignment Problem.  In an LAP we are trying to assign objects in one set to an equal number of objects in another set, where we are given the cost of assigning any pair of objects.  In QAP we still have two sets of objects – and I will stick with the original formulation of the problem (from the 50’s) and call the first step of objects “facilities” and the second set “locations”.  Imagine that we are going to build a bunch of computer labs in different rooms of a building, for example.  Suppose that each lab has a different purpose, and the rooms are wired up and we’re trying to figure out where to put each lab.  The quadratic assignment problem is the problem of finding an assignment
of an equal number of facilities to locations that minimizes the transportation cost between facilities.  In my example, the transportation cost is a product of how far the rooms are from each other and how often one travels from one lab to the other.

QAPs are interesting problems in part because they are so difficult to solve.  Finding the optimal assignment is an NP-hard problem.  In upcoming posts I will write down the problem more formally, and explore some heuristics for QAP, as well as some tools for finding exact solutions for QAP.

## Linear assignment problems, part III

Here’s the code. See the previous posting for background.

```using System;

namespace NathanBrixius.AssignmentLib
{
/// <summary>Linear assignment problem solver.<summary>
/// <remarks>
/// This code is a port of Fortran code by R.E. Burkard and U. Derigs.
/// Assignment and matching problems: Solution methods with fortran programs,
/// volume 184 of Lecture Notes in Economics and Mathematical Systems.
/// Springer, Berlin, 1980.  The original code is available on the QAPLIB site.
/// </remarks>
public class LinearAssignment
{
public static void Main(params string[] args)
{
double[][] L = new double[];
for (int i = 0; i < L.Length; i++)
{
L[i] = new double;
for (int j = 0; j < L[i].Length; j++)
{
L[i][j] = i + j;
}
}
int[] p = new int[L.Length];
int[] ip = new int[L.Length];
double[] ys = new double[L.Length];
double[] yt = new double[L.Length];
Solve(L, p, ip, ys, yt);
}

/// <summary>Evaluate a permutation against the given LAP.<summary>
/// <param name="L">LAP cost matrix.</param>
/// <param name="p">Permutation.</param>
/// <returns>The cost of the permutation.</returns>
public static double Evaluate(double[][] L, int[] p)
{
double cost = 0.0;
for (int i = 0; i < p.Length; i++)
{
cost += L[i][p[i]];
}
return cost;
}

/// <summary>Calculate the reduced costs matrix for the given LAP and solution.<summary>
/// <param name="L">LAP cost matrix.</param>
/// <param name="U">Reduced cost matrix (filled in by the method)</param>
/// <param name="ys">Dual variables.</param>
/// <param name="yt">Dual variables.</param>
public static void ReducedCosts(
double[][] L,
double[][] U,
double[] ys,
double[] yt)
{
for (int i = 0; i < L.Length; i++)
{
for (int j = 0; j < L[i].Length; j++)
{
U[i][j] = L[i][j] - (yt[j] + ys[i]);

}
}
}

/// <summary>Solve the specified linear assignment problem.<summary>
/// <param name="L">LAP cost matrix.</param>
/// <param name="p">The solution permutation.</param>
/// <param name="ip">The inverse of the solution permutation.</param>
/// <param name="ys">Dual variables.</param>
/// <param name="yt">Dual variables.</param>
/// <returns></returns>
public static double Solve(
double[][] L,
int[] p,
int[] ip,
double[] ys,
double[] yt
)
{
if (L.Length == 0)
{
return 0.0;
}
return SolveCore(L, ip, p, ys, yt, 1.0e-7);
}

/// <summary>Linear sum assignment problem with double precision costs.<summary>
/// <param name="L">Cost matrix.</param>
/// <param name="ip">Inverse of optimal assignment.</param>
/// <param name="p">Optimal assignment.</param>
/// <param name="ys">Optimal dual (row) variables.</param>
/// <param name="yt">Optimal dual (column) variables.</param>
/// <param name="tolerance">Machine accuracy.</param>
/// <returns>The optimal cost.</returns>
private static double SolveCore(
double[][] L,
int[] ip,
int[] p,
double[] ys,
double[] yt,
double tolerance
)
{
int[] vor = new int[L.Length];
bool[] label = new bool[L.Length];
double[] dMinus = new double[L.Length];
double[] dPlus = new double[L.Length];

// Start: Construction of an initial (partial) assignment.
// (partial --> p[i]/ip[j] may be BadIndex for some i,j)
double maxValue = 1.0e15;
for (int i = 0; i < p.Length; i++)
{
ys[i] = 0;
yt[i] = 0;
}

double ui = 0.0;
int jOpt = -1; // optimal j index
for (int i = 0; i < L.Length; i++)
{
for (int j = 0; j < L[i].Length; j++)
{
if ((j == 0) || (L[i][j] - ui < tolerance))
{
ui = L[i][j];
jOpt = j;
}
}

// guarantee: jOpt >= 0
ys[i] = ui;
{
ip[jOpt] = i;
p[i] = jOpt;
}
}

for (int j = 0; j < ip.Length; j++)
{
{
yt[j] = maxValue;
}
}

for (int i = 0; i < ys.Length; i++)
{
ui = ys[i];
for (int j = 0; j < yt.Length; j++)
{
double vj = yt[j];
if (vj > tolerance)
{
double cc = L[i][j] - ui;
if (cc + tolerance < vj)
{
yt[j] = cc;
vor[j] = i;
}
}
}
}

for (int j = 0; j < vor.Length; j++)
{
int i = vor[j];
{
p[i] = j;
ip[j] = i;
}
}

for (int i = 0; i < p.Length; i++)
{
{
ui = ys[i];
for (int j = 0; j < ip.Length; j++)
{
{
double lij = L[i][j];
if ((lij - ui - yt[j] + tolerance) <= 0)
{
p[i] = j;
ip[j] = i;
}
}
}
}
}

// Construction of the optimal assignment
double d = 0.0;
for (int u = 0; u < p.Length; u++)
{
{
// Shortest path computation
for (int i = 0; i < vor.Length; i++)
{
vor[i] = u;
label[i] = false;
dPlus[i] = maxValue;
dMinus[i] = L[u][i] - ys[u] - yt[i];
}
dPlus[u] = 0;

bool done = false;
while (!done)
{
d = maxValue;
for (int i = 0; i < dMinus.Length; i++)
{
if (!label[i] && (dMinus[i] + tolerance < d))
{
d = dMinus[i];
index = i;
}
}

// note: index >= 0
{
label[index] = true;
w = ip[index];
dPlus[w] = d;
for (int i = 0; i < label.Length; i++)
{
if (!label[i])
{
double vgl = d + L[w][i] - ys[w] - yt[i];
if (dMinus[i] > vgl + tolerance)
{
dMinus[i] = vgl;
vor[i] = w;
}
}
}
}
else
{
done = true;
}
}

// Augmentation
done = false;
while (!done)
{
w = vor[index];
ip[index] = w;
int tempIndex = p[w];
p[w] = index;
if (w == u)
{
done = true;
}
else
{
index = tempIndex;
}
}

// Transformation
for (int i = 0; i < dPlus.Length; i++)
{
if (dPlus[i] != maxValue)
{
ys[i] += d - dPlus[i];
}

if (dMinus[i] + tolerance < d)
{
yt[i] += dMinus[i] - d;
}
}
}
}

// Computation of the optimal value
return Evaluate(L, p);
}
}
}```

## Linear assignment problems, part II

In a previous post I talked a bit about the linear assignment problem:  finding the minimum cost matching in a bipartite graph.  In my next post I will give a simple self-contained C# program for solving LAP based on some Fortran code from my colleague Rainer Burkard.  First a few notes and disclaimers!

The input to the code is a matrix which represents the adjacency matrix of the graph.  In my application (solving quadratic assignment problems) I am concerned with input matrices that are relatively small and dense.   The code is a port of Fortran code by Burkard and Ulrich Derigs:  Assignment and matching problems: Solution methods with Fortran programs, volume 184 of Lecture Notes in Economics and Mathematical Systems.  Springer, Berlin, 1980.  The original code is available on the QAPLIB site.  Rainer has given me permission to post this code, although if you look at the terms on the original you will see that it is freely distributable.  This code is offered under the same conditions and I make no claims about its correctness, value, etc.  But it works okay for me.

My port is relatively “conservative” in that I have tried to retain the original structure and flow as much as possible because I don’t want to introduce any bugs.  I want to keep the implementation clean, so like many numerical optimization codes there isn’t much input validation.  I have chosen to use arrays rather than trying to encapsulate the input and output into matrix and vector classes.  I kept it simple so you can adapt it to your needs.

Lastly, I wanted to mention that this is not the most efficient LAP solver available.  At the time I wrote my thesis the fastest code over the domain I was working on was the code of Jonker and Volgenant, which is available on this site.  I notice that the code is used by others as well, see here and here.  This code is really quite slick and I had a great experience working with it all those years ago.  It’s not hard to port it over to C#, so go for it!

## Linear assignment problems, part I

I’m going to devote a few posts to my thesis topic, the quadratic assignment problem.  Quadratic assignment problems are a type of matching problem that are NP hard and have some interesting applications.  I’ll start out by describing linear assignment problems, which are simpler but still fun.

As the name suggests, in a linear assignment problem objects in one set are to be assigned to an equal number of objects in a second set. The cost of assigning objects in the first set to objects in the second set are given, and the objective is to find the assignment that minimizes the total assignment cost.  For example, the problem may be to assign workers to jobs so as to minimize the total amount of working time.  A more formal description of the problem can be found on this page (and I recommend the book, BTW).

LAP can also be posed as the problem of finding the minimum cost matching on a bipartite graph.  For this reason, in algorithms texts this problem is also sometimes referred to as “minimum cost weighted bipartite matching”.  LAPs are linear programs and so they can be solved in polynomial time.  As the wikipedia article on the problem states, there are algorithms that take advantage of the special structure of LAP to solve them more efficiently than general purpose LP solvers.

Next time I’ll go through one possible solution approach for LAP.