Local Search
Local search is a technique that takes as an input a feasible solution to a combinatorial optimization problem, and tries to produce a better solution. The idea, at a high level, is to search locally in a neighborhood of the input point for a solution that immediately improves the objective value. If no such solution can be found, our input is locally optimal and is returned. Otherwise, we recursively call our local search algorithm on the superior point found.
More formally, suppose we are given some optimization problem
\[
\begin{aligned}
&\min & f(\mathbf x)\\
&\text{subject to}& \mathbf x &\in S
\end{aligned}
\]
where
and
. In local search, for every
, we must give some
to be the neighborhood of
. Then given some initial feasible solution
, for
, we let
\[
\mathbf x_{n} = \mathop{\rm arg\, min}_{\mathbf x \in N_{\mathbf x_{n-1}}} \left\{ f(\mathbf x)\right\}
\]
To run a local search algorithm you begin with
, and iterate the above process until it converges, i.e.
.
The key to making local search work is how we choose
. We need be able to efficiently optimize over
, otherwise our local search heuristic will be no easier to solve than our original problem. However, if we define
to be too small a set, then we may get stuck in a local optimum that isn't very good.
A few warnings about local search:
- Without special problem structure, local search provides no guarentee on the quality of the solution produced relative to the global optimum
- Generally speaking, we do not have a good upper bound on the number of iterations until local search converges. We have the trivial bounds of
, which may be exponentially large, and (if the objective is guaranteed to be integer),
, which can still often be exponential (assuming the problem has weights which are input in binary).
An example
TODO
Two-Opt for TSP
Two Opt is one of the simplest of many local search techniques for TSP. The idea is very simple: given a tour, choose any two (non-adjacent) edges, and swap the end points to make a new tour. For example, in the figure below, we begin with the tour given by the node ordering 1, 2, 3, 4, 5, 6, 7, 8, and we we pick the edge from 2 to 3 and the edge from 5 to 6 and swap, to produce the node ordering 1, 2, 5, 4, 3, 6, 7, 8.
For a given tour, the neighborhood is the set of tours that can be made with such a swap (if
, there are
). To optimize over all such neighborhoods, you simply enumerate them and calculate their cost relative to the existing tour (this can be done in
as really you only need add and subtract the cost of the four edges you are adjusting). As discussed here, there are a variety of heuristics that can be used to speed up this process.
The origin of the name is as follows. For any tour, we can represent it as a vector
where for each
, we have
if
is in the tour and
otherwise. If
is the set of edges in the current tour being considered by Two-Opt, and
is the set of all feasible tours, then the set of tours considered by Two-Opt this round is exactly
\[
\begin{aligned}
\sum_{e \in E \setminus T} x_e &\leq 2,\\
\mathbf x &\in \mathcal T.
\end{aligned}
\]
More generally, one can consider swapping the end points of up to
edges, giving the local search algorithm
-Opt. However, with a naïve implementation, this will require
time per call on
vertices (times the recursive depth!), the cost prohibitive.
Implementing Two-Opt for User Generated Solutions
We have provided you with a simple Two-Opt implementation, TwoOptSolver
. As a minor optimization, it caches ever tour it has ever checked the neighborhood of, so it can abort the search as soon as it has been determined that it will end up in a local optimum that has already been visited. It has a simple interface
Method Name |
Return Type |
Arguments |
Description |
searchNeighborhoodEdgeList |
Set<E> |
List<E> startTour |
Applies Two-Opt recurssively to startTour, returns a locally optimal tour, or null if an already visited tour is hit. |
searchNeighborhoodEdgeSet |
Set<E> |
Set<E> startTour |
Same as above, but must pay a one time
to convert into the List<E> format. |
We are going to apply Two-Opt to the solutions we generate through Christofides algorithm when we set the MIP start and when we run our Heuristic Callback. First, we add a field and initialize it in the constructor:
private ChristofidesSolver<V,E> christofides;//this was old
private TwoOpt<V,E> twoOpt;//<- add this line!
and change this piece of existing code:
if(options.contains(Option.christofidesApprox)){
System.out.println("Beginning Christofides...");
List<E> christofidesApproxTour = christofides.approximateBestTour(new HashSet<E>());
Set<E> submittedTour = new HashSet<E>(christofidesApproxTour);
double[] edgeVals = inverseEdgesUsed(submittedTour);
cplex.addMIPStart(edgeVariablesAsArray, edgeVals);
}
to read:
if(options.contains(Option.twoOpt)){
twoOpt = new TwoOpt<V,E>(tspInstance);
}
if(options.contains(Option.christofidesApprox)){
System.out.println("Beginning Christofides...");
List<E> christofidesApproxTour = christofides.approximateBestTour(new HashSet<E>());
Set<E> submittedTour;
if(options.contains(Option.twoOpt)){
submittedTour = twoOpt.searchNeighborhoodEdgeList(christofidesApproxTour);
}
else{
submittedTour = new HashSet<E>(christofidesApproxTour);
}
double[] edgeVals = inverseEdgesUsed(submittedTour);
cplex.addMIPStart(edgeVariablesAsArray, edgeVals);
}
Now Two-Opt will be run on the initially generated solution. Note that we do not need to make sure that the output of Two-Opt is non-null, since we cannot hit the same tour twice on our first run. Now try and modify ChristofidesCallback
yourself!
Toggle Solution
Replace the code
List<E> approxTourAsList = christofides.approximateBestTour(suggest);
if(approxTourAsList == null){
return;
}
Set<E> submittedTour = new HashSet<E>(approxTourAsList);
double submittedCost = tspInstance.cost(submittedTour);
System.out.println("heuristic solution cost: " + submittedCost);
double[] suggestedSolution = inverseEdgesUsed(submittedTour);
this.setSolution(edgeVariablesAsArray, suggestedSolution);
from ChristofidesCallback
by
List<E> approxTourAsList = christofides.approximateBestTour(suggest);
if(approxTourAsList == null){
return;
}
Set<E> submittedTour;
if(options.contains(Option.twoOpt)){
submittedTour = twoOpt.searchNeighborhoodEdgeList(approxTourAsList);
if(submittedTour == null){
System.out.println("Two-Opt revisted same solution");
return;
}
}
else{
submittedTour= new HashSet<E>(approxTourAsList);
}
double submittedCost = tspInstance.cost(submittedTour);
System.out.println("heuristic solution cost: " + submittedCost);
double[] suggestedSolution = inverseEdgesUsed(submittedTour);
this.setSolution(edgeVariablesAsArray, suggestedSolution);
Profiling on TSPLIB
Small instances
Using QuadraticWaiting as our BackOff functions for both user cut generation and heuristic solution generation and rerunning our code with two-opt in place, we get the following results:
We see that Two-Opt gives a significantly better initial solution. Running times and number of nodes visited are also reduced.
Large instances
We see similar results for large instances. However, as the majority of the solution time is spent closing the final few percentage points, and our heuristics don't help much in this regime, there still seems to be significant room for improvement.
Perhaps more impressively, we now can get a near optimal solution very quickly. On d493, in less than one minute and on only 5 nodes, we found a solution with 0.77% of optimal, and on d657, in 90 seconds and 142 nodes, we get within 2.02%, which is often good enough. However, on large problems, like d657, we still have no way of reaching optimality in a reasonable amount of time.
Complete Source
If all else fails, your code should like the code below.
Toggle Complete Source
public class TspIpSolver<V,E> {
public static enum Option{
lazy,userCut, randomizedUserCut, christofidesApprox, christofidesHeuristic,twoOpt,incumbent;
}
private EnumSet<Option> options;
private IloCplex cplex;
private TspInstance<V,E> tspInstance;
private final ImmutableBiMap<E,IloIntVar> edgeVariables;
private ImmutableSet<E> edgesInOpt;
private double optVal;
private IloIntVar[] edgeVariablesAsArray;
private ChristofidesSolver<V,E> christofides;
private TwoOpt<V,E> twoOpt;
public TspIpSolver(TspInstance<V,E> tspInstance) throws IloException{
this(tspInstance,EnumSet.of(Option.lazy, Option.userCut, Option.christofidesApprox, Option.christofidesHeuristic));
}
public TspIpSolver(TspInstance<V,E> tspInstance, EnumSet<Option> options) throws IloException{
this.options = options;
this.tspInstance = tspInstance;
this.cplex = new IloCplex();
UndirectedGraph<V,E> graph = tspInstance.getGraph(); //for convenience, we will be using this a lot
this.edgeVariables = Util.makeBinaryVariables(cplex, graph.getEdges());
edgeVariablesAsArray = edgeVariables.inverse().keySet().toArray(new IloIntVar[]{});
//the degree constraints
for(V vertex: graph.getVertices()){
cplex.addEq(Util.integerSum(cplex, edgeVariables,
graph.getIncidentEdges(vertex)), 2);
}
//the objective
cplex.addMinimize(Util.integerSum(
cplex, edgeVariables, graph.getEdges(),tspInstance.getEdgeWeights()));
if(options.contains(Option.lazy)){
cplex.use(new IntegerCutCallback());
}
if(options.contains(Option.userCut)){
cplex.use(new MinCutCallback());
if(options.contains(Option.randomizedUserCut)){
System.err.println("Warning, userCut and randomizedUserCut both selected, ignoring randomizedUserCut");
}
}
else if(options.contains(Option.randomizedUserCut)){
cplex.use(new RandomMinCutCallback());
}
if(options.contains(Option.christofidesApprox)||options.contains(Option.christofidesHeuristic)){
christofides = new ChristofidesSolver<V,E>(graph,tspInstance.getEdgeWeights());
}
if(options.contains(Option.twoOpt)){
twoOpt = new TwoOpt<V,E>(tspInstance);
}
if(options.contains(Option.christofidesApprox)){
System.out.println("Beginning Christofides...");
List<E> christofidesApproxTour = christofides.approximateBestTour(new HashSet<E>());
Set<E> submittedTour;
if(options.contains(Option.twoOpt)){
submittedTour = twoOpt.searchNeighborhoodEdgeList(christofidesApproxTour);
}
else{
submittedTour = new HashSet<E>(christofidesApproxTour);
}
double[] edgeVals = inverseEdgesUsed(submittedTour);
cplex.addMIPStart(edgeVariablesAsArray, edgeVals);
}
if(options.contains(Option.christofidesHeuristic)){
cplex.use(new ChristofidesCallback());
}
}
public void solve() throws IloException{
if(!cplex.solve()){
throw new RuntimeException();
}
optVal = cplex.getObjValue();
edgesInOpt = ImmutableSet.copyOf(edgesUsed(cplex.getValues(edgeVariablesAsArray)));
cplex.end();
}
public ImmutableSet<E> getEdgesInOpt(){
return edgesInOpt;
}
public double getOptVal(){
return optVal;
}
/**
* assumes edgeVarVals[i] in {0,1}, up to some tolerance.
* @param edgeVarVals
* @return the edges where the variable takes the value one.
*/
private Set<E> edgesUsed(double[] edgeVarVals){
Set<E> ans = new HashSet<E>();
for(int e = 0; e < edgeVarVals.length; e++){
if(Util.doubleToBoolean(edgeVarVals[e])){
ans.add(edgeVariables.inverse().get(edgeVariablesAsArray[e]));
}
}
return ans;
}
/**
* Given a subset of the edges, produces a binary vector indicating which edges were included using the same ordering as edgeVariablesAsArray
* @param edgesUsed a subset of the set of edges in the grpah
* @return an array of zeros and ones where the ith entry corresponds to the ith edge variable from edgeVariablesAsArray
*/
private double[] inverseEdgesUsed(Set<E> edgesUsed){
double[] edgeVals = new double[this.edgeVariablesAsArray.length];
for(int i =0; i < edgeVals.length; i++){
edgeVals[i] = edgesUsed.contains(edgeVariables.inverse().get(edgeVariablesAsArray[i])) ? 1 : 0;
}
return edgeVals;
}
private Map<E,Double> getNonZeroEdgeWeights(double[] edgeValues){
Map<E,Double> edgeWeights = new HashMap<E,Double>();
for(int i = 0; i < edgeValues.length; i++){
if(edgeValues[i] > Util.epsilon){
edgeWeights.put(edgeVariables.inverse().get(edgeVariablesAsArray[i]), edgeValues[i]);
}
}
return edgeWeights;
}
/**
* assumes edgeVarVals[i] in [0,1]. Different than edgesUsed because an exception will not be thrown in
* 0 < edgeVarVals[i] < 1.
* @param edgeVarVals
* @return the edges where the variable takes the value one, up to some tolerance.
*/
private Set<E> edgesAtOne(double[] edgeVarVals){
Set<E> ans = new HashSet<E>();
for(int e = 0; e < edgeVarVals.length; e++){
if(edgeVarVals[e] >= 1-Util.epsilon){
ans.add(edgeVariables.inverse().get(edgeVariablesAsArray[e]));
}
}
return ans;
}
private class IntegerCutCallback extends LazyConstraintCallback{
private double minCutVal;
public IntegerCutCallback(){
minCutVal = 1.99;
}
@Override
protected void main() throws IloException {
Set<E> edgesUsed = edgesUsed(this.getValues(edgeVariablesAsArray));
Set<Set<V>> connectedComponents = tspInstance.getConnectedComponents(edgesUsed);
if(connectedComponents.size() > 1){
for(Set<V> connectedComponent: connectedComponents){
this.add(cplex.ge(Util.integerSum(cplex, edgeVariables, tspInstance.cutEdges(connectedComponent)),2));
}
}
}
}
private class MinCutCallback extends UserCutCallback{
private BackOffFunction backOff;
private double minCutVal;
public MinCutCallback(){
minCutVal = 1.99;
this.backOff = new BackOffFunction.QuadraticWaiting();
}
@Override
protected void main() throws IloException {
//do not attempt to add any cutset constraints until CPLEX is done adding constraints
if(!this.isAfterCutLoop()){
return;
}
double[] edgeVals = this.getValues(edgeVariablesAsArray);
Map<E,Double> edgeWeights = getNonZeroEdgeWeights(edgeVals);
UndirectedGraph<V,E> graph = tspInstance.getGraph();
MinCutSolver<V,E> minCutSolver = new MinCutSolver<V,E>(graph,edgeWeights);
Set<Set<E>> cutSets = new HashSet<Set<E>>();
Set<Set<V>> connectedComponents = minCutSolver.checkConnectedComponents();
if(connectedComponents.size() > 1){
for(Set<V> connectedComponent: connectedComponents){
cutSets.add(tspInstance.cutEdges(connectedComponent));
}
System.out.println("trivial cuts exist, found " + cutSets.size());
}
else if(this.getNnodes64() == 0 || this.backOff.makeAttempt()){
V source = graph.getVertices().iterator().next();
double best = Double.MAX_VALUE;
System.out.print("looking for cut... ");
for(V sink : graph.getVertices()){
if(sink != source){
Set<V> cutNodes = minCutSolver.computeMinCut(source, sink);
Set<E> cutEdges = tspInstance.cutEdges(cutNodes);
double cutVal = Util.calcSum(cutEdges,edgeWeights);
best = Math.min(best, cutVal);
if(cutVal < minCutVal){
cutSets.add(cutEdges);
}
}
}
System.out.println("found " + cutSets.size() + " cuts, best " + best);
}
for(Set<E> cutEdges: cutSets){
this.add(cplex.ge(Util.integerSum(cplex, edgeVariables, cutEdges), 2));
}
}
}
private class RandomMinCutCallback extends UserCutCallback{
private int numAttempts;
private Random rand;
private BackOffFunction backOff;
private double minCutVal;
public RandomMinCutCallback(){
numAttempts = 20;
rand = new Random();
backOff = new BackOffFunction.QuadraticWaiting();
minCutVal = 1.99;
}
@Override
protected void main() throws IloException {
//do not attempt to add any cutset constraints until CPLEX is done adding constraints
if(!this.isAfterCutLoop()){
return;
}
double[] edgeVals = this.getValues(edgeVariablesAsArray);
Map<E,Double> edgeWeights = getNonZeroEdgeWeights(edgeVals);
UndirectedGraph<V,E> graph = tspInstance.getGraph();
MinCutSolver<V,E> minCutSolver = new MinCutSolver<V,E>(graph,edgeWeights);
Set<Set<E>> cutSets = new HashSet<Set<E>>();
Set<Set<V>> connectedComponents = minCutSolver.checkConnectedComponents();
if(connectedComponents.size() > 1){
for(Set<V> connectedComponent: connectedComponents){
cutSets.add(tspInstance.cutEdges(connectedComponent));
}
System.out.println("trivial cuts exist, found " + cutSets.size());
}
else if(this.getNnodes64() == 0 || this.backOff.makeAttempt()){
double best = Double.MAX_VALUE;
System.out.print("looking for cut... ");
List<V> nodes = new ArrayList<V>(graph.getVertices());
for(int i = 0; i <numAttempts; i++){
V source = nodes.get(rand.nextInt(nodes.size()));
V sink = nodes.get(rand.nextInt(nodes.size()));
if(source != sink){
Set<V> cutNodes = minCutSolver.computeMinCut(source, sink);
Set<E> cutEdges = tspInstance.cutEdges(cutNodes);
double cutVal = Util.calcSum(cutEdges,edgeWeights);
best = Math.min(best, cutVal);
if(cutVal < minCutVal){
cutSets.add(cutEdges);
}
}
}
System.out.println("found " + cutSets.size() + " cuts, best " + best);
}
for(Set<E> cutEdges: cutSets){
this.add(cplex.ge(Util.integerSum(cplex, edgeVariables, cutEdges), 2));
}
}
}
private class ChristofidesCallback extends HeuristicCallback{
private Set<Set<E>> previouslySuggested;
private BackOffFunction backoff;
public ChristofidesCallback(){
previouslySuggested = new HashSet<Set<E>>();
backoff = new BackOffFunction.QuadraticWaiting();
}
@Override
protected void main() throws IloException {
if(this.getMIPRelativeGap()< 0.01|| this.getNnodes64() == 0){
return;
}
double[] edgeVals = this.getValues(edgeVariablesAsArray);
Set<E> suggest = edgesAtOne(edgeVals);
if(previouslySuggested.contains(suggest)){
return;
}
else{
previouslySuggested.add(suggest);
}
if(!backoff.makeAttempt()){
return;
}
List<E> approxTourAsList = christofides.approximateBestTour(suggest);
if(approxTourAsList == null){
return;
}
Set<E> submittedTour;
if(options.contains(Option.twoOpt)){
submittedTour = twoOpt.searchNeighborhoodEdgeList(approxTourAsList);
if(submittedTour == null){
System.out.println("Two-Opt revisted same solution");
return;
}
}
else{
submittedTour= new HashSet<E>(approxTourAsList);
}
double submittedCost = tspInstance.cost(submittedTour);
System.out.println("heuristic solution cost: " + submittedCost);
double[] suggestedSolution = inverseEdgesUsed(submittedTour);
this.setSolution(edgeVariablesAsArray, suggestedSolution);
}
}
}