Mean return time of king on chess board

Question

I am trying to find the stationary distribution of the markov chain of a king on a chess board, where a king can uniformly choose between any of its legal moves.

I have potential approaches, first, finding the transition matrix P for the entire board and then solving for π*P = π, where pi is the stationary distribution. I wrote a Python program that does this

import numpy as np
def is_on_board(i, j):
    return i >= 0 and j >= 0 and i < 8 and j < 8
P = [[0 for _ in range(64)] for _ in range(64)]
print(P)
for i in range(8):
    for j in range(8):
        curr = i8+j
        deg = 0
        for x in [-1, 0, 1]:
            for y in [-1, 0, 1]:
                if x != 0 or y != 0:
                    next = (i+x)  8 + (j+y)
                    deg += is_on_board(i+x, j+y)
        print((i,j), deg)
        for x in [-1, 0, 1]:
            for y in [-1, 0, 1]:
                if (x != 0 or y != 0) and is_on_board(i+x, j+y):
                    next = (i+x) * 8 + (j+y)
                    print(curr, next, 1/deg)
                    P[curr][next] = 1/deg
print(P)
import numpy as np
def calculate_eigenvectors(matrix):
    eigenvalues, eigenvectors = np.linalg.eig(matrix)
    return eigenvalues, eigenvectors
matrix = P
eigenvalues, eigenvectors = calculate_eigenvectors(np.array(matrix))
for i in range(len(eigenvalues)):
    eigenvalue = eigenvalues[i]
    eigenvector = eigenvectors[:, i]
    print("Eigenvalue:", eigenvalue)
    print("Eigenvector:", eigenvector)
    print()
    print("Check")
    print(np.dot(np.array(matrix), eigenvector))
    print(eigenvalue * eigenvector)
    break

Running this, I see that the eigenvector I get for the first eigenvalue (which is 1), is [1/8 1/8 .... 1/8]. Since the values should actually sum to 1, I can just divide this by 8 and get pi = [1/64 1/64 1/64 ... 1/64].

However, another approach is to find the mean hitting time of each position on the board, and then use the Fundamental Theorem of Markov Chains to find π_i = 1/T_ii where T_ii is the expected return time of the position i on the board.

I came across this paper here http://www.stevenheilman.org/~heilman/teach/505bexam1soln.pdf which finds a different mean hitting time for each point on the board, and thus would contradict the earlier solution I got.

Would appreciate some help with this problem.

Thanks!

Answer 1

The problem

As you write in your post the stationary distribution $\pi$ is the (normalized) left eigenvector of the transition matrix $P$. But in your code you compute the right eigenvector of $P$. Note that the uniform distribution is always a right eigenvector to eigenvalue $1$ of $P$, because the rows of $P$ must sum to $1$.

You can fix this by simply transposing the transition matrix $P$ before computing its eigenvalues (the resulting eigenvector will also need to be normalized). I tested your code with the fix and it works.

I post my own version of the code that i used to debug your code below. Maybe you or someone else is interested.

My version of the code

import numpy as np 
np.set_printoptions(precision=3)

We need to a function that checks if a pair of indices $(i,j)$ is still on the board:

def is_on_board(i,j):
    return 0<=i<8 and 0<=j<8

The main code to compute the transition matrix $P$ and the degree matrix is the following:

deg_matrix  = np.zeros(shape=(8,8),dtype=np.int16)
transition_matrix = np.zeros(shape=(8,8,8,8))
#compute the degree matrix and the transition matrix
for i in range(8):
    for j in range(8):
        for s in (-1,+1):
            if is_on_board(i+s,j): 
                transition_matrix[i,j,i+s,j]= 1 #add legal moves up down
            if is_on_board(i,j+s):
                transition_matrix[i,j,i,j+s] =1  #add legal moves left right
            if is_on_board(i+s,j+s):
                transition_matrix[i,j,i+s,j+s]=1#legal moves on one diagonal
            if is_on_board(i+s,j-s):
                transition_matrix[i,j,i+s,j-s]=1 #legal moves on other diagonal
            normf = np.sum(transition_matrix[i,j])
            deg_matrix[i,j]= normf
        transition_matrix[i,j] = transition_matrix[i,j]/normf
now transition_matrix[i,j,k,l] is the probability of the king to move from i,j to k,l

We compare the obtained degree matrix with the one presented in the linked document.

print(deg_matrix)

The output is:

[[3 5 5 5 5 5 5 3]
 [5 8 8 8 8 8 8 5]
 [5 8 8 8 8 8 8 5]
 [5 8 8 8 8 8 8 5]
 [5 8 8 8 8 8 8 5]
 [5 8 8 8 8 8 8 5]
 [5 8 8 8 8 8 8 5]
 [3 5 5 5 5 5 5 3]]

Which is exactly the same as in the document.

To compute the eigenvectors of the transition matrix we will need a square matrix. To obtain one we can flatten the transition matrix.

t_matrix_flat = transition_matrix.reshape((8*8,8*8))
# now t_matrix_flat[i+8j, k+8l] is the probability of the king to move from i,j to k,l

The stationary distribution $\pi$ is the left eigenvector of the flattened matrix to eigenvalue $1$.

eig_vals, eig_vecs = np.linalg.eig(t_matrix_flat.T) 
#left eigenvectors are right eigenvectors of the transpose
pi = eig_vecs[:,0].real #eigenvector of the transition matrix with biggest eigenvalue
pi= pi/np.sum(pi) #we have to normalize it

Printing the result in chessboard form:

print(pi.reshape((8,8)))# the stationary distribution calculated from the transition matrix

The output is:

[[0.007 0.012 0.012 0.012 0.012 0.012 0.012 0.007]
 [0.012 0.019 0.019 0.019 0.019 0.019 0.019 0.012]
 [0.012 0.019 0.019 0.019 0.019 0.019 0.019 0.012]
 [0.012 0.019 0.019 0.019 0.019 0.019 0.019 0.012]
 [0.012 0.019 0.019 0.019 0.019 0.019 0.019 0.012]
 [0.012 0.019 0.019 0.019 0.019 0.019 0.019 0.012]
 [0.012 0.019 0.019 0.019 0.019 0.019 0.019 0.012]
 [0.007 0.012 0.012 0.012 0.012 0.012 0.012 0.007]]

Now we can compare $\pi$ to the solution in the linked document from your post.

print(deg_matrix/420)# the stationary distribution according to the linked document

The output is:

[[0.007 0.012 0.012 0.012 0.012 0.012 0.012 0.007]
 [0.012 0.019 0.019 0.019 0.019 0.019 0.019 0.012]
 [0.012 0.019 0.019 0.019 0.019 0.019 0.019 0.012]
 [0.012 0.019 0.019 0.019 0.019 0.019 0.019 0.012]
 [0.012 0.019 0.019 0.019 0.019 0.019 0.019 0.012]
 [0.012 0.019 0.019 0.019 0.019 0.019 0.019 0.012]
 [0.012 0.019 0.019 0.019 0.019 0.019 0.019 0.012]
 [0.007 0.012 0.012 0.012 0.012 0.012 0.012 0.007]]