page_rank.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "### Introduction\n",
    "\n",
    "PageRank (developed by Larry Page and Sergey Brin) revolutionized web search by generating a\n",
    "ranked list of web pages based on the underlying connectivity of the web. The PageRank algorithm is\n",
    "based on an ideal random web surfer who, when reaching a page, goes to the next page by clicking on a\n",
    "link. The surfer has equal probability of clicking any link on the page and, when reaching a page with no\n",
    "links, has equal probability of moving to any other page by typing in its URL. In addition, the surfer may\n",
    "occasionally choose to type in a random URL instead of following the links on a page. The PageRank is\n",
    "the ranked order of the pages from the most to the least probable page the surfer will be viewing.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "# Let's load the libraries\n",
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "np.set_printoptions(suppress=True)\n",
    "\n",
    "def generate_internet(n) :\n",
    "    c = np.full([n,n], np.arange(n))\n",
    "    c = (abs(np.random.standard_cauchy([n,n])/2) > (np.abs(c - c.T) + 1)) + 0\n",
    "    c = (c+1e-10) / np.sum((c+1e-10), axis=0)\n",
    "    return c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "### PageRank as a linear algebra problem\n",
    "Let's imagine a micro-internet, with just 6 websites (**A**vocado, **B**ullseye, **C**atBabel, **D**romeda, **e**Tings, and **F**aceSpace).\n",
    "Each website links to some of the others.\n",
    "\n",
    "The design principle of PageRank is that important websites will be linked to by important websites.\n",
    "This somewhat recursive principle will form the basis of our thinking.\n",
    "\n",
    "Imagine we have 100 *Procrastinating Pat*s on our micro-internet, each viewing a single website at a time.\n",
    "Each minute the Pats follow a link on their website to another site on the micro-internet.\n",
    "After a while, the websites that are most linked to will have more Pats visiting them, and in the long run, each minute for every Pat that leaves a website, another will enter keeping the total numbers of Pats on each website constant.\n",
    "The PageRank is simply the ranking of websites by how many Pats they have on them at the end of this process.\n",
    "\n",
    "We represent the number of Pats on each website with the vector,\n",
    "$$\\mathbf{r} = \\begin{bmatrix} r_A \\\\ r_B \\\\ r_C \\\\ r_D \\\\ r_E \\\\ r_F \\end{bmatrix}$$\n",
    "And say that the number of Pats on each website in minute $i+1$ is related to those at minute $i$ by the matrix transformation\n",
    "\n",
    "$$ \\mathbf{r}^{(i+1)} = L \\,\\mathbf{r}^{(i)}$$\n",
    "with the matrix $L$ taking the form,\n",
    "$$ L = \\begin{bmatrix}\n",
    "L_{A→A} & L_{B→A} & L_{C→A} & L_{D→A} & L_{E→A} & L_{F→A} \\\\\n",
    "L_{A→B} & L_{B→B} & L_{C→B} & L_{D→B} & L_{E→B} & L_{F→B} \\\\\n",
    "L_{A→C} & L_{B→C} & L_{C→C} & L_{D→C} & L_{E→C} & L_{F→C} \\\\\n",
    "L_{A→D} & L_{B→D} & L_{C→D} & L_{D→D} & L_{E→D} & L_{F→D} \\\\\n",
    "L_{A→E} & L_{B→E} & L_{C→E} & L_{D→E} & L_{E→E} & L_{F→E} \\\\\n",
    "L_{A→F} & L_{B→F} & L_{C→F} & L_{D→F} & L_{E→F} & L_{F→F} \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "where the columns represent the probability of leaving a website for any other website, and sum to one.\n",
    "The rows determine how likely you are to enter a website from any other, though these need not add to one.\n",
    "The long time behaviour of this system is when $ \\mathbf{r}^{(i+1)} = \\mathbf{r}^{(i)}$, so we'll drop the superscripts here, and that allows us to write,\n",
    "$$ L \\,\\mathbf{r} = \\mathbf{r}$$\n",
    "\n",
    "which is an eigenvalue equation for the matrix $L$, with eigenvalue 1 (this is guaranteed by the probabalistic structure of the matrix $L$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "# Define L\n",
    "L = np.array([[0,   1/2, 1/3, 0, 0,   0 ],\n",
    "              [1/3, 0,   0,   0, 1/2, 0 ],\n",
    "              [1/3, 1/2, 0,   1, 0,   1/2 ],\n",
    "              [1/3, 0,   1/3, 0, 1/2, 1/2 ],\n",
    "              [0,   0,   0,   0, 0,   0 ],\n",
    "              [0,   0,   1/3, 0, 0,   0 ]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "In principle, we could use a linear algebra library, as below, to calculate the eigenvalues and vectors.\n",
    "And this would work for a small system. But this gets unmanagable for large systems.\n",
    "And since we only care about the principal eigenvector (the one with the largest eigenvalue, which will be 1 in this case), we can use the *power iteration method* which will scale better, and is faster for large systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([16.        ,  5.33333333, 40.        , 25.33333333,  0.        ,\n",
       "       13.33333333])"
      ]
     },
     "execution_count": 5,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eVals, eVecs = la.eig(L) # Gets the eigenvalues and vectors\n",
    "order = np.absolute(eVals).argsort()[::-1] # Orders them by their eigenvalues\n",
    "eVals = eVals[order]\n",
    "eVecs = eVecs[:,order]\n",
    "\n",
    "r = eVecs[:, 0] # Sets r to be the principal eigenvector\n",
    "100 * np.real(r / np.sum(r)) # Make this eigenvector sum to one, then multiply by 100 Procrastinating Pats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "We can see from this list, the number of Procrastinating Pats that we expect to find on each website after long times.\n",
    "Putting them in order of *popularity* (based on this metric), the PageRank of this micro-internet is:\n",
    "\n",
    "**C**atBabel, **D**romeda, **A**vocado, **F**aceSpace, **B**ullseye, **e**Tings\n",
    "\n",
    "Let's now try to get the same result using the Power-Iteration method.\n",
    "This method will be much better at dealing with large systems.\n",
    "\n",
    "First let's set up our initial vector, $\\mathbf{r}^{(0)}$, so that we have our 100 Procrastinating Pats equally distributed on each of our 6 websites."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([16.66666667, 16.66666667, 16.66666667, 16.66666667, 16.66666667,\n",
       "       16.66666667])"
      ]
     },
     "execution_count": 6,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = 100 * np.ones(6) / 6 # Sets up this vector (6 entries of 1/6 × 100 each)\n",
    "r # Shows it's value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Next, let's update the vector to the next minute, with the matrix  L."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([13.88888889, 13.88888889, 38.88888889, 27.77777778,  0.        ,\n",
       "        5.55555556])"
      ]
     },
     "execution_count": 7,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = L @ r # Apply matrix L to r\n",
    "r # Show it's value\n",
    "# Re-run this cell multiple times to converge to the correct answer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "We can automate applying this matrix multiple times as follows,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([16.        ,  5.33333333, 40.        , 25.33333333,  0.        ,\n",
       "       13.33333333])"
      ]
     },
     "execution_count": 8,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = 100 * np.ones(6) / 6 # Sets up this vector (6 entries of 1/6 × 100 each)\n",
    "for i in np.arange(100) : # Repeat 100 times\n",
    "    r = L @ r\n",
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Or even better, we can keep running until we get to the required tolerance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "18 iterations to convergence.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([16.00149917,  5.33252025, 39.99916911, 25.3324738 ,  0.        ,\n",
       "       13.33433767])"
      ]
     },
     "execution_count": 9,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = 100 * np.ones(6) / 6 # Sets up this vector (6 entries of 1/6 × 100 each)\n",
    "lastR = r\n",
    "r = L @ r\n",
    "i = 0\n",
    "while la.norm(lastR - r) > 0.01 :\n",
    "    lastR = r\n",
    "    r = L @ r\n",
    "    i += 1\n",
    "print(str(i) + \" iterations to convergence.\")\n",
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "See how the PageRank order is established fairly quickly, and the vector converges on the value we calculated earlier after a few tens of repeats."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "### Damping Parameter\n",
    "The system we just studied converged fairly quickly to the correct answer.\n",
    "Let's consider an extension to our micro-internet where things start to go wrong.\n",
    "\n",
    "Intuitively, only *FaceSpace*, which is in the bottom half of the page rank, links to this website amongst the two others it links to,\n",
    "so we might expect *Geoff's* site to have a correspondingly low PageRank score.\n",
    "\n",
    "Build the new $L$ matrix for the expanded micro-internet, and use Power-Iteration on the Procrastinating Pat vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    " # We'll call this one L2, to distinguish it from the previous L.\n",
    "L2 = np.array([[0,   1/2, 1/3, 0, 0,   0, 0 ],\n",
    "               [1/3, 0,   0,   0, 1/2, 0, 0 ],\n",
    "               [1/3, 1/2, 0,   1, 0,   1/3, 0 ],\n",
    "               [1/3, 0,   1/3, 0, 1/2, 1/3, 0 ],\n",
    "               [0,   0,   0,   0, 0,   0, 0 ],\n",
    "               [0,   0,   1/3, 0, 0,   0, 0 ],\n",
    "               [0,   0,   0,   0, 0,   1/3, 1 ]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "131 iterations to convergence.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 0.03046998,  0.01064323,  0.07126612,  0.04423198,  0.        ,\n",
       "        0.02489342, 99.81849527])"
      ]
     },
     "execution_count": 11,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = 100 * np.ones(7) / 7 # Sets up this vector (6 entries of 1/6 × 100 each)\n",
    "lastR = r\n",
    "r = L2 @ r\n",
    "i = 0\n",
    "while la.norm(lastR - r) > 0.01 :\n",
    "    lastR = r\n",
    "    r = L2 @ r\n",
    "    i += 1\n",
    "print(str(i) + \" iterations to convergence.\")\n",
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "That's no good! *Geoff* seems to be taking all the traffic on the micro-internet, and somehow coming at the top of the PageRank.\n",
    "This behaviour can be understood, because once a Pat get's to *Geoff's* Website, they can't leave, as all links head back to Geoff.\n",
    "\n",
    "To combat this, we can add a small probability that the Procrastinating Pats don't follow any link on a webpage, but instead visit a website on the micro-internet at random.\n",
    "We'll say the probability of them following a link is $d$ and the probability of choosing a random website is therefore $1-d$.\n",
    "We can use a new matrix to work out where the Pat's visit each minute.\n",
    "$$ M = d \\, L + \\frac{1-d}{n} \\, J $$\n",
    "where $J$ is an $n\\times n$ matrix where every element is one.\n",
    "\n",
    "If $d$ is one, we have the case we had previously, whereas if $d$ is zero, we will always visit a random webpage and therefore all webpages will be equally likely and equally ranked.\n",
    "For this extension to work best, $1-d$ should be somewhat small - though we won't go into a discussion about exactly how small.\n",
    "\n",
    "Let's retry this PageRank with this extension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "d = 0.5\n",
    "M = d * L2 + (1-d)/7 * np.ones([7, 7])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8 iterations to convergence.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([13.68217054, 11.20902965, 22.41964343, 16.7593433 ,  7.14285714,\n",
       "       10.87976354, 17.90719239])"
      ]
     },
     "execution_count": 15,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = 100 * np.ones(7) / 7 # Sets up this vector (6 entries of 1/6 × 100 each)\n",
    "lastR = r\n",
    "r = M @ r\n",
    "i = 0\n",
    "while la.norm(lastR - r) > 0.01 :\n",
    "    lastR = r\n",
    "    r = M @ r\n",
    "    i += 1\n",
    "print(str(i) + \" iterations to convergence.\")\n",
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "This is certainly better, the PageRank gives sensible numbers for the Procrastinating Pats that end up on each webpage.\n",
    "This method still predicts Geoff has a high ranking webpage however.\n",
    "This could be seen as a consequence of using a small network. We could also get around the problem by not counting self-links when producing the L matrix (an if a website has no outgoing links, make it link to all websites equally).\n",
    "We won't look further down this route, as this is in the realm of improvements to PageRank, rather than eigenproblems.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Now we produce a function that can calculate the PageRank for an arbitrarily large probability matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "def pageRank(linkMatrix, d) :\n",
    "    n = linkMatrix.shape[0]\n",
    "\n",
    "    M = d * linkMatrix + (1-d) / n * np.ones([n, n])\n",
    "    r = 100 * np.ones(n) / n\n",
    "    last_r = r \n",
    "    r = M @ r\n",
    "    i = 0\n",
    "\n",
    "    while la.norm(last_r - r) > 0.01:\n",
    "        last_r = r\n",
    "        r = M @ r\n",
    "        i += 1\n",
    "\n",
    "    return r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Test it,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "generate_internet(5)\n",
    "\n",
    "L = generate_internet(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.69659305,  1.69659305,  1.69659305,  2.54571186, 30.50751672,\n",
       "       35.58567689, 13.55601964,  9.32210963,  1.69659305,  1.69659305])"
      ]
     },
     "execution_count": 21,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pageRank(L, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.69491526,  1.69491526,  1.69491526,  2.54237289, 30.50847456,\n",
       "       35.59322031, 13.55932203,  9.3220339 ,  1.69491526,  1.69491526])"
      ]
     },
     "execution_count": 22,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eVals, eVecs = la.eig(L) # Gets the eigenvalues and vectors\n",
    "order = np.absolute(eVals).argsort()[::-1] # Orders them by their eigenvalues\n",
    "eVals = eVals[order]\n",
    "eVecs = eVecs[:,order]\n",
    "\n",
    "r = eVecs[:, 0]\n",
    "100 * np.real(r / np.sum(r))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (system-wide)",
   "language": "python",
   "metadata": {
    "cocalc": {
     "description": "Python 3 programming language",
     "priority": 100,
     "url": "https://www.python.org/"
    }
   },
   "name": "python3",
   "resource_dir": "/ext/jupyter/kernels/python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}