Task 1

Subtask a.

In [1]:
import numpy as np
from numpy import linalg as la
np.set_printoptions(precision=3,suppress=True)
In [2]:
A1=np.array([[1,0,0,1/2,0,0],[0,0,1,1/2,1/3,1/2],[0,1,0,0,0,0],[0,0,0,0,1/3,0],[0,0,0,0,0,1/2],[0,0,0,0,1/3,0]])
print(A1)
la.det(A1)
[[1.    0.    0.    0.5   0.    0.   ]
 [0.    0.    1.    0.5   0.333 0.5  ]
 [0.    1.    0.    0.    0.    0.   ]
 [0.    0.    0.    0.    0.333 0.   ]
 [0.    0.    0.    0.    0.    0.5  ]
 [0.    0.    0.    0.    0.333 0.   ]]
Out[2]:
0.0

$A$ is a stochastic matrix because all its columns sum up to 1. We apply the Power method to see where the Markov chain converges.

  • if it starts at page 6 and takes an even number of iterations
In [3]:
x_0=np.array([0,0,0,0,0,1])
x_n=x_0
for i in range(22):
    x_n = np.dot(A1,x_n)
x_n
Out[3]:
array([0.1, 0.2, 0.7, 0. , 0. , 0. ])
  • if it starts at page 6 and takes an odd number of iterations
In [4]:
x_0=np.array([0,0,0,0,0,1])
x_n=x_0
for i in range(21):
    x_n = np.dot(A1,x_n)
x_n
Out[4]:
array([0.1, 0.7, 0.2, 0. , 0. , 0. ])
  • if it starts at page 4 and takes an even number of iterations
In [5]:
x_0=np.array([0,0,0,1,0,0])
x_n=x_0
for i in range(22):
    x_n = np.dot(A1, x_n)
x_n
Out[5]:
array([0.5, 0. , 0.5, 0. , 0. , 0. ])
  • if it starts at page 4 and takes an odd number of iterations
In [6]:
x_0=np.array([0,0,0,1,0,0])
x_n = x_0
for i in range(21):
    x_n = np.dot(A1,x_n)
x_n
Out[6]:
array([0.5, 0.5, 0. , 0. , 0. , 0. ])

So the Markov process does not converge.

Remember that: If a Markov chain converges to a steady-state vector $\vec x$, if $\lambda_1=1$ is a dominat eigenvalue of $A$.

Let's have a look at the eigenvalues of the transition matrix:

In [7]:
la.eig(A1)
Out[7]:
(array([ 1.   ,  1.   , -1.   ,  0.   ,  0.408, -0.408]),
 array([[ 1.   ,  0.   ,  0.   , -0.408, -0.308, -0.173],
        [ 0.   ,  0.707, -0.707,  0.   , -0.251,  0.141],
        [ 0.   ,  0.707,  0.707, -0.408, -0.615, -0.346],
        [ 0.   ,  0.   , -0.   ,  0.816,  0.364,  0.487],
        [ 0.   ,  0.   , -0.   ,  0.   ,  0.446, -0.597],
        [ 0.   ,  0.   , -0.   ,  0.   ,  0.364,  0.487]]))

There are three eigenvalues with absolute value equal to 1, ie, $|\lambda_i|=1, i=1,2,3$. In general, a stochastic matrix has a dominant eigenvalue equal to one but not in this case.

A stochastic matrix has the following properties. The largest absolute value of a stochastic matrix is at most 1 by Gershgorin circle theorem (not discussed in class). Additionally, every stochastic matrix has an "obvious" column eigenvector associated to the eigenvalue 1: the vector ${\boldsymbol {1}}$, whose coordinates are all equal to 1.

On the other hand, Perron theorem applied to stochastic matrices tells us that if the stochastic matrix is positive then it has a dominant eigenvalue $\lambda = 1$. More generally, Frobenius theorem tells us that if the stochastic matrix is nonnegative and irreducible then again it has a dominant eigenvalue $\lambda = 1$. However, in general stochastic matrices need not be positive or irreducible.

In the next subtask we try to modify the transition matrix such that it has all entries positive (no zeros) such that Perron's theorem applies and the corresponding process converges.

A Markov process with transition matrix $A$ is said to be regular if all the entries of some power of $A$ are positive. It can be shown that if this happens the $A$ has a dominant eigenvalue equal to 1 as well.

Subtask b.

In [8]:
A2 = 1/6*np.ones((6,6))
In [9]:
x_0=np.array([0,0,0,1,0,0])
x_n = x_0
for i in range(21):
    x_n = np.dot(A2,x_n)
x_n
Out[9]:
array([0.167, 0.167, 0.167, 0.167, 0.167, 0.167])

Hence, every node is eqaully likely and we gain no information.

In [10]:
ell, P = la.eig(A2)
ell, P
Out[10]:
(array([ 1., -0.,  0.,  0.,  0.,  0.]),
 array([[ 0.408,  0.   ,  0.767, -0.   , -0.   , -0.   ],
        [ 0.408,  0.894,  0.331, -0.134,  0.   ,  0.   ],
        [ 0.408, -0.224, -0.275,  0.89 ,  0.   ,  0.   ],
        [ 0.408, -0.224, -0.275, -0.252, -0.577, -0.577],
        [ 0.408, -0.224, -0.275, -0.252,  0.789, -0.211],
        [ 0.408, -0.224, -0.275, -0.252, -0.211,  0.789]]))

The matrix is positive hence it is a dominant eigenvalue equal to 1. The corresponding eigenvector is a vectors of entries or 1 or normalized as in the matrix P above.

If we can calculate the eigenvalues because computationally feasibile as in these small examples, we can also find the steady-state vectors by applying the theory seen on slide 20.

In [11]:
P_inv = la.inv(P) # this is also computationally demanding
z_0 = P_inv @ x_0
x_n = z_0[0] * P[:,0]
print(x_n)
[0.167 0.167 0.167 0.167 0.167 0.167]
In [12]:
### Subtask c.
In [13]:
A = 0.85 * A1 + 0.15 * A2
In [14]:
x_0=np.array([0,0,0,1,0,0])
x_n = x_0
for i in range(21):
    x_n = np.dot(A,x_n)
x_n
Out[14]:
array([0.277, 0.324, 0.285, 0.036, 0.041, 0.036])

Hence page 2 has the highest ranking

In [15]:
ell, P = la.eig(A)
ell, P
Out[15]:
(array([ 1.   ,  0.85 , -0.85 ,  0.347, -0.347, -0.   ]),
 array([[ 0.522,  0.816,  0.   ,  0.308, -0.173,  0.408],
        [ 0.618, -0.408,  0.707,  0.251,  0.141,  0.   ],
        [ 0.574, -0.408, -0.707,  0.615, -0.346,  0.408],
        [ 0.071, -0.   , -0.   , -0.364,  0.487, -0.816],
        [ 0.078, -0.   ,  0.   , -0.446, -0.597,  0.   ],
        [ 0.071, -0.   , -0.   , -0.364,  0.487,  0.   ]]))
In [16]:
P_inv = la.inv(P) # this is also computationally demanding
z_0 = P_inv @ x_0
x_n = z_0[0] * P[:,0]
print(x_n)
[0.27  0.32  0.297 0.036 0.041 0.036]

The slight discrepancies might be due to rounding erros or to an insufficient number of iterations, 21, in the power method.

Task 2

In [17]:
with open("../assets/top250movies.txt", encoding="utf-8") as f:
    lines = f.readlines()
In [18]:
db = {} 
for line in lines:
    entries = line.strip().split("/")
    db[entries[0]] = entries[1:]

To handle the issue of weights due to repeated co-presences of actors in movies, we can first create a multi directed graph and then convert it in a directed graph with weights on arcs. A multi digraph is a graph that allows multiple arcs between nodes. Alternatively, as hinted by the text of the exercise we could create a adjacency dictionary where for every actor we list the actors that are reached by the first actor (ie, more expensive actors) allowing repeated entries. Then we can construct the digraph using the adjacency dictionary. However, graphs in networkx automatically construct adjacency lists and since library functions are to be preferred in Python because more efficient, we prefer to use the first alternative with multi digraph then converted in digraph.

In [19]:
import networkx as nx

MDG = nx.MultiDiGraph()
for k in db:
    for i in range(len(db[k])):
        actor = db[k][i]
        MDG.add_edges_from([(cheaper_actor,actor) for cheaper_actor in db[k][(i+1):]])
In [20]:
MDG.number_of_nodes(), MDG.number_of_edges()
Out[20]:
(14882, 886259)
In [21]:
DG = nx.DiGraph()
for node, outgoing_neighbors in MDG.adjacency():
    for neighbor, arc_dict in outgoing_neighbors.items():
        value = len(arc_dict.values())
        DG.add_edge(node, neighbor, weight = value)
In [22]:
DG.number_of_nodes(), DG.number_of_edges()
Out[22]:
(14882, 880639)
In [23]:
PR = nx.pagerank(DG, alpha=0.7)
In [24]:
sorted(PR, key=PR.get, reverse=True)[0:5]
Out[24]:
['Leonardo DiCaprio', 'Robert De Niro', 'Tom Hanks', 'Jamie Foxx', 'Al Pacino']