Objective

An abstract for our project is available on our home page. On this page we are going to discuss the R simulation that computes the mean travel time of \(m+1\) identically distributed Hamiltonian circuits sampled uniformly at random. \(m\) and \(n\) will be parameters to our simulation where \(n\) is the number of trials until equilibrium from which the \(n^{th}\) trial and the next \(m\) trials will be used to approximate the true expectation travel time, \[ \sum_{i} f(i)\mu(i) \approx \frac{1}{m+1} \sum_{k=0}^{m} f(X_k)\]
Towards this goal we first label 20 cities numerically from 1:20

U <-c(1:20) #vertices

In graph terminology these are our vertices and will make up the elements within our state space. Next we resolve the issue of forbidden travel and the missing component of our Hamiltonian graph, edges.

Distance Matrix

D <- matrix(runif(n=20*20, min=0, max =1), ncol = 20) #non-stochastic uniformly distributed distance matrix

We then manually specify 5 pairs \((i,j)\), \(i,j \in U\) for which travel is forbidden(i.e. the edge between \(i\), \(j\) is weighted 0). In the context of our simulation that means a \(0\) entry in the distance matrix.

D[4,6] = 0
D[9,13]= 0
D[1,19]= 0 
D[10,17]= 0  
D[3,8] = 0 #NO TRAVEL NO DISTANCE

We are almost done setting up our graph except for a couple outstanding requirements.First the edge connecting a pair of cities is more than likely not identical. Second, travel from city \(i\) to city \(i\) for \(i \in U\) is greater than \(0\) and most importantly we still haven’t assigned the other half of the travel restrictions. To satisfy these requirements, we transform \(D\) into a symmetric matrix by assigning the lower triangular of \(D\) the entries of the upper triangular of the transpose of D and explicitly set the diagonal entries to zero. Thus as you can see there are 185 pairs of vertices with positive real distance, \((19*20)/2 - 5\) from observing \(D\) or \({20 \choose 2}-5\).

D[lower.tri(D)] = t(D)[lower.tri(D)] #symmetric matrix
diag(D) <- rep(0,20) #diagonals are 0

Notice that the forbidden pairs were entered in the form of D[i,j] where i<j or else when we go to make the matrix symmetric the upper triangular entry is still non-zero. This can be remedied by manually specifying all 10 entries \(0\), or specifying D[i,j], j<i and assigning the upper-triangular of D with the upper-triangular of the transpose of D like this

D[upper.tri(D)] = t(D)[upper.tri(D)] #symmetric matrix

or in a more proper application handling user input.


Algorithm

Thus with the vertices set up and the weighted-edges specified, we have initialized our graph. We can now consider our \(Q\). The most prudent idea as suggested is a method which picks two points uniformly at random from a valid sequence and interchanges their position. The code chunk below is portion of the function that accomplishes this task and our broader objective. For a more comprehensive glance at the code refer to the Appendix.

tspMCMC <- function(n,m){
v <- sample(U,size = 20, replace=FALSE) #without repetition sample first sequence
while(indicatorPermissible(v)==0){#make sure its valid
  v <- sample(U,size=20, replace=FALSE) 
}
S <- list(v) #record it in the state space
...
...
...
while(i <= m+n){
  P <- sample(U,size=1) #basic sampling 
  Q <- sample(U,size=1) #likewise for proposed state, 
  index_p = match(P, unlist(S[i])) #the position of P in the sequence
  index_q = match(Q, unlist(S[i])) #likewise for the second element in the pair
  k <- unlist(S[i])  
  z <- c(replace(k, c(index_p,index_q), c(Q,P))) #interchange P and Q 
  if(indicatorPermissible(z) == 1){#if the new sequence is valid
    if(i >= n){#has equilibrium been reached yet?
      travel_time[j] = distance_travelled(S[i]) #store travel time
      j <- j + 1
    }
    S[i+1] <- list(z) #adds the new sequence to collection of sampled state spaces
    i <- i + 1 
  }
  else{
    if(i<=n){#this is done to make sure we dont double count travel times
      S[i+1] <- S[i]
      i <- i + 1
    }
    rejection_counter = rejection_counter+1 #see how many time our proposed transition introduced forbidden travel
  }
}
}

There is alot to unpack but very concisely we take two arguments \(n\) (integer) and \(m\) (integer). We then randomly sample an initial sequence \(X_1\) without-repetition which we check is valid and add to our record of sampled sequences. From there we sample two points, extracts their place in the sequence and use the last recorded valid sequence, \(X_{i}\) to propose a new sequence. Since we are using a Metropolis with \(\mu = 1\) (i.e. no rejection), we are only responsible for checking that the proposed sequence \(X_{i+1}\) is also valid. This is where indicatorPermissable (\(X_{i+1}\)) comes in. It’s the responsibility of the indicatorPermissable function to return \(0\) if the proposed interchange will create a sequence where travel between two forbidden cities is introduced. In the case where the indicator returns \(1\), we accept the new sequence as valid and add it to the record (i.e list) of valid circuits and increment \(i\).[1] Furthermore if we have already reached equilibrium or \(n\) trials, we also start computing the travel times by passing the sequence \(X_{i}\) to a function called distance_travelled(\(X_{i}\)). This function takes the sequence \(X_i\) not \(X_{i+1}\) so we can compute the travel time for \(X_n\) once equilibrium is reached and tallies the distance \(i \to i+1\) by starting at i and consulting D[i,i+1] \(\forall i \in (1:19)\). The special case of the endpoint transitions \(20 \to 1\) is manually treated.

[1]As you can see in the code chunk above, our control sequence of choice is a while loop. The reason being, unlike a for-loop we have more control on our index \(i\).For example, until the simulation reaches equilibrium we aren’t entirely concerned with repeat sequences. When we go to compute travel-times(\(i \geq n\)) however we control for repeat samples by ignoring trials for which indicatorPermissable returns 0 and trying again at index \(i\). The motivation being that our approximation especially at lower bounds will be skewed by many repeat sequences.


Observations

Now let us test the algorithm using the same arbitrary forbidden cities from above and a \(Uni(0,1)\) distributed distance matrix, additionally choose \(n=150\) and \(m=250\) as a starting point for our bounds.

tspMCMC(150,250)
## Average Return Time: 10.6136 units

Furthermore re-initializing our graph but keeping the distance matrix fixed for multiple arbitrary choices of n and m.

tspMCMC(150,250)
tspMCMC(250,500)
tspMCMC(500,1000)
tspMCMC(1000,2000)
## tspMCM( 150 , 250 ) 
## Rejection Counter:  41 
## Average Return Time: 9.986814 units
## tspMCM( 250 , 500 ) 
## Rejection Counter:  88 
## Average Return Time: 10.09586 units
## tspMCM( 500 , 1000 ) 
## Rejection Counter:  168 
## Average Return Time: 10.08226 units
## tspMCM( 1000 , 2000 ) 
## Rejection Counter:  314 
## Average Return Time: 10.09558 units
## Computational Length: 0.3000782 seconds


As you can see in these observations, our mean travel times stays within a band of one unit even as we double or quadruple our bounds \(n\) and \(m\) and while \(D\) is fixed. This variation is inherently linked with the type of positive real density we use to sample \(D\). We note that we would observe wilder swings if we used \(D \sim Uni(0,200)\) as an example. Furthermore we conclude by settling on bounds \(n = 1000\) and \(m=2000\) for practical reasons. Even though it is computationally feasible to pass parameters with \(n\) and \(m\) in the millions, there are \(20! = ~2.432902e+18\) sequences rendering such an idea frivolous by the sheer magnitude of our state space \(S\).

Appendix

Below is the code.

 

By Howell Lu, Adeeb Rouhani , Sebhat Yacob