## Quadratic Assignment Problems for Bartenders

(Thanks to Jeff Linderoth and Jean-Pierre Goux for coming up with the “bar30” QAP that is the inspiration for this story.)

So this Guy walks into a bar.

Guy was working on his business degree, but he needed a part time job. He’d applied all over town, and was his last shot. At the back of the dimly lit bar, a silhouette. “Can, I help you, son?” An Irish brogue. Clearly a Master bartender.

“Looking for a job.”

“Do you have experience tending bar?”

“Absolutely.”

“Don’t try and fool me, son.”

“Well, maybe once at a wedding.”

“I like your style, kid. C’mon back.”

Guy slid around the back of the bar. On top of the bar were two cardboard boxes containing 30 bottles of various sorts, and a huge pile of receipts. On the counter behind him there was a wooden 6 x 5 lattice: 30 slots for 30 bottles.  “You want a job? Your first job is to get these bottles into these here holes. But do it right! We’re going to be busy tonight, and I don’t want to be dancing all over the place. Get those bottles in the right place!” Without another word, the Master went in to a small office in the back, to do who knows what.

Guy, an earnest fellow, thought about his task for a moment. “What does he mean, ‘do it right’? I guess he means place the bottles that he uses at the same time close to each other. That way he won’t have to dance around to get his job done.”

Guy started to ask himself how he’d ever figure out which bottles the Master uses together. Just then he noticed the huge pile of receipts. “I can look at these receipts to see what kinds of drinks people are ordering. Then I can look in my bartender’s guide to see which spirits I mix to make each drink. If I keep a log of all this then I can figure out how often, say, vodka is mixed with Red Bull, or 7-Up and grenadine.”  Guy smiled to himself. “I can do this. If I put the bottles in the holes, then I can figure out how far on average the Master will have to move. For each pair of bottles, I multiply the frequency that the bottles are used together in a drink with the distance between the bottles. If I add that up for all the pairs, I will know how far the Master will have to move around on average.”

He thought for a moment.

“Crap. How am I ever going to find the best arrangement? How many different ways can I stick the bottles into the holes, anyway? Well, if there were two slots and two bottles, it would be easy: 1 2 or 2 1. If there were three, then there would be 3 x 2 = 6 ways: 1 2 3, 1 3 2, 2 1 3, 2 3 1, 3 1 2, 3 2 1. So if there are 30 then its 30 x 29 x … x 2.” Musing to himself he realized that was 2.6525285981219105863630848e+32 combinations, fewer than the number of atoms in the universe, but still an exceptionally large number. Not even his iPhone 4 could carry out such a computation by brute force.

At this point, a reasonable person would simply do their best. They’d do something simple, like find the most frequently used bottle and stick it in the middle, and then put its frequently used friends as close as possible. But Guy was not a reasonable man. Nothing meant more to him than getting – and keeping – this job, and he already knew from experience that nothing gets past the Master. But Guy was stuck. At that very moment, a tall rakishly handsome man holding a black book entered the room. You guessed it, it was me.

Guy explained his problem to me. I looked him straight in the eyes. “Dude, you’ve got a QAP.”

“A what?”

“A quadratic assignment problem. You want to put each bottle in each hole. That’s an assignment. You say you want to do it the best way you can. It sounds like ‘best’ means move as little as possible, and the amount of movement depends on the frequency of use times the distance between the holes. When you are multiplying unknowns together, that’s quadratic. And you’ve definitely got a problem.”

“Okay, so what do I do?”

“Divide and conquer, baby. Take a bottle and stick it in a hole.”

“How do I know which one? And where? What if I get it wrong?”

“Just stick a bottle in a hole.”

“Now what do you have?”

“Well, I’ve got the Smirnoff in the upper-left-hand corner.”

“Yes, and 29 empty holes and 29 more bottles to place. Think about it: that’s another QAP, just a smaller one. Let’s assume you can find the best solution for the smaller problem. Now imagine putting the Smirnoff in all of the other holes, and solving those size 29 QAPs. If you know how to solve size 29 QAPs, then you will have your answer.”

“I know where you’re going with this, and I appreciate the advice, but I don’t see how that helps. I don’t know how to solve size 29 QAPs, except by solving a bunch of size 28 QAPs. And I can’t keep doing this Russian doll stuff because there are too many combinations.”

“True. But suppose you knew how rule out arrangements before you placed all the bottles. For example, say you try your best and come up with an arrangement where the Master consumes 800 calories a night. You know you can do at least that good. So suppose you have stuck 10 bottles in holes, but you’ve done a really bad job. In fact, you tally up the known cost with the 10 bottles and you’re already over 800 calories, even with 20 bottles left to go. Then you know that there is no way that any solution involving that exact placement of the first 10 bottles is going to be the best. Ruling out partial placements is called ‘fathoming‘.”

“Cool. But it still seems like we’re going to have to try out an awful lot of combinations.  I mean, how often can we really rule things out?”

“Well, you probably won’t rule out too many that way. What you need to do is figure out better rules for fathoming. You need to deduce, given a partial placement of bottles, the minimum number of calories that any complete placement will take. And you want that deduction – that lower bound – to be as good as possible. If it is a good lower bound, then you will be able to eliminate partial solutions much sooner – maybe even after placing 4 or 5 bottles. If you can do that, then you may end up with a reasonable number of possibilities to check.”

“I have no idea how to come up with a good lower bound.”

“Here, read Chapter 2.”

I made myself an Arnold Palmer as Guy proceeded to read Chapter 2 from the black book.

“That’s one way to do it. There are others. Some of them go way back, to like, the 60s. The key is that you’ve got to do it every time you place a bottle, so whatever you do, you’d better do it fast.”

“Fast. I can handle fast.”

“One more tip.”

“What’s that?”

“When you are picking bottles to place – branching – be careful.”

“Why is that?”

“Because the bound depends on the partial placement of the bottles. If you take the two bottles that are mixed together the most frequently – say vodka and Red Bull – I notice we are near a university – and put them in opposite corners, you will probably get a pretty large lower bound. That’s good if it helps you rule out partial solutions quickly.  What you want to do when you pick a bottle is make a selection that is likely to increase the lower bounds in the next round. Then you can rule out lots of possibilities. Sometimes the bounding procedure itself will give you hints. Do you know about dual values? Shadow prices?”

“Hell, no.”

“Never mind. Just read Section 3.5. And do what it says.” I hand the book to Guy, and start to edge towards the door.

“One more thing. You realize you don’t actually have to try putting the Smirnoff in all 30 holes, right?”

“Seriously?”

“If you’re putting the Smirnoff in the upper left hand corner, you don’t need to try putting it in the upper right hand corner. Whatever solution you come up with on the left side case you can flip and use for the right side case. Symmetry.”

“Symmetry. OK…”

“Section 3.6. Although I hear folks are into orbital branching these days…anyway…I’ve got to go. You can keep the book, I’ve got plenty. Hey – one more thing. And seriously, this is the last thing.”

“What’s that?”

“This ain’t shuffleboard on a Carnival Cruise. You’re going to need some help, even with all of these tricks. But most of what I told you about can be done in parallel.”

“How’s that?”

“Remember when I told you to try taking a bottle and sticking it in each of the possible holes? If you had a friend, and he tried half, and you tried half…. It’s more complicated than that, but that’s the idea.”

“Got it.” Guy looked through a doorway, where he saw another small bar with a virtually identical setup. “Do you think….?”

“Absolutely. What are friends for.”

“Thanks. Hey, now I’ve got a question for you. I’m glad it worked out for me, but QAP sounds a little bit, you know, contrived. I mean, how often do they come up in real life, anyway?”

“Well, they don’t bump into you on the street, but contrived is a matter of opinion…I mean, this whole scene reminds me of the movie ‘Cocktail’. I mean, there’s an Irish bartender and everything. Isn’t that contrived too? And speaking of Tom Cruise movies, have you heard the theory about Top Gun?

“No, I have not.”

“Top Gun is actually a thinly veiled Quadratic Assignment Problem tutorial. Inside the jets: instrument panels. How to lay out the instruments on the panel? QAP. How to do the wiring on the circuit boards for the fancy targeting systems? QAP. How to balance the turbines? QAP. How to allocate crews to flight decks? QAP. They could hardly make it more obvious. Days of Thunder? Same deal.”

“Cocktail, Top Gun, Days of Thunder. I suppose you’re going to tell me that Interview With the Vampire is all about QAP too.”

“You don’t even want to know.”

## Mixed integer linearizations for QAP using Solver Foundation

I read an interesting paper by Zhang, Beltran-Royo and Ma on the excellent Optimization Online site this morning. As an exercise I wrote a Solver Foundation program to reproduce some of the results.  In previous posts I showed how to write a C# program to solve Quadratic Assignment Problems (QAP). I wrote everything from scratch, including the Gilmore-Lawler bound computation. In this post (and the ones that follow) I am going to write everything using Solver Foundation Services in an effort to show how much easier it is!

If you are coming at this from the QAP angle, I recommend reading my previous QAP posts above. If you are coming at this from the Solver Foundation angle, look at page 4 of the paper (equations 2.12 – 2.15 in particular) and look at how easy it is to write down this rather complicated model using Solver Foundation Services. If you want a short Solver Foundation, check this out, or check out our new and improved documentation.

The paper talks about ways to use mixed integer linear programs (such as those built into Solver Foundation) to solve linearizations of QAP. One popular linearization is the one proposed by Kaufman and Broeckx. This formulation introduces n^2 continuous variables in addition to the standard integer variables x_ij (which represent the optimal assignment of facilities to locations). If you look at equations 2.12 – 2.15 on page 4, you will see that the goal function involves auxiliary real variables z.  2.13 is the most important constraint.  2.14 just enforces the non-negativity of the z variables, and we also need to make sure that the standard constraints on x hold: the row and column sums are equal to 1.

The code to compute the Kaufman-Broeckx linearization is below. I hardcoded a sample problem in A and B. You can replace it to read a standard problem from QAPLIB. Notice that the goals and constraints are basically transcriptions of the equations! The equations are complicated, but the way they map to the SFS APIs is straightforward. The other issue is how to use data binding to establish the values of the “q” and “a” parameters. This is done with the help of an auxiliary Cell class, and helper methods that return the “q” and “a” items. The SetBinding calls are pretty straightforward once you’ve got those.

If I were better at LINQ I could have simplified some of the code in GetCells and GetSums. But I’m not, so I didn’t.  Viewing the output shows that the goal function is 50 – which appears to be correct (that’s the optimal value for the nug05 problem). Looking at the x values you can determine the optimal permutation according to the linearization.

`using System;`
`using System.Collections.Generic;`
`using System.Linq;`
`using System.Text;`
`using Microsoft.SolverFoundation.Services;`
` `
`namespace QapLinearization {`
`  class Program {`
` `
`    // Here is the data for nug05.dat from QAPLIB.`
`    private double[][] A = new double[][] {`
`        new double[] { 0, 5, 2, 4, 1 },`
`        new double[] { 5, 0, 3, 0, 2 },`
`        new double[] { 2, 3, 0, 0, 0 },`
`        new double[] { 4, 0, 0, 0, 5 },`
`        new double[] { 1, 2, 0, 5, 0 }`
`      };`
`    private double[][] B = new double[][] {`
`        new double[] { 0, 1, 1, 2, 3 },`
`        new double[] { 1, 0, 2, 1, 2 },`
`        new double[] { 1, 2, 0, 1, 2 },`
`        new double[] { 2, 1, 1, 0, 1 },`
`        new double[] { 3, 2, 2, 1, 0 }`
`      };`
` `
`    static void Main(string[] args) {`
`      Program p = new Program();`
`      p.KaufmanBroeckx();`
`    }`
` `
`    public void KaufmanBroeckx() {`
`      SolverContext context = SolverContext.GetContext();`
`      Model model = context.CreateModel();`
` `
`      Set N = new Set(Domain.Integer, "N");`
` `
`      Parameter q = new Parameter(Domain.Real, "q", N, N, N, N);`
`      q.SetBinding(GetCells(), "Value", "I", "J", "K", "L");`
`      Parameter a = new Parameter(Domain.Real, "a", N, N);`
`      a.SetBinding(GetSums(), "Value", "I", "J");`
`      model.AddParameters(q, a);`
` `
`      Decision x = new Decision(Domain.IntegerRange(0, 1), "x", N, N);`
`      Decision z = new Decision(Domain.RealNonnegative, "z", N, N);`
`      model.AddDecisions(x, z);`
` `
`      model.AddGoal("g", GoalKind.Minimize,`
`        Model.Sum(Model.ForEach(N, i =>`
`          Model.Sum(Model.ForEach(N, j => z[i, j]))`
`          ))`
`          );`
` `
`      model.AddConstraint("c",`
`        Model.ForEach(N, i =>`
`          Model.ForEach(N, j =>`
`            z[i, j] >= Model.Sum(Model.ForEach(N, k =>`
`              Model.Sum(Model.ForEach(N, l => q[i, j, k, l] * x[k, l] - a[i, j] * (1 - x[i, j])))`
`              ))`
`        ))`
`      );`
` `
`      model.AddConstraint("rowsum",`
`        Model.ForEach(N, i =>`
`          Model.Sum(Model.ForEach(N, j => x[i, j])) == 1`
`      ));`
` `
`      model.AddConstraint("colsum",`
`        Model.ForEach(N, j =>`
`          Model.Sum(Model.ForEach(N, i => x[i, j])) == 1`
`      ));`
` `
`      var solution = context.Solve();`
`      Console.WriteLine(solution.GetReport());`
`    }`
` `
`    // This is a helper class for data binding.`
`    private class Cell {`
`      public int I { get; set; }`
`      public int J { get; set; }`
`      public int K { get; set; }`
`      public int L { get; set; }`
`      public double Value { get; set; }`
`    }`
` `
`    // This returns the q_ijkl as described on page 2 of the paper.`
`    private IEnumerable<Cell> GetCells() {`
`      for (int i = 0; i < A.Length; i++) {`
`        for (int j = 0; j < A[i].Length; j++) {`
`          for (int k = 0; k < B.Length; k++) {`
`            for (int l = 0; l < B[k].Length; l++) {`
`              yield return new Cell { I = i, J = j, K = k, L = l, Value = A[i][k] * B[j][l] };`
`            }`
`          }`
`        }`
`      }`
`    }`
` `
`    // This returns the a_ij as described on page 4 of the paper.`
`    private IEnumerable<Cell> GetSums() {`
`      for (int i = 0; i < A.Length; i++) {`
`        for (int j = 0; j < A[i].Length; j++) {`
`          double sum = 0;`
`          for (int k = 0; k < B.Length; k++) {`
`            for (int l = 0; l < B[k].Length; l++) {`
`              if (B[k][j] != 0) {`
`                sum += A[i][k] * B[j][l];`
`              }`
`            }`
`          }`
`          yield return new Cell { I = i, J = j, Value = sum };`
`        }`
`      }`
`    }`
`  }`
`}`

## Quadratic Assignment Problems: solution enumeration

In the last few posts we have given code for computing the Gilmore Lawler bound for quadratic assignment problems, and described two different branching techniques.   In this post we show how to enumerate solutions of a small QAP.  The possible solutions to a QAP range over permutations over unassigned facilities and locations.   If lower bounds do not increase after repeated branching, we may eventually assign nearly all the facilities to locations.  If there are only 2 or 3 remaining, we may as well just enumerate all the possible solutions and pick the best.   In practice, we hope that we will rarely have to do this, otherwise the branch-and-bound tree will be too large to process efficiently.  Here is a procedure that enumerates over all solutions of a size 3 QAP.  Iterating over permutations is a standard textbook problem, so maybe you can write cleaner code than I did here!

/// <summary>
/// Assume that all facilities have been assigned to locations,
/// except three. EnumerateSolutions enumerates all six remaining assignments
/// returns the best.
/// </summary>
private double EnumerateSolutions(Qap q, BranchAndBoundNode node, int[] pBest) {
double objBest;
int nPerm = perm3.GetLength(0);
int n = perm3.GetLength(1);
double[] obj = new double[nPerm];
int[] r = node.p.UnusedIndices().ToArray();
int[] c = node.pInv.UnusedIndices().ToArray();
int[] p = (int[])node.p.Clone();
Debug.Assert(r.Length == n, “Only call me for size “ + n + ” QAPs.”);

// Try all the possible permutations.
for (int i = 0; i < nPerm; i++) {
for (int j = 0; j < n; j++) {
p[r[j]] = c[perm3[i, j]];
}
obj[i] = q.Evaluate(p);
}

// Get the best permutation and objective.
int best = obj.ArgMin(out objBest);
p.CopyTo(pBest, 0);
for (int j = 0; j < n; j++) {
pBest[r[j]] = perm3[best, j];
}
return objBest;
}

If you combine the code from the following posts you will see that we almost have a working B&B solver:

In my next post I will fill in a couple of the gaps, and then I will re-post all of the code.  Then we’ll think about what we can actually do with a QAP solver!

## Quadratic Assignment Problems: strong and weak branching

Let’s pick up where we left off last time and write a score-based Branch delegate.  It will be used as the basis for both of our “real” branching functions.

private BranchingDecision BranchCore(BranchAndBoundNode node, double bound, double[][] S) {
RowColumnSums(S, node.Size);
double rowVal, colVal;
int rowBest = _rowSum.ArgMax(out rowVal);
int colBest = _colSum.ArgMax(out colVal);
if (rowVal > colVal) {
return new BranchingDecision(true, rowBest);
}
else {
return new BranchingDecision(false, colBest);
}
}

private void RowColumnSums(double[][] U, int size) {
if (_rowSum == null) {
_rowSum =
new double[_qap.Size];
_colSum =
new double[_qap.Size];
}
_rowSum.ConstantFill(
Double.MinValue);
_colSum.ConstantFill(
Double.MinValue);

for (int i = 0; i < size; i++) {
for (int j = 0; j < size; j++) {
_rowSum[i] += U[i][j];
_colSum[j] += U[i][j];
}
}
}

If we can create a reasonable score matrix and pass it into BranchCore, then we’ve got something.  Two common techniques for QAP are:

1. Weak branching. Let S_ij be the reduced cost U_ij that we got from computing the lower bound.  (U_ij is a lower bound on the amount that the bound on the subproblem will increase.  So large U_ij means the subproblem will be a lot easier.)
2. Strong branching. Let S_ij be the bound for subproblem ij.  The bigger the bound, the easier the subproblem.

Weak branching is trivial to implement for the Gilmore-Lawler (or GLB) bound.  As I mentioned a long time ago, the reduced cost matrix essentially comes for free.

private BranchingDecision WeakBranch(BranchAndBoundNode node, double bound, double[][] U) {
return BranchCore(node, bound, U);
}

For strong branching we need to compute bounds for each subproblem.  We overwrite the U matrix and pass it into BranchCore.  Since the U matrix is used to fathom subproblems later on in the code, we subtract the parent bound from the bound of each subproblem.

private BranchingDecision StrongBranch(BranchAndBoundNode node, double bound, double[][] U) {
Qap r = new Qap(node.Size – 1);
double[][] Ur = MatrixUtilities.NewMatrix(node.Size – 1, node.Size – 1);
for (int i = 0; i < node.Size; i++) {
for (int j = 0; j < node.Size; j++) {
Qap.Reduce(node.Qap, r, i, j);
U[i][j] =
GLB.Bound(r, Ur) – bound;
}
}
return BranchCore(node, bound, U);
}

Strong branching in general produces better guesses that weak branching.  However, it comes at a computational cost.  Therefore I want to introduce one last complication…er…feature.  Since strong branching is smarter but more costly, we should only use it where it is most needed.  I want to have a set of rules that governs when to use strong or weak branching.  Branching decisions are most important at the top of the tree, or on nodes where the bound sucks.  We can measure the suckiness of the bound for a node by computing the relative gap:  the ratio of the gap between the bound and the incumbent solution, and the corresponding gap at the root.   So our branching rule data structure will have Depth and Gap properties, as well as properties that store the Branch delegate information.

/// <summary>A description of a branching rule: when and how.
/// </summary>
public class BranchingRule {
/// <summary>The ID for this rule.
/// </summary>
public int Index { get; set; }
/// <summary>Maximum depth for which the rule applies.
/// </summary>
public int Depth { get; set; }
/// <summary>Maximum relative gap for which the rule applies.
/// </summary>
public double Gap { get; set; }
/// <summary>Strong or Weak.
/// </summary>
public BranchType Type { get; set; }
/// <summary>StrongBranch() or WeakBranch().
/// </summary>
internal Branch Rule { get; set; }

public BranchingRule(int index, double gap, int depth, BranchType type) {
Index = index;
Depth = depth;
Gap = gap;
Type = type;
}

public override string ToString() {
return “Depth = “ + Depth + “, Gap = “ + Gap + “, Type = “ + Type.ToString();
}
}

We will have an array of BranchingRule called _branchingRule.  Given a node, we will select the first rule in the array whose depth is larger than the node depth, and whose gap is larger than the node gap.

private Branch SelectBranchingRule(BranchAndBoundNode node) {
double gap = _results.RelativeGap(node.LowerBound);
BranchingRule criterion = _branchingRules.First(entry => (node.Level <= entry.Depth) && (gap >= entry.Gap));
node.RuleIndex = criterion.Index;
return criterion.Rule;
}

I am cheating here because I have not defined the local variable _results.  We’ll get to that in a few posts.  For now, just assume that RelativeGap is defined as I described above: (Objective – bound) / (Objective – RootBound).

At long last, that covers branching.  There are many variations, but they are essentially variations of strong or weak branching.  At this point we only have two more main topics to cover: what to do when you get to the bottom of the tree, and how do you capture and update results.

## Quadratic assignment problems: branching

(This is part of a long-running series on quadratic assignment problems in C#.  Click here to catch up.)

Branching means to divide the search space represented by a node into subnodes.  QAP is about making assignments of facilities to locations, so two possible approaches come to mind:

1. Branch into two subproblems: one where an assignment i -> j is allowed, and another where it is disallowed.
2. Pick a facility and try assigning it to all possible locations, creating N subproblems.  (Or pick a location and assign all facilities to it…)

The second approach, called polytomic branching, turns out to be more effective.  The problem with the first (“single assignment branching”) is that the subproblems aren’t that much easier than the parent problem.  You’ll have to branch a whole bunch of times before you are able to determine anything.

We still need to decide whether to fix a facility (“row branching”) or a location (“column branching”).  And after that we need to figure out which particular facility or location is best.  At this point you may be wondering, “does it even matter?”  The answer is that, yes, branching decisions are earth shatteringly important.  Where’s what I mean: if you make bad branching choices it could take a thousand times longer to solve the problem, or more!  This is part of the reason why QAP (and branch-and-bound) are so interesting.  Bad branching decisions lead to subproblems that are just about as hard as the original problem.  It’s kind of like this. Here is a data structure that represents the results of a polytomic branching decision.

/// <summary>A branching decision, the result of applying a branching rule.
/// </summary>
internal class BranchingDecision {
/// <summary>If true, enumerate over all unassigned rows.
/// </summary>
public bool IsRowBranch { get; set; }
/// <summary>The row or column index.
/// </summary>
public int Index { get; set; }
/// <summary>Create a new instance.
/// </summary>
public BranchingDecision(bool isRowBranch, int index) {
IsRowBranch = isRowBranch;
Index = index;
}
}

Ultimately we want to write methods that produce BranchingDecisions.  The following delegate represents such a method (we will write two methods in the next blog post):

/// <summary>Makes a branching decision.
/// </summary>
internal delegate BranchingDecision Branch(BranchAndBoundNode node, double bound, double[][] U);

If you accept the claim that branching choices make a big difference, then it is clear that we want to branch so we get the easiest subproblems possible.  That is, nodes that will be processed by the branch-and-bound algorithm as quickly as possible.  We don’t know how to measure that directly (why?), so we have to guess.  A general framework for polytomic branching is to compute a score for each of the N^2 subproblems, and store it in a matrix.  S_ij is the score for the subproblem created by fixing facility i to location j.  Then compute the row and column sums of S, and branch on the row or the column with maximum sum.  Making a branching choice reduces to computing scores on subproblems.

The next step will be to write methods that compute scores and return BranchingDecisions.  We’ll look at two common ways to do that next time.  We’re about two or three posts away from having a pretty decent branch-and-bound solver for QAP.

## Simple matrix, vector, permutation utilities

I am going to be rolling out the rest of my branch-and-bound algorithm in the next few posts.  To make that easier, in this post I introduce some common matrix, vector, and permutation methods.  It turns out that for technical computing applications, C#’s extension methods (introduced in 3.0) are awesome.  With vectors it’s great because you retain the control of having direct array access, but you get the nice object-oriented notation.

I’ve left out comments and error-checking, and it’s not super-optimized, but you get the idea.  None of these methods end up being the bottleneck.  Most of the methods are obvious, but I do want to comment on the permutation methods.  The goal of QAP is to find an optimal assignment of facilities to locations, represented as a permutation.  In the course of our branch-and-bound algorithm, we’ll be making partial assignments – that is, only some of the entries in the permutation will be filled in.  By convention, p[i] == -1 will mean that facility i is unassigned.  An index where p[i] < 0 is called “unused”.  When we branch, we’ll want to pick an unused facility (or location), and try assigning all unused locations (or facilities) to it.   The use of C#’s iterator and yield concepts really help here.

That said, here’s the code.

```public static class MatrixUtilities {
#region Vector
public static void ConstantFill(this T[] data, T val)   {     for (int i = 0; i < data.Length; i++) {       data[i] = val;     }   }   public static int ArgMin(this double[] data, out double best) {     if (data.Length == 0) {       best = Double.MinValue;       return -1;     }     int iBest = 0;     best = data;     for (int i = 1; i < data.Length; i++) {       if (best > data[i]) {         best = data[i];         iBest = i;       }     }     return iBest;   } ```
`  public static int ArgMax(this double[] data, out double best) {     if (data.Length == 0) { best = Double.MinValue; return -1; }     int iBest = 0;     best = data;     for (int i = 1; i < data.Length; i++) {       if (best < data[i]) {         best = data[i];         iBest = i;       }     }     return iBest;   } #endregion #region Permutation public static void Swap(this int[] p, int i, int j) {  int temp = p[i];   p[i] = p[j];   p[j] = temp; } public static int FindUnused(this int[] data, int index) {   foreach (int unused in data.UnusedIndices()) {     if (index-- <= 0) {       return unused;     }   }   return -1; } public static IEnumerable UnusedIndices(this int[] data) {   for (int i = 0; i < data.Length; i++) {     if (data[i] < 0) {       yield return i;     }   } } public static string PermutationToString(this int[] p, bool oneBased) {   StringBuilder build = new StringBuilder(p.Length * 4);   int width = ((int)Math.Log10(p.Length)) + 2;   build.Append("[");   for (int i = 0; i < p.Length; i++) {     int index = oneBased ? p[i] + 1 : p[i];     build.Append(index.ToString().PadLeft(width));   }   build.Append("]");   return build.ToString(); } #endregion #region Matrix public static string MatrixToString(this double[][] A) {   if (A != null) {     StringBuilder build = new StringBuilder(A.Length * A.Length * 3);     for (int i = 0; i < A.Length; i++) {       if (A[i] != null) {         for (int j = 0; j < A[i].Length; j++) {           build.AppendFormat("{0,4}", A[i][j]);           build.Append(" ");         }         build.AppendLine();       }     }     return build.ToString();   }   return null; } public static double[][] NewMatrix(int m, int n) {   double[][] M = new double[m][];   for (int i = 0; i < M.Length; i++) {     M[i] = new double[n];   }   return M; } `
`public static double[][] NewMatrix(int n) {   return NewMatrix(n, n); } #endregion } `

## 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 readonly Qap _qap;
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.

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

## 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) {
string content = reader.ReadToEnd();
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.