4-opt* – general idea, part 2

The second crucial part of 4-opt* idea is to examine only a specific cases of potentially improving moves.

Any k-opt move that improves a tour move must fulfill the condition:

addLength < delLength

Sum of length of links added must be less than sum of links removed from tour. It is therefore necessary (though insufficient) that at least one of added links is shorter than one of removed links. Lengths of all four links removed in 4-opt* procedure are (see schemes below):

delLen_1 = distance(tour[i], tour[(i+1) mod N])
delLen_2 = distance(tour[j], tour[(j+1) mod N])
delLen_3 = distance(tour[k], tour[(k+1) mod N])
delLen_4 = distance((k+1) mod N], tour[(k+2) mod N])

We take the length of the longest of them:

delLenMax = max(delLen_1, max (delLen_2, max(delLen_3, delLen_4)))

We determine that to reconnect the cities we always add the link to city in position k+1 either from city in position i+1 or from city in j and that we take the shorter of these two links:

addLen_1 = distance(tour[(i+1) mod N], tour[(k+1) mod N])
addLen_2 = distance(tour[(j) mod N], tour[(k+1) mod N])

addLenMin = min(addLen1, addLen2)

Possible cases of reconnection and gain from them are tested only if:

addLenMin < delLenMax

This test guarantees that at least this added link is shorter than at least one of the links removed. If the condition is not met, then cases of 4-opt* move are not examined. This way we examine only one of two classes of moves:

Fixed links in 4-opt* moves

Class A	Class B

So for given i, j, k we consider only one class of reconnections: A or B. Either class have only 8 cases of reconnection (see below), as opposed to 48 cases in standard 4-opt. Moreover, since the scheme of class B is symmetrical to the scheme of class A, then cases of reconnection for class B will be also symmetrical to cases for class A.

4-opt* move: reconnection cases of class A

2-opt move
3-opt moves
4-opt moves

4-opt* – general idea, part 1

Since full 4-opt is rather complex and time comsuming, Renaud, Boctor and Laporte (1996) proposed a tour improving procedure caled 4-opt*, based on using a specific subset of 4-opt moves. Its main goal is to extend searching for improving moves beyond 2-opt moves, but keep this searching fast.

Elements of 4-opt move	Elements of 4-opt* move

As we can see, the first difference between standard 4-opt and 4-opt* is that while the scheme of general 4-opt move involves four non-overlapping segments with eight endpoints, the scheme of 4-opt* involves only three segments and seven points; fourth segment is reduced to one point between endpoints of two other segments. The layout for standard 4-opt requires four parameters (i, j, k, m), but for 4-opt* it is only three parameters (i, j, k). This allows us to reduce number of corresponding nested loops in searching for an improvement, each of O(N). (There is additional advantage: reduction of possible cases of reconnection. In 4-opt move they result from its four parameters, but 4-opt* cases result from only three parameters. Therefore number of cases (types of reconnection) of a move in 4-opt* is reduced compared to the number of cases in standard 4-opt. For example double bridge cannot be made in 4-opt*).

Moreover, value of j is obtained from: j = i + u, where u <= w and w is some arbitrary number chosen by the developer or a user and usually relatively small compared to N (for example: 10, 15 or 20). Therefore for given u the procedure considers only O(N²) choices, changing values of i and k. This way we consider links at the ends of two not overlaping chains: from position i to i+u+1, and from position k to k+2.

The 4-opt* algorithm starts with u=1 and considers all choices of such two chains of cities. If an improvement is found then it is immediately applied and searching is repeated. When no improvement is found then u is incremented by 1 and the same process is repeated. When the process of improvimg for u==w is finished, then u is reset to 1 and the same procedure is applied once again, until there is no improvement in such full round.

improved = false
while not improved:
  u = 1
  improved = false
  while u <= w:
    block two_loops:
      for i in ... :
        j = i + u
        ...
        for k in ... :
          ...
          # additional speed-up
          ...
          for case_of_move in ... :
            gain = Gain_From_4optAsterisk(i, j, k, case_of_move)
            if gain > 0:
              Make_4optAsterisk_Move(i, j, k, case_of_move)
              improved = true
              break two_loops
    #end_block two_loops
    u = u + 1
  #end_while u loop

4-opt and Double Bridge

Replacing 3-optic moves in the local optimization with 4-opt moves increases the time complexity the algorithm from O(N³) to O(N⁴), or, in the case of using lists of candidates (neighbors), from O(N * m²) to O (N * m³), where m is number of neighbors in the list.

However, the problem does not lie only in much greater time of execution. Applying 4-opt moves gives tours only slightly shorter than those from 3-optimization. Moreover, while for 3-opt there are 8 types (cases) of move, for 4-opt this number increases to 48, so considering all of them becomes much more complicated. As in 3-opt, one of them gives the tour identical with the tour before the move. Analogically to 3-opt, many of them are equal to a single 2-opt or a single 3-opt move; there are only 25 pure 4-opt moves. But something new appears in 4-opt moves: most known 4-opt move, called double bridge move.

The simplicity and symmetry of the double bridge scheme make it particularly easy to recognize and remember:

Double bridge, as the simplest non-sequential move (idea of so called sequenctial moves is described later), plays important role even in algorithms that do not use full 4-optimization. By definition given tour after full 3-optimization cannot be further improved by 3-opt, but it does not mean that there are no similar but shorter tours. Instead of trying to search or build such tours from scratch, we often use 3-optimized tour and perturb it to obtain a new good initial tour, that we can then optimize, perhaps getting better result. Double bridge is used as simple and effective perturbation.

Double bridge and other non-sequential 4-opt moves are discussed in detail later.

Dynasearch – the idea of set of independent moves

So far we have been dealing with the simplest methods in which a single optimization step is a single improving move of given type (for example 2-opt, 2.5-opt, 3-opt). In more advanced methods, a single step is an improving set of moves. An example of such algorithms is the dynasearch method.

In the dynasearch approach we use dynamic programming algorithm to look for the best set consisting of independent improving moves. Two moves are called independent when, regardless of the order in which they are applied, the resulting tour is the same. For example two 2-opt moves are independent when they do not overlap (their ranges do not overlap):

independent moves	non-independent moves

This means that none of the endpoints of the second move cannot be between the endpoints of the first move.

So for pairs in which the two positions are in ascending order:

0 <= i1 < ... < j1 <= N-1, where i1+1 != j1 and j1+1 != i1

and

0 <= i2 < ... < j2 <= N-1, where i2+1 != j2 and j2+1 != i2

the two 2-opt moves are independent if:

j1 < i2 or j2 < i1

2-opt dynasearch on average is faster than traditional best improvement approach and gives slightly better results. Unfortunatelly, running times are far longer than those for 3-opt with neighbor lists, which gives better results and is much simpler to implement.

Beware

But what when the endpoints of a move are not in ascending order of indices? Let us consider the following:

If the first move is defined by the pair (3, N-2) then the moves are not independent, because all links involved in second move would be affected by changes made by the first move (case A). But if the first move is defined by the pair (N-2, 3) then the two moves are independent (case B):

case A reversing the segment without the city on position 0	case B reversing the segment with the city on position 0

How is it possible that in one case the moves are independent and in the other they are not, when they seem to produce the same tour, differing only in orientation?

The main difference and the problem with case A is that the first move by reversing a segment simultaneously changes the tour positions (indices) of the cities it contains, so after this move i2, i2+1, j2, j2+1 no longer point at the same cities and links they pointed at before the move. So if we found the second move and set up the values for i2, j2 independently from any other move, not assuming that the move (i1, j1) would be apllied before, then after this first move these values are wrong. This does not happen in case B.

But on the other hand, unfortunatelly in case B we cannot apply the rule:

j1 < i2 or j2 < i1

to check independence of the two moves, because the inequity i1 < j1 is not satisfied.

More general condition for independent 2-opt moves is:

Two moves are independend when every city involved in one move comes before every city involved in second move or vice versa.

Alternative for DLB

Instead of don't look bits speedup we can use an array indexed by city numbers, in which each element contains length of tour after improving from given city. The algorithm using has the following scheme:

var
    tourLen: tourLen  # current length of tour
    tourLenArr: array[City_Number, tourLen]

  # calculate initial length of tour
  tourLen = distance(tour[N-1], tour[0])
  for pos in 0 .. N-2:
    tourLen = tourLen + distance(tour[pos], tour[pos+1])

  # initialise array of tour lengths with some arbitrary
  # value different from current value of 'tourLen'
  for city in 0 .. N-1:
    tourLenArr[city] = 1 # 

  # main loop
  for pos in 0 .. N-1:
    baseCity = tour[pos]
    if (tourLen == tourLenArr[baseCity]):
       continue

    improved = Improve(baseCity)  # MUST update 'tourLen' value!

    tourLenArr[baseCity] = tourLen
  #end_for

Stop search after encountering the successor

Less popular speed-up is using search cut after encountering a successor of given city when iterating on its neighbors. The idea is based on the observation that lead to 2-opt fixed radius search.

We consider neighbours of city X1 to find a good candidate for Y1, the first of endpoints of second link to remove. If we make a 2-opt move, then Y1 would be the new tour successor of X1, so it should not be in distance from X1 greater then is X2, the current successor. Otherwise, the new link would be longer than the current one. Since we consider neighbors of X1 in ascending order of their distance from this city, encountering X2 during search means that all subsequent neighbors would give longer links to X1.

...
  for neighbor_1 in 0 .. numberOfNeigbors-1:
    Y2 = neighbor[X1][neighbor_1]
    j  = (position[Y2] + N - 1) mod N
    Y1 = tour[j]

    if (Y1 == X2):  # speed-up
      break
...

Analogous condition can be used in 3-opt, also for Z1, although justification in the case of 3-opt is not so strong, even when forcing only Y1 != X2.

Note that this condition rejects some valid 3-opt moves, Node Shifts that are used in 2.5-opt:

This speed-up is very effective and can be used instead of popular don’t look bits solution.

The full code for simple 3-opt and the speed-up in one loop only is as follows:

proc One_City_3_Opt_NC(tour: var Tour_Array;
                       basePos: Tour_Index;
                       neighbor: Neighbor_Matrix;
                       position: var City_Position_In_Tour;
                       numberOfNeigbors: int): bool =
  type
    Move3 = object
      i, j, k: Tour_Index
      optCase: Reconnection_3_optCase
      gain: Length_Gain
  var
    improved: bool
    i, j, k: Tour_Index
    X1, X2, Y1, Y2, Z1, Z2: City_Number
    k_6, k_7:  Tour_Index
    Z1_6, Z2_6, Z1_7, Z2_7: City_Number
    gainExpected: Length_Gain
    goodMove: Move3
  
  improved = false
  goodMove.gain = 0
  block three_loops:
    
    for direction in [forward, backward]:
      if direction == forward:
        i  = basePos
      else:
        i  = (N + basePos - 1) mod N
      X1 = tour[i]
      X2 = tour[(i+1) mod N]
 
      for neighbor_1 in 0 .. numberOfNeigbors-1:
        # new edges in optCase=6: *X1-Y2*, Y1-Z1, X2-Z2
        # new edges in optCase=7: *X1-Y2*, X2-Z1, Y1-Z2
        Y2 = neighbor[X1][neighbor_1]
        j  = (position[Y2] + N - 1) mod N
        Y1 = tour[j]

        if (Y1 != X1) and (Y1 != X2):
          gainExpected = Gain_From_2_Opt(X1, X2, Y1, Y2)
          if gainExpected > goodMove.gain:
            improved = true
            goodMove = Move3(i: i,  j: j,  k: j,
                             optCase: opt3_case_1, gain: gainExpected)

        if (Y1 == X2):  # NOTE
          break

        for neighbor_2 in 0 .. numberOfNeigbors-1:
          # new edges in optCase=6: X1-Y2, *Y1-Z1*, X2-Z2
          Z1_6 = neighbor[Y1][neighbor_2]
          k_6  = position[Z1_6]
          Z2_6 = tour[(k_6 + 1) mod N]
          # new edges in optCase=7: X1-Y2, X2-Z1, *Y1-Z2*
          Z2_7 = neighbor[Y1][neighbor_2]
          k_7  = (position[Z2_7] + N - 1) mod N
          Z1_7 = tour[k_7]
          
          if Between(i, j, k_6):
            gainExpected = Gain_From_3_Opt(X1, X2, Y1, Y2, Z1_6, Z2_6,
                                           opt3_case_6)
            if gainExpected > goodMove.gain:
              improved = true
              goodMove = Move3(i: i,  j: j,  k: k_6,
                               optCase: opt3_case_6, gain: gainExpected)

          if Between(i, j, k_7):
            gainExpected = Gain_From_3_Opt(X1, X2, Y1, Y2, Z1_7, Z2_7,
                                           opt3_case_7)
            if gainExpected > goodMove.gain:
              improved = true
              goodMove = Move3(i: i,  j: j,  k: k_7,
                               optCase: opt3_case_7, gain: gainExpected)

  if improved:
    Make_3_Opt_Move(tour,
                    goodMove.i,  goodMove.j,  goodMove.k,
                    goodMove.optCase,  position)
  result = improved

proc LS_3_Opt_NC(tour: var Tour_Array;
                 neighbor: Neighbor_Matrix;
                 neighborListLen: int = N-1) =
  var
    locallyOptimal: bool = false
    improved: bool
    baseCity: City_Number
    position: City_Position_In_Tour

  let numberOfNeigbors = min(neighborListLen, len(neighbor[0]))
  position = Create_Position_In_Tour(tour)
  
  while not locallyOptimal:
    locallyOptimal = true
    for basePos in 0 .. N-1:
      baseCity = tour[basePos]
      improved = One_City_3_Opt_NC(tour, basePos, neighbor, position,
                                   numberOfNeigbors)
      if improved:
        locallyOptimal = false

Some simple statistics

Relative time, best length, average length and worst length for all tours created by NN heuristic and improved by some 2-opt and 3-opt algorithms. Random planar problems, N=200.

Problem #	Method	Time	Best	Avg	Worst
1	NN	100	100.00	106.52	114.46
	NN + 2opt	4721	88.17	90.89	94.52
	NN + 2opt Take Best	5119	87.53	89.09	91.25
	NN + 3opt-ND(24)	5851	86.44	87.94	89.65
	NN + 3opt-NC(24)	3589	85.76	87.89	89.86
2	NN	100	100.00	104.27	110.57
	NN + 2opt	3623	90.30	92.45	95.40
	NN + 2opt Take Best	3889	89.07	90.67	92.78
	NN + 3opt-ND(24)	5551	88.31	89.59	92.20
	NN + 3opt-NC(24)	2959	88.36	89.46	91.41
3	NN	100	100.00	110.24	122.10
	NN + 2opt	3839	91.03	94.61	97.58
	NN + 2opt Take Best	4450	90.34	91.70	94.47
	NN + 3opt-ND(24)	5808	88.72	90.63	92.96
	NN + 3opt-NC(24)	3634	88.82	90.71	92.56
4	NN	100	100.00	106.87	114.35
	NN + 2opt	3955	87.53	90.64	94.81
	NN + 2opt Take Best	4521	87.10	89.33	91.84
	NN + 3opt-ND(24)	5785	85.53	87.35	90.17
	NN + 3opt-NC(24)	3423	85.99	87.44	89.56
5	NN	100	100.00	106.86	113.89
	NN + 2opt	3723	90.16	93.61	96.94
	NN + 2opt Take Best	4421	89.79	92.18	95.17
	NN + 3opt-ND(24)	5982	88.61	90.61	92.40
	NN + 3opt-NC(24)	3691	88.57	90.66	94.46
6	NN	100	100.00	106.80	114.76
	NN + 2opt	4853	85.03	88.16	92.81
	NN + 2opt Take Best	5087	85.07	86.95	91.02
	NN + 3opt-ND(24)	5882	83.19	84.86	87.59
	NN + 3opt-NC(24)	3823	83.69	85.11	88.18
7	NN	100	100.00	104.23	110.78
	NN + 2opt	4487	86.00	88.67	92.58
	NN + 2opt Take Best	5287	84.92	86.93	88.92
	NN + 3opt-ND(24)	5882	83.92	85.59	88.45
	NN + 3opt-NC(24)	3923	83.63	85.60	88.26
8	NN	100	100.00	109.86	116.72
	NN + 2opt	5285	86.01	89.92	93.72
	NN + 2opt Take Best	5951	86.04	89.06	91.54
	NN + 3opt-ND(24)	5985	84.57	86.88	90.18
	NN + 3opt-NC(24)	3523	84.83	86.74	89.65
9	NN	100	100.00	104.99	110.15
	NN + 2opt	4321	87.50	90.42	93.91
	NN + 2opt Take Best	4721	87.01	88.96	90.70
	NN + 3opt-ND(24)	5951	85.97	87.41	89.39
	NN + 3opt-NC(24)	4055	85.51	87.37	89.05
10	NN	100	100.00	106.83	115.85
	NN + 2opt	4321	90.29	92.68	96.71
	NN + 2opt Take Best	4721	89.26	90.50	92.66
	NN + 3opt-ND(24)	5717	89.08	90.35	92.09
	NN + 3opt-NC(24)	3823	89.15	90.24	91.92

3-opt with min-max lenght of link – example

proc One_City_3_Opt_NDM(tour: var Tour_Array;
                        basePos: Tour_Index;
                        neighbor: Neighbor_Matrix;
                        position: var City_Position_In_Tour;
                        numberOfNeigbors: int;
                        DontLook: var DLB_Array;
                        minLenInSet: Length;
                        maxLenInTour: var Length): bool =
  type
    Move3 = object
      i, j, k: Tour_Index
      optCase: Reconnection_3_optCase
      gain: Length_Gain
  var
    improved: bool
    i, j, k: Tour_Index
    ii, jj, kk: Tour_Index
    X1, X2, Y1, Y2, Z1, Z2: City_Number
    A1, A2, B1, B2, C1, C2: City_Number
    gainExpected: Length_Gain
    goodMove: Move3
    d1, d2: Length
  
  improved = false
  goodMove.gain = 0
  block three_loops:
    
    for direction in [forward, backward]:
      if direction == forward:
        i  = basePos
      else:
        i  = (N + basePos - 1) mod N
      X1 = tour[i]
      X2 = tour[(i+1) mod N]
 
      for neighbor_1 in 0 .. numberOfNeigbors-1:
        Y2 = neighbor[X1][neighbor_1]
        j  = (position[Y2] + N - 1) mod N
        Y1 = tour[j]

        if (Y1 != X1) and (Y1 != X2):
          gainExpected = Gain_From_2_Opt(X1, X2, Y1, Y2)
          if gainExpected > goodMove.gain:
            improved = true
            goodMove = Move3(i: i,  j: j,  k: j,
                             optCase: opt3_case_1, gain: gainExpected)
        else:
          continue
          
        d2 = distance(X1, Y1) + minLenInSet
        if d2 + minLenInSet > distance(X1, X2) + 2 * maxLenInTour:
           break
        d1 = distance(X1, X2) + distance(Y1, Y2) + maxLenInTour
        if d2 + minLenInSet > d1:
           continue

        for neighbor_2 in 0 .. numberOfNeigbors-1:
          Z2 = neighbor[X1][neighbor_2]
          k  = (position[Z2] + N - 1) mod N
          Z1 = tour[k]

          if (Z2 == Y2) or (Z1 == Y2) or
             (Z1 == X1) or (Z1 == X2) or
             (Y1 == X1) or (Y1 == X2):
              continue

          if d2 + distance(Y1, Z1) > d1:
            break

          if Between(i, j, k):
            jj = j
            kk = k
            (A1, A2, B1, B2, C1, C2) = (X1, X2, Y1, Y2, Z1, Z2)
          else:
            jj = k
            kk = j
            (A1, A2, B1, B2, C1, C2) = (X1, X2, Z1, Z2, Y1, Y2)

          for optCase in opt3_case_4 .. opt3_case_7:
            gainExpected = Gain_From_3_Opt(A1, A2, B1, B2, C1, C2,
                                           optCase)
            if gainExpected > goodMove.gain:
              improved = true
              goodMove = Move3(i: i,  j: jj,  k: kk,
                               optCase: optCase, gain: gainExpected)

  if improved:
    A1 = tour[goodMove.i]
    A2 = tour[(goodMove.i + 1) mod N]
    B1 = tour[goodMove.j]
    B2 = tour[(goodMove.j + 1) mod N]
    C1 = tour[goodMove.k]
    C2 = tour[(goodMove.k + 1) mod N]
    Set_DLB_off(DontLook, [A1, A2, B1, B2, C1, C2])
    Make_3_Opt_Move(tour,
                    goodMove.i,  goodMove.j,  goodMove.k,
                    goodMove.optCase,  position)

    case goodMove.optCase
    of opt3_case_1:
      maxLenInTour = max( maxLenInTour, max(distance(A1,B1),
                          distance(A2,B2)) )
    of opt3_case_4:
      maxLenInTour = max( maxLenInTour, max(distance(A1,B1),
                          max(distance(A2,C1),
                              distance(B2,C2)) ) )
    of opt3_case_5:
      maxLenInTour = max( maxLenInTour, max(distance(A1,C1),
                          max(distance(A2,B2),
                              distance(B1,C2)) ) )
    of opt3_case_6:
      maxLenInTour = max( maxLenInTour, max(distance(A1,B2),
                          max(distance(A2,C2),
                              distance(B1,C1)) ) )
    of opt3_case_7:
      maxLenInTour = max( maxLenInTour, max(distance(A1,B2),
                          max(distance(A2,C1),
                              distance(B1,C2)) ) )
    else:
      echo "Should not happen: ", goodMove.optCase
      maxLenInTour = distance(tour[0], tour[N-1])
      for pos in 0 .. N-2:
        maxLenInTour = max(maxLenInTour,
                           distance(tour[pos], tour[pos+1]))
    #end goodMove.optCase

  result = improved

proc LS_3_Opt_NDM(tour: var Tour_Array;
                  neighbor: Neighbor_Matrix;
                  neighborListLen: int = N-1) =
  ## Optimizes the given tour using 3-opt with neighbor lists
  ## Don't Look Bits (DLB) and min-max Links
  var
    locallyOptimal: bool = false
    improved: bool
    DontLook: DLB_Array
    baseCity: City_Number
    position: City_Position_In_Tour
    maxLenInTour: Length

  let numberOfNeigbors = min(neighborListLen, len(neighbor[0]))
  position = Create_Position_In_Tour(tour)

#[
  # NOTE: minLenInSet is constant for the dataset
  minLenInSet = high(Length)
  for city_A in 0 .. N-1:
    for city_B in (city_A+1) .. N-1:
      minLenInSet = min(minLenInSet, distance(city_A, city_B))
]#

  maxLenInTour = distance(tour[0], tour[N-1])
  for pos in 0 .. N-2:
    maxLenInTour = max(maxLenInTour, distance(tour[pos], tour[pos+1]))

  Set_DLB_off(DontLook, tour)
  
  while not locallyOptimal:
    locallyOptimal = true
    for basePos in 0 .. N-1:
      baseCity = tour[basePos]
      if isDLB_on(DontLook, baseCity):
        continue
      improved = One_City_3_Opt_NDM(tour, basePos, neighbor, position,
                                    numberOfNeigbors, DontLook,
                                    minLenInSet, maxLenInTour)
      if not improved:
        Set_DLB_on(DontLook, baseCity)
      else:
        locallyOptimal = false

Minimum and maximum length of link as a speed-up

Any 2-opt or 3-opt move that improves a tour move must fulfill the condition:

addLength < delLength

Sum of length of links added must be less than sum of links removed from tour.

In the case of 2-opt the inquity can be written as:

$$d(X1, Y1) + d(X2, Y2) \lt d(X1, X2) + d(Y1, Y2)$$

We can speed-up the search for improving moves by avoiding considering all elements of this inequity.

One method of doing this is fixed radius search. We search for Y1 (or for Y2) inside a the circle centered at X1 and radius to equal the distance between X1 and X2. Moves with both Y's outside the circle do not improve the tour. Unfortunatelly, this method requires symmetrical searching, which could be inconvenient, especiallly when used with some other speed-ups.

However still: in an improving move one or two longer links are removed and then one or two shorter links are added. Suppose that we know the two distances:

$$d(X1, X2), d(X1, Y1)$$

then the inquity has the form:

distance(X1, Y1) + link_added_len < distance(X1, X2) + link_deleted_len

Situation is the better, the bigger link_deleted_len is and the smaller link_added_len is. The maximum value for link_deleted_len is equal to maxLenInTour, which is the length of the longest link in current tour. The mininum value for link_added_len is equal to minLenInSet, which is the length of the shortest of all possible links beetwen cities in the problem (for simplicity we do not exclude links already in current tour, this way minLenInSet is a constant value for the dataset).

var
  minLenInSet, maxLenInTour: Length

  # NOTE: minLenInSet is constant for the dataset
  minLenInSet = high(Length)
  for city_A in 0 .. N-1:
    for city_B in (city_A+1) .. N-1:
      minLenInSet = min(minLenInSet, distance(city_A, city_B))

  maxLenInTour = distance(tour[0], tour[N-1])
  for pos in 0 .. N-2:
    maxLenInTour = max(maxLenInTour, distance(tour[pos], tour[pos+1]))

So if for some Y1 we have got:

distance(X1, Y1) + minLenInSet > distance(X1, X2) + maxLenInTour

then there is no possibility that any 2-opt move with X1, X2 and Y1 could improve the tour, no matter which of the two Y1's tour neighbours we would like to choose for Y2.

Similarilly, we can ommit candidates for 3-opt moves. Consider for example:

distance(X1, Y1) + distance(Y1, Z1) + minLenInSet >
    distance(X1, X2) + distance(Y1, Y2) + maxLenInTour

distance(X1, Y1) + 2 * minLenInSet > 
    distance(X1, X2) + distance(Y1, Y2) + maxLenInTour

distance(X1, Y1) + 2 * minLenInSet >
    distance(X1, X2) + 2 * maxLenInTour

3-opt with Neighbor Lists and DLB – simple

Compare the code below with the code for 2-opt with neighbor lists and DLB and notice that this algorithm is simplified, not truly bi-directional version.

proc One_City_3_Opt_ND(tour: var Tour_Array;
                       basePos: Tour_Index;
                       neighbor: Neighbor_Matrix;
                       position: var City_Position_In_Tour;
                       numberOfNeigbors: int;
                       DontLook: var DLB_Array): bool =
  type
    Move3 = object
      i, j, k: Tour_Index
      optCase: Reconnection_3_optCase
      gain: Length_Gain
  var
    improved: bool
    i, j, k: Tour_Index
    X1, X2, Y1, Y2, Z1, Z2: City_Number
    k_6, k_7:  Tour_Index
    Z1_6, Z2_6, Z1_7, Z2_7: City_Number
    gainExpected: Length_Gain
    goodMove: Move3
  
  improved = false
  goodMove.gain = 0
  block three_loops:
    
    for direction in [forward, backward]:
      if direction == forward:
        i  = basePos
      else:
        i  = (N + basePos - 1) mod N
      X1 = tour[i]
      X2 = tour[(i+1) mod N]
 
      for neighbor_1 in 0 .. numberOfNeigbors-1:
        # new edges in optCase=6: *X1-Y2*, Y1-Z1, X2-Z2
        # new edges in optCase=7: *X1-Y2*, X2-Z1, Y1-Z2
        Y2 = neighbor[X1][neighbor_1]
        j  = (position[Y2] + N - 1) mod N
        Y1 = tour[j]

        if (Y1 != X1) and (Y1 != X2):
          gainExpected = Gain_From_2_Opt(X1, X2, Y1, Y2)
          if gainExpected > goodMove.gain:
            improved = true
            goodMove = Move3(i: i,  j: j,  k: j,
                             optCase: opt3_case_1, gain: gainExpected)

        for neighbor_2 in 0 .. numberOfNeigbors-1:
          # new edges in optCase=6: X1-Y2, *Y1-Z1*, X2-Z2
          Z1_6 = neighbor[Y1][neighbor_2]
          k_6  = position[Z1_6]
          Z2_6 = tour[(k_6 + 1) mod N]
          # new edges in optCase=7: X1-Y2, X2-Z1, *Y1-Z2*
          Z2_7 = neighbor[Y1][neighbor_2]
          k_7  = (position[Z2_7] + N - 1) mod N
          Z1_7 = tour[k_7]
          
          if Between(i, j, k_6):
            gainExpected = Gain_From_3_Opt(X1, X2, Y1, Y2, Z1_6, Z2_6,
                                           opt3_case_6)
            if gainExpected > goodMove.gain:
              improved = true
              goodMove = Move3(i: i,  j: j,  k: k_6,
                               optCase: opt3_case_6, gain: gainExpected)

          if Between(i, j, k_7):
            gainExpected = Gain_From_3_Opt(X1, X2, Y1, Y2, Z1_7, Z2_7,
                                           opt3_case_7)
            if gainExpected > goodMove.gain:
              improved = true
              goodMove = Move3(i: i,  j: j,  k: k_7,
                               optCase: opt3_case_7, gain: gainExpected)

  if improved:
    Set_DLB_off(DontLook,
                  [ tour[goodMove.i], tour[(goodMove.i + 1) mod N],
                    tour[goodMove.j], tour[(goodMove.j + 1) mod N],
                    tour[goodMove.k], tour[(goodMove.k + 1) mod N] ])
    Make_3_Opt_Move(tour,
                    goodMove.i,  goodMove.j,  goodMove.k,
                    goodMove.optCase,  position)
  result = improved

Some simple statistics

Relative time, best length, average length and worst length for all tours created by NN heuristic and improved by some 2-opt and 3-opt algorithms. Random planar problems, N=100.

Problem #	Method	Time	Best	Avg	Worst
1	NN	100	100.00	111.12	122.44
	NN + 2opt	1360	88.94	92.58	98.91
	NN + 2opt Take Best	1140	89.47	91.75	95.40
	NN + 3opt	79166	87.99	89.61	92.82
	NN + 3opt-ND(24)	4586	87.99	89.86	94.22
2	NN	100	100.00	107.46	119.61
	NN + 2opt	1146	90.85	94.70	100.03
	NN + 2opt Take Best	1046	88.59	92.09	97.39
	NN + 3opt	69473	87.22	89.89	94.23
	NN + 3opt-ND(24)	4586	87.22	89.87	93.04
3	NN	100	100.00	106.81	117.56
	NN + 2opt	681	88.92	93.34	97.90
	NN + 2opt Take Best	881	89.38	92.28	96.22
	NN + 3opt	66112	87.85	89.53	93.19
	NN + 3opt-ND(24)	4106	87.68	89.07	92.68
4	NN	100	100.00	108.58	116.96
	NN + 2opt	1460	85.60	88.52	92.28
	NN + 2opt Take Best	1253	85.40	86.87	89.81
	NN + 3opt	73226	83.14	84.80	87.22
	NN + 3opt-ND(24)	4586	83.15	85.01	87.43
5	NN	100	100.00	110.59	124.43
	NN + 2opt	1146	92.73	95.42	98.30
	NN + 2opt Take Best	940	91.43	93.58	97.02
	NN + 3opt	72393	90.91	92.32	94.66
	NN + 3opt-ND(24)	4166	90.91	92.27	93.67
6	NN	100	100.00	107.15	118.12
	NN + 2opt	1253	89.35	92.71	98.25
	NN + 2opt Take Best	939	88.74	90.81	92.73
	NN + 3opt	69373	87.94	89.63	91.73
	NN + 3opt-ND(24)	4373	88.20	90.00	92.09
7	NN	100	100.00	110.20	121.48
	NN + 2opt	1253	92.30	95.31	98.94
	NN + 2opt Take Best	939	91.53	93.94	97.52
	NN + 3opt	65306	90.41	91.74	94.38
	NN + 3opt-ND(24)	4273	90.41	91.63	93.30
8	NN	100	100.00	105.35	112.99
	NN + 2opt	1253	88.20	92.02	97.83
	NN + 2opt Take Best	833	88.69	90.72	95.39
	NN + 3opt	66980	87.66	88.77	94.75
	NN + 3opt-ND(24)	4273	87.66	89.34	92.54
9	NN	100	100.00	107.15	115.57
	NN + 2opt	781	92.08	95.79	100.95
	NN + 2opt Take Best	875	90.94	94.06	98.73
	NN + 3opt	71681	89.76	91.35	94.37
	NN + 3opt-ND(24)	4393	89.76	91.32	93.54
10	NN	100	100.00	106.56	113.89
	NN + 2opt	1253	87.07	90.87	94.98
	NN + 2opt Take Best	1040	86.66	88.63	91.83
	NN + 3opt	71353	86.01	87.64	90.86
	NN + 3opt-ND(24)	4273	86.01	87.41	91.10

2-opt with Neighbor Lists and Best from given City – alternative

Note

This example of code does not contain any new ideas and can be skipped by a reader.

proc One_City_2_Opt_N_Best_From_City(tour: var Tour_Array;
                                     basePos: Tour_Index;
                                     neighbor: Neighbor_Matrix;
                                     position: var City_Position_In_Tour;
                                     numberOfNeigbors: int): bool =
  ## Makes the best of improving moves starting from given city
  type
    Move2 = object
      i, j: Tour_Index
      gain: Length_Gain
  var
    improved: bool
    i, i_pred, i_succ, j, j_pred, j_succ: Tour_Index
    X0, X0_pred, X0_succ, Y0, Y0_pred, Y0_succ: City_Number
    gainExpected: Length_Gain
    bestMove: Move2

  improved = false
  bestMove.gain = 0
  
  i   = basePos
  X0  = tour[i]
  i_pred = (N + i - 1) mod N
  i_succ = (i + 1) mod N
  X0_pred = tour[i_pred]
  X0_succ = tour[i_succ]

  for neighbor_number in 0 .. numberOfNeigbors-1:
    Y0 = neighbor[X0][neighbor_number]    
    j  = position[Y0]
    j_pred = (N + j - 1) mod N
    j_succ = (j + 1) mod N
    Y0_pred = tour[j_pred]
    Y0_succ = tour[j_succ]

    if (X0_succ != Y0) and (Y0_succ != X0):
      gainExpected  = Gain_From_2_Opt(X0, X0_succ, Y0, Y0_succ)
      if gainExpected > bestMove.gain:
        bestMove = Move2(i: i,  j: j,  gain: gainExpected)
        improved = true
    if (X0_pred != Y0) and (Y0_pred != X0):
      gainExpected = Gain_From_2_Opt(X0_pred, X0, Y0_pred, Y0)
      if gainExpected > bestMove.gain:
        bestMove = Move2(i: i_pred,  j: j_pred,  gain: gainExpected)
        improved = true
  if improved:
    Make_2_Opt_Move(tour, bestMove.i, bestMove.j, position)
  result = improved

2-opt with Neighbor Lists and Best from given City

proc One_City_2_Opt_N_Best_From_City(tour: var Tour_Array;
                                     basePos: Tour_Index;
                                     neighbor: Neighbor_Matrix;
                                     position: var City_Position_In_Tour;
                                     numberOfNeigbors: int): bool =
  ## Makes the best of improving moves starting from given city
  type
    Move2 = object
      i, j: Tour_Index
      gain: Length_Gain
  var
    improved: bool
    pos_Y2, i, j: Tour_Index
    X1, X2, Y1, Y2: City_Number
    gainExpected: Length_Gain
    bestMove: Move2

  improved = false
  bestMove.gain = 0
  block two_loops:
    for direction in [forward, backward]:
      if direction == forward:
        i  = basePos
        X1 = tour[i]
        X2 = tour[(i+1) mod N]
      else:
        i  = (N + basePos - 1) mod N
        X2 = tour[i]
        X1 = tour[(i+1) mod N]  # == tour[basePos]

      for neighbor_number in 0 .. numberOfNeigbors-1:
        Y1 = neighbor[X1][neighbor_number]
        if direction == forward:
          j = position[Y1]
          Y2 = tour[(j+1) mod N]
        else:
          j  = (N + position[Y1] - 1) mod N  # pos[Y1] == j+1
          Y2 = tour[j]

        if (X2 == Y1) or (Y2 == X1):
          continue

        gainExpected = Gain_From_2_Opt(X1, X2, Y1, Y2)
        if  gainExpected > bestMove.gain:
          bestMove = Move2(i: i,  j: j,  gain: gainExpected)
          improved = true
  if improved:
    Make_2_Opt_Move(tour, bestMove.i, bestMove.j, position)
  result = improved

proc LS_2_Opt_N_BC(tour: var Tour_Array;
                   neighbor: Neighbor_Matrix;
                   neighborListLen: int = N-1) =
  ## Optimizes tour using 2-opt with neighbor lists and fixed radius,
  ## for each city makes the best of improving moves that starts from it
  var
    locallyOptimal: bool = false
    improved: bool
    position: City_Position_In_Tour

  let numberOfNeigbors = min(neighborListLen, len(neighbor[0]))
  position = Create_Position_In_Tour(tour)

  while not locallyOptimal:
    locallyOptimal = true
    for basePos in 0 .. N-1:
      improved = One_City_2_Opt_N_Best_From_City(tour, basePos,
                                                 neighbor, position,
                                                 numberOfNeigbors)
      if improved:
        locallyOptimal = false

Some simple statistics

Relative time, best length, average length and worst length for all tours created by NN heuristic and improved by some 2-opt algorithms. Random planar problems, N=200.

Problem #	Method	Time	Best	Avg	Worst
1	NN	100	100.00	109.58	116.87
	NN + 2opt	4421	89.66	92.50	96.98
	NN + 2opt Take Best	4953	87.93	90.77	92.79
	NN + 2opt-NR(24)	265	88.94	91.32	95.01
	NN + 2opt-ND(24)	231	88.72	92.77	96.39
	NN + 2opt-NDR(24)	134	89.36	92.05	95.94
	NN + 2opt-N-BC(24)	497	89.20	92.18	96.75
	NN + 2opt-NR-BC(24)	234	88.42	90.90	95.06
2	NN	100	100.00	106.22	113.13
	NN + 2opt	4289	88.27	91.20	95.07
	NN + 2opt Take Best	5053	88.53	90.12	92.64
	NN + 2opt-NR(24)	265	87.78	90.26	93.83
	NN + 2opt-ND(24)	265	88.96	91.51	94.93
	NN + 2opt-NDR(24)	131	87.83	90.79	93.90
	NN + 2opt-N-BC(24)	499	87.83	90.74	93.97
	NN + 2opt-NR-BC(24)	231	87.11	89.91	93.20
3	NN	100	100.00	110.81	122.10
	NN + 2opt	3503	90.58	94.09	99.85
	NN + 2opt Take Best	4259	90.31	93.01	96.21
	NN + 2opt-NR(24)	201	90.76	93.92	97.50
	NN + 2opt-ND(24)	175	90.50	94.74	100.37
	NN + 2opt-NDR(24)	125	91.10	94.47	99.50
	NN + 2opt-N-BC(24)	403	90.69	94.09	99.06
	NN + 2opt-NR-BC(24)	177	90.26	93.54	97.36
4	NN	100	100.00	108.19	115.51
	NN + 2opt	4155	89.16	92.48	97.91
	NN + 2opt Take Best	4821	88.63	91.13	94.83
	NN + 2opt-NR(24)	331	88.18	91.57	96.27
	NN + 2opt-ND(24)	234	89.36	93.08	96.79
	NN + 2opt-NDR(24)	131	88.74	91.82	96.09
	NN + 2opt-N-BC(24)	597	88.44	92.33	96.40
	NN + 2opt-NR-BC(24)	199	88.23	91.66	95.12
5	NN	100	100.00	105.47	111.66
	NN + 2opt	3100	88.44	91.20	94.82
	NN + 2opt Take Best	3125	87.93	90.83	93.19
	NN + 2opt-NR(24)	198	87.55	90.33	93.73
	NN + 2opt-ND(24)	198	88.04	91.80	95.37
	NN + 2opt-NDR(24)	123	88.07	90.98	95.03
	NN + 2opt-N-BC(24)	322	88.06	91.21	94.08
	NN + 2opt-NR-BC(24)	173	87.45	90.11	94.13
6	NN	100	100.00	106.99	115.27
	NN + 2opt	5119	88.61	91.70	95.36
	NN + 2opt Take Best	5553	86.74	89.36	92.40
	NN + 2opt-NR(24)	297	87.67	90.39	93.91
	NN + 2opt-ND(24)	234	87.99	91.14	94.90
	NN + 2opt-NDR(24)	131	87.63	90.75	93.93
	NN + 2opt-N-BC(24)	500	87.49	90.55	93.80
	NN + 2opt-NR-BC(24)	231	86.69	89.70	91.95
7	NN	100	100.00	109.06	119.98
	NN + 2opt	4853	87.43	90.86	95.27
	NN + 2opt Take Best	4753	86.14	89.41	94.34
	NN + 2opt-NR(24)	265	85.74	89.66	95.14
	NN + 2opt-ND(24)	265	87.58	91.12	96.98
	NN + 2opt-NDR(24)	134	86.64	90.15	96.67
	NN + 2opt-N-BC(24)	565	86.83	90.74	95.90
	NN + 2opt-NR-BC(24)	331	85.83	89.68	95.53
8	NN	100	100.00	106.79	113.15
	NN + 2opt	4953	87.26	91.30	95.33
	NN + 2opt Take Best	5651	87.65	89.53	92.11
	NN + 2opt-NR(24)	231	88.07	90.27	93.96
	NN + 2opt-ND(24)	234	88.31	91.60	95.67
	NN + 2opt-NDR(24)	131	88.13	90.99	95.28
	NN + 2opt-N-BC(24)	500	87.81	90.11	94.61
	NN + 2opt-NR-BC(24)	231	87.36	89.20	92.07
9	NN	100	100.00	104.58	109.98
	NN + 2opt	4089	87.86	91.80	95.52
	NN + 2opt Take Best	5019	87.18	89.66	91.63
	NN + 2opt-NR(24)	399	87.96	90.73	94.13
	NN + 2opt-ND(24)	231	88.53	91.94	94.86
	NN + 2opt-NDR(24)	365	88.24	91.06	94.24
	NN + 2opt-N-BC(24)	929	88.03	91.00	94.26
	NN + 2opt-NR-BC(24)	531	87.73	90.29	93.77
10	NN	100	100.00	106.10	113.14
	NN + 2opt	6682	88.46	90.90	94.87
	NN + 2opt Take Best	5585	87.82	89.72	91.80
	NN + 2opt-NR(24)	297	87.94	90.73	95.79
	NN + 2opt-ND(24)	234	87.70	91.51	94.74
	NN + 2opt-NDR(24)	165	87.96	91.18	96.28
	NN + 2opt-N-BC(24)	563	87.88	90.56	93.03
	NN + 2opt-NR-BC(24)	234	86.44	89.79	92.67

2-opt with Neighbor Lists and optional DLB and FR

Note

This example of code does not contain any new ideas and can be skipped by a reader.

proc One_City_2_Opt_NDR(tour: var Tour_Array;
                        basePos: Tour_Index;
                        neighbor: Neighbor_Matrix;
                        position: var City_Position_In_Tour;
                        numberOfNeigbors: int;
                        DontLook: var DLB_Array;
                        useDLB: bool;
                        useFixedRadius: bool): bool =
  var
    improved: bool
    pos_Y2, i, j: Tour_Index
    X1, X2, Y1, Y2: City_Number
    d_X1_X2, d_X1_Y1: Length

  improved = false
  block two_loops:
    for direction in [forward, backward]:
      if direction == forward:
        i  = basePos
        X1 = tour[i]
        X2 = tour[(i+1) mod N]
      else:
        i  = (N + basePos - 1) mod N
        X2 = tour[i]
        X1 = tour[(i+1) mod N]  # == tour[basePos]
      d_X1_X2 = distance(X1, X2)

      for neighbor_number in 0 .. numberOfNeigbors-1:
        # after 2-opt move new tour would contain direct link X1-Y1,
        # so we look for city Y1 among cities that are close to city X1
        Y1 = neighbor[X1][neighbor_number]
        if direction == forward:
          j = position[Y1]
          Y2 = tour[(j+1) mod N]
        else:
          j  = (N + position[Y1] - 1) mod N  # pos[Y1] == j+1
          Y2 = tour[j]

        if (X2 == Y1) or (Y2 == X1):
          # 2-opt is for non-adjacent links only
          # by constuction X1 != Y1
          continue
          
        d_X1_Y1 = distance(X1, Y1)
          
        if useFixedRadius:
          if d_X1_Y1 > d_X1_X2:
            break

        # the same as: Gain_From_2_Opt(X1, X2, Y1, Y2) > 0
        if (d_X1_X2 + distance(Y1, Y2)) -
           (d_X1_Y1 + distance(X2, Y2)) > 0:
          if useDLB:
            Set_DLB_off(DontLook, [X1, X2, Y1, Y2])
          Make_2_Opt_Move(tour, i, j, position)
          improved = true
          break two_loops
  result = improved

proc LS_2_Opt_NDR(tour: var Tour_Array;
                  neighbor: Neighbor_Matrix;
                  neighborListLen: int = N-1;
                  useDLB: bool = true;
                  useFixedRadius: bool = true) =
  ## Optimizes the given tour using 2-opt with neighbor lists
  ## and with optional Don't Look Bits (DLB) and fixed radis
  var
    locallyOptimal: bool = false
    improved: bool
    DontLook: DLB_Array
    baseCity: City_Number
    position: City_Position_In_Tour

  let numberOfNeigbors = min(neighborListLen, len(neighbor[0]))
  position = Create_Position_In_Tour(tour)

  if useDLB:
    Set_DLB_off(DontLook, tour)  
  while not locallyOptimal:
    locallyOptimal = true
    for basePos in 0 .. N-1:
      baseCity = tour[basePos]
      if useDLB and isDLB_on(DontLook, baseCity):
        continue
      improved = One_City_2_Opt_NDR(tour, basePos, neighbor, position,
                                    numberOfNeigbors, DontLook,
                                    useDLB, useFixedRadius)
      if improved:
        locallyOptimal = false
      else:
        if useDLB:
          Set_DLB_on(DontLook, baseCity)