Created
November 11, 2023 23:51
-
-
Save dmishin/1a0f42bbdb606e63ac30516bef75e9e6 to your computer and use it in GitHub Desktop.
Use linear programming to find counterexample for a conjecture regarding the Traveling Salesman Problem
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
from scipy.optimize import linprog | |
import numpy as np | |
import itertools | |
#Searching for the counterexample for the following conjecture: | |
# in the shortest hamiltonian path, there is at least one | |
# node, connected to its 12 nearest neighbors | |
# https://www.reddit.com/r/math/comments/17rxc28/travelling_salesman_algorithm/ | |
# Use linear programming to find a counterexample | |
#FIrst define the number of nodes: | |
N = 5 | |
#Now the variables: they are pairwise distances between nodes. | |
#THere are N*(N-1)/2 of them | |
distance_vars = {} | |
idx = 0 | |
for i in range(N): | |
for j in range(i+1, N): | |
distance_vars[(i,j)] = idx | |
distance_vars[(j,i)] = idx | |
idx += 1 | |
M = idx | |
print(f"Total number of distance variables: {M}") | |
#Define objective as a sum of all distances, to be sure that the problem is always bounded | |
objective = np.array([1.0 for _ in range(M)]) | |
#Constraints are of several sets | |
#First constraints that tell us that the path 0-1-2-...-(N-1)-0 is hamiltonian. | |
#that is, d_01 + d_12 + d_23 < all other sums | |
A = [] | |
B = [] | |
#Enumerate all transpositions of the nodes 2..(N-1) | |
#Construct the duration of the initial cycle | |
initial_cycle_coeffs = np.zeros(M) | |
for i in range(N): | |
j = (i+1)%N | |
initial_cycle_coeffs[distance_vars[i,j]] = 1 | |
print("Constructed the function for the shortest cycle:", initial_cycle_coeffs) | |
#we would start enumerating all cycles from 0 | |
for idx, perm in enumerate(itertools.permutations(list(range(1, N)))): | |
#Skip the first element: it is the zero permutation, base value | |
if idx == 0: continue | |
#Every cycle would appear twice: in two different directions. | |
#Exclude one of them. | |
if perm[-1]<perm[0]: continue | |
cycle = (0,) + perm | |
cycle_coeffs = np.zeros(M) | |
for i in range(N): | |
ci = cycle[i] | |
cj = cycle[(i+1)%N] | |
cycle_coeffs[distance_vars[ci,cj]] = 1 | |
#Condition is that this cycle length is strictly bigger than the shortest cycle | |
# base_cycle_len - cycle_len <= -1 | |
A.append(initial_cycle_coeffs - cycle_coeffs) | |
B.append(-1) | |
#OK, we now have conditions telling that our cycle is indeed hamiltonian | |
#Now construct the counterexample conditions. | |
#For each node in the shortest cycle 0-1-2-...-(N-1)-0, | |
#(at least) one of 2 distances in the cycle | |
#must be bigger than one of the other N-3 distances | |
#This gives 2*(N-3) choices for one node, and (2*(N-3))^N choices total | |
# | |
# to make solution nondegenerate, we require that the difference is bigger than 1.0 | |
#and add the counterexample constrints to them | |
# construct choices for each node. | |
def node_counterexample_choices(n): | |
nprev = (n-1)%N | |
nnext = (n+1)%N | |
#Variables for 2 lengths in the cycle | |
len_forward = distance_vars[n, nnext] | |
len_backward = distance_vars[n, nprev] | |
#Variables for other distances, not used in the cycle | |
other_lens = [distance_vars[n, (n+i)%N] for i in range(2,N-1)] | |
assert len(other_lens) == N-3 | |
print(f"For node {n}, counterexample vars are: in_cycle={[len_forward,len_backward]}, out cycle={other_lens}") | |
return list(itertools.product([len_forward, len_backward], other_lens)) | |
all_counterexample_choices = [node_counterexample_choices(n) for n in range(N)] | |
#Define bounds for all variables | |
bounds = [(1.0, np.inf) for _ in range(M)] | |
cnt = 0 | |
for choice in itertools.product(*all_counterexample_choices): | |
#Choice is a list of pairs of 2 variables (variable for edge len in a loop) (variable for some other edge len) | |
#Take the hamiltonianity constraints... | |
A1 = A[:] | |
B1 = B[:] | |
print("##################################") | |
print("Counterexample choice:", choice) | |
#And add counterexample constraints to them | |
for var_in_cycle, var_out_cycle in choice: | |
assert var_in_cycle != var_out_cycle | |
#Sonstraint must say that the edge length in cycle is longer than edge outside. | |
# var_in_cycle >= var_out_cycle + 1.0 | |
# var_out_cycle - var_in_cycle <= -1.0 | |
a = np.zeros(M) | |
a[var_in_cycle] = -1 | |
a[var_out_cycle] = 1 | |
A1.append(a) | |
B1.append(-1.0) | |
print(f" constraint {a}.x <= -1") | |
#Now trying to solve the problem! | |
opt = linprog(c=objective, | |
A_ub = A1, | |
b_ub = B1, | |
bounds = bounds) | |
if opt.success: | |
print("COunterexample found!") | |
for i in range(N): | |
for j in range(i+1,N): | |
v = distance_vars[i,j] | |
print(f" d({i},{j})={opt.x[v]}") | |
break | |
else: | |
print("No counterexample") | |
cnt += 1 | |
else: | |
print(f"Tried {cnt} constructions, no counterexample found!") | |
exit(0) | |
D = np.zeros((N,N)) | |
for i in range(N): | |
for j in range(N): | |
if j == i: continue | |
v = distance_vars[i,j] | |
D[i,j] = opt.x[v] | |
D[j,i] = opt.x[v] | |
print("Distance matrix:") | |
print(repr(D)) | |
print("Checking cycle lengths") | |
for perm in itertools.permutations(list(range(1,N))): | |
cycle = (0,) + perm | |
cycle_len = sum(D[cycle[i], cycle[(i+1)%N]] for i in range(N)) | |
print(cycle,"len=",cycle_len) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment