Description
The purpose of this project1
is to learn the importance of numerical linear algebra algorithms to solve fundamental
linear algebra problems that occur in search engines.
Page-Rank Algorithm – Linear Systems
One of the reasons why GoogleTM is such an effective search engine is the Page-RankTM algorithm developed by
Google’s founders, Larry Page and Sergey Brin, when they were graduate students at Stanford University. PageRank
is determined entirely by the link structure of the World Wide Web. It is recomputed about once a month and does
not involve the actual content of any Web pages or individual queries. Then, for any particular query, Google finds
the pages on the Web that match that query and lists those pages in the order of their PageRank. Imagine surfing the
Web, going from page to page by randomly choosing an outgoing link from one page to get to the next. This can lead
to dead ends at pages with no outgoing links, or cycles around cliques of interconnected pages. So, a certain fraction
of the time, simply choose a random page from the Web. This theoretical random walk is known as a Markov chain
or Markov process. The limiting probability that an infinitely dedicated random surfer visits any particular page is its
PageRank. A page has high rank if other pages with high rank link to it.
Let W be the set of Web pages that can be reached by following a chain of hyperlinks starting at some root page
and let n be the number of pages in W. For Google, the set W highly fluctuates over time and is not exactly known,
but was estimated to be about 30 billion in 20162
. Let G be the n-by-n connectivity matrix of a portion of the Web,
that is gij = 1 if there is a hyperlink to page i from page j and zero otherwise. The matrix G can be huge, but it is
very sparse. Its j-th column shows the links on the j-th page. The number of nonzeros in G is the total number of
hyperlinks in W. Let ri and cj be the row and column sums of G:
ri =
X
j
gij and cj =
X
i
gij . (1)
The quantities rj and cj are the in-degree and out-degree of the j-th page. Let p be the probability that the random
walk follows a link. A typical value is p = 0.85. Then 1 − p is the probability that some arbitrary page is chosen and
δ = (1 − p)/n is the probability that a particular random page is chosen. Let A be a n-by-n matrix with elements
aij =
p gij/cj + δ cj 6= 0,
1/n cj = 0.
(2)
Notice that A comes from scaling the connectivity matrix by its column sums. The j-th column is the probability of
jumping from the j-th page to the other pages on the Web. If the j-th page is a dead end – i.e., it has no out-links –
then we assign a uniform probability of 1/n to all the elements in its column. Most of the elements of A are equal
to δ, the probability of jumping from one page to another without following a link. If n = 4 · 109
and p = 0.85,
then δ = 3.75 · 10−11. Matrix A is the transition probability matrix of the Markov chain. Its elements are all strictly
between zero and one and its column sums are all equal to one. An important result from matrix theory, known as the
Perron-Frobenius theorem, applies to such matrices: a nonzero solution of the equation
x = Ax (3)
1This project is originally based on a SIAM book chapter from Numerical Computing with Matlab by Clever B. Moler.
2
see ”Estimating search engine index size variability” by van den Bosch et al., Scientometrics,
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4833824/
1
Numerical Computing, Fall Semester 2022
Lecturer: Prof. O. Schenk
Assistants: E. Vecchi, L. Gaedke-Merzhauser ¨
exists and is unique to within a scaling factor. If this scaling factor is chosen so that
X
i
xi = 1, (4)
then x is the state vector of the Markov chain and corresponds to Google’s PageRank. The elements of x are all
positive and less than one. Vector x is the solution to the singular, homogeneous linear system
(I − A)x = 0 (5)
For modest n, an easy way to compute x in Julia is to start with some approximate solution, such as the PageRanks
from the previous month, or
x = ones(n, 1) / n
Then simply repeat the assignment statement
x = A * x
until the difference between successive vectors is within a specified tolerance. This is known as the power method
and is about the only possible approach for very large n. In practice, the matrices G and A are never actually formed.
One step of the power method would be done by one pass over a database of Web pages, updating weighted reference
counts generated by the hyperlinks between pages. The best way to compute PageRank in Julia is to take advantage
of the particular structure of the Markov matrix. Here is an approach that preserves the sparsity of G. The transition
matrix can be written
A = pGD + ez|
, (6)
where D is the diagonal matrix formed from the reciprocals of the out-degrees,
djj =
1/cj cj 6= 0
0 cj = 0.
(7)
e is the n-vector of all ones, and z is the vector with components
zj =
δ cj 6= 0
1/n cj = 0.
(8)
The rank-one matrix ez|
accounts for the random choices of Web pages that do not follow links. The equation
x = Ax (9)
can be written as
(I − pGD)x = γe (10)
where
γ = z
|x. (11)
We do not know the value of γ because it depends on the unknown vector x, but we can temporarily take γ = 1. As
long as p is strictly less than one, the coefficient matrix I − pGD is nonsingular and the equation
(I − pGD)x = e (12)
can be solved for x. Then the resulting x can be rescaled so that
X
i
xi = 1. (13)
Notice that vector z is not involved in this calculation. The following Julia statements implement this approach.
2
Numerical Computing, Fall Semester 2022
Lecturer: Prof. O. Schenk
Assistants: E. Vecchi, L. Gaedke-Merzhauser ¨
c = vec(sum(G, dims=1));
k = findall(!iszero, c);
D = sparse(k, k, map(x -> 1/x, c[k]), n, n);
e = ones(n, 1);
sI = sparse(I, n, n);
x = (sI – p * G * D) \ e;
x = x/sum(x);
The power method [1, 2] can also be implemented in a way that does not actually form the Markov matrix and so
preserves sparsity. Compute
G = p * G * D;
z = ((1- p) * map(!iszero, c) + map(iszero, c)) / n;
z = reshape(z, (1,500));
Start with
x = e / n;
Then repeat the statement
x = G * x + e * (z * x);
until x settles down to several decimal places.
It is also possible to use an algorithm known as inverse iteration [1, 3] where one defines
x = e / n;
A = p * G * D + e * z;
and then the following statement is repeated until convergence given a suitable value α.
x = (alpha * I – A) \ x;
x = x / sum(x);
For details of when and why this method works, please refer to the references cited above. However, we would like
to briefly observe what happens when we take α = 1 in the expression (αI − A). Since the resulting matrix (I − A)
is theoretically singular, with exact computation some diagonal elements of the upper triangular factor of (I − A)
would, in principle, be zero and this computation should fail. But with roundoff error, the computed matrix (I −A) is
probably not exactly singular. Even if it is singular, roundoff during Gaussian elimination will most likely prevent any
exact zero diagonal elements. We know that Gaussian elimination with partial pivoting always produces a solution
with a small residual, relative to the computed solution, even if the matrix is badly conditioned. The vector obtained
with the backslash operation, (I − A)\x, usually has very large components. If it is rescaled by its sum, the residual
is scaled by the same factor and becomes very small. Consequently, the two vectors x and Ax equal each other to
within roundoff error. In this setting, solving the singular system with Gaussian elimination blows up, but it blows up
in exactly the right direction.
Figure 1 is the graph for a tiny example, with n = 6 instead of n = 4 · 109
. The pages on the Web are identified by
strings known as uniform resource locators, or URLs. Most URLs begin with https because they use the hypertext
transfer protocol. In Julia, we can store the URLs as an array of strings in a cell array. The example in Figure 1 can
be written as follows by using a 6-by-1 cell array:
U = [“https://www.alpha.com”
“https://www.beta.com”
“https://www.gamma.com”
“https://www.delta.com”
“https://www.rho.com”
“https://www.sigma.com”]
In order to access the string contained in a cell we can simply write U[k], with k = 1, . . . , 6 in this small web graph.
Thus, U[1] would, e.g., correspond to the string ”https://www.alpha.com”.
We can generate the connectivity matrix by specifying the pairs of indices (i, j) of the nonzero elements. Because
there is a link to beta.com from alpha.com, the (2, 1) element of G is nonzero. The nine connections are described by
3
Numerical Computing,
Figure 1: A tiny web graph.
i = [ 2 6 3 4 4 5 6 1 1]
j = [ 1 1 2 2 3 3 3 4 6]
A sparse matrix is stored in a data structure that requires memory only for the nonzero elements and their indices.
This is hardly necessary for a 6-by-6 matrix with only 27 zero entries, but it becomes crucially important for larger
problems. The statements
n = 6;
G = sparse(vec(i), vec(j), vec(ones(Bool, 1, 9)));
Matrix(G)
generate the sparse representation of an n-by-n matrix with ones in the positions specified by the vectors i and j and
display its full representation.
0 0 0 1 0 1
1 0 0 0 0 0
0 1 0 0 0 0
0 1 1 0 0 0
0 0 1 0 0 0
1 0 1 0 0 0
The statement
c = sum(G, dims=1)
computes the column sums
2 2 3 1 0 1
Notice that c(5) = 0 because the 5th page, labeled rho, has no out-links. The statements
x = (sI – p * G * D) \ e;
x = x / sum(x);
solve the sparse linear system to produce
julia> x
0.3210
0.1705
0.1066
0.1368
0.0643
0.2007
4
Numerical Computing, Fall Semester 2022
Lecturer: Prof. O. Schenk
Assistants: E. Vecchi, L. Gaedke-Merzhauser ¨
Figure 2: Page Rank for the tiny web graph
The bar graph of x is shown in Figure 2. If the URLs are sorted in PageRank order and listed along with their in- and
out-degrees, the result is
index page-rank in out url
1 0.321017 2 2 https://www.alpha.com
6 0.200744 2 1 https://www.sigma.com
2 0.170543 1 2 https://www.beta.com
4 0.136793 2 1 https://www.delta.com
3 0.106592 1 3 https://www.gamma.com
5 0.0643118 1 0 https://www.rho.com
We see that alpha has a higher PageRank than delta or sigma, even though they all have the same number of
in-links. A random surfer will visit alpha over 32% of the time and rho only about 6% of the time.
For this tiny example with p = .85, the smallest element of the Markov transition matrix is δ = .15/6 = .0250.
A =
0.0250 0.0250 0.0250 0.8750 0.1667 0.8750
0.4500 0.0250 0.0250 0.0250 0.1667 0.0250
0.0250 0.4500 0.0250 0.0250 0.1667 0.0250
0.0250 0.4500 0.3083 0.0250 0.1667 0.0250
0.0250 0.0250 0.3083 0.0250 0.1667 0.0250
0.4500 0.0250 0.3083 0.0250 0.1667 0.0250
We can notice that all the columns of matrix A sum to one, thus agreeing with the definition of left stochastic matrix.
This mini-project includes the Julia file surfer.jl. A statement like
U, G = surfer(’https://www.xxx.zzz’,n)
starts at a specified URL and tries to surf the Web until it has visited n pages. If successful, it returns an n-by-1 cell
array of URLs and an n-by-n sparse connectivity matrix. Surfing the Web automatically is a dangerous undertaking
and this function must be used with care. Some URLs contain typographical errors and illegal characters. There is a
list of URLs to avoid that includes .gif files and Web sites known to cause difficulties. Most importantly, surfer can
get completely bogged down trying to read a page from a site that appears to be responding, but that never delivers the
complete page. When this happens, it may be necessary to have the computer’s operating system ruthlessly terminate
Julia. With these precautions in mind, you can use surfer to generate your own PageRank examples. The statement
U, G = surfer(’https://www.inf.ethz.ch’, 500);
5
Numerical Computing, Fall Semester 2022
Lecturer: Prof. O. Schenk
Assistants: E. Vecchi, L. Gaedke-Merzhauser ¨
accesses the home page of the Department of Computer Science at the Swiss Federal Institute of Technology in Zurich
(ETH) and generates a 500-by-500 test case. The file ETH500.mat included in the code folder was generated by
using this function in September 2014. In the following, it is indicated as the ETH500 data set. The statement
spy(G)
produces a spy plot (Figure 3) that shows the nonzero structure of the connectivity matrix. The statement
pagerank(U, G)
computes the page ranks, produces a bar graph (Figure 4) of the ranks, and prints the most highly ranked URLs in
PageRank order. The highly ranked pages are
julia> pagerank(U,G)
index page-rank in out url
57 0.145755 292 1 http://www.zope.org
466 0.049996 19 1 http://purl.org/dc/elements/1.1
39 0.028766 311 0 http://www.ethz.ch
58 0.021863 291 0 http://www.infrae.com
53 0.021202 288 0 http://www.cd.ethz.ch/services/web
101 0.016601 234 6 http://www.hk.ethz.ch/index_EN
72 0.016473 235 0 http://www.ethz.ch/index_EN
486 0.007374 14 2 http://www.handelszeitung.ch
463 0.006896 17 0 http://ns.adobe.com/xap/1.0
Figure 3: Spy plot of the Department of Computer Science (ETH Zurich) web graph.
Solve the following problems:
Preliminary: Read “A First Course on Numerical Methods”
Read Chapter 8 [1] carefully, in order to gain a better understanding of the topic summarized above.
1. Theory [20 points]
(a) What are an eigenvector, an eigenvalue and an eigenbasis?
(b) What assumptions should be made to guarantee convergence of the power method?
6
Numerical Computing, Fall Semester 2022
Lecturer: Prof. O. Schenk
Assistants: E. Vecchi, L. Gaedke-Merzhauser ¨
Figure 4: PageRank of the Department of Computer Science (ETH Zurich) web graph.
(c) What is the shift and invert approach?
(d) What is the difference in cost of a single iteration of the power method, compared to the inverse iteration?
(e) What is a Rayleigh quotient and how can it be used for eigenvalue computations?
2. Other webgraphs [10 points]
Use surfer.jl and pagerank.jl to compute PageRanks for some subset of the Web of your choice. Do you
see any interesting structure in the results (e.g., cliques, see next question)? Report on three PageRanks for different
webgraphs by showing the connectivity matrix, the PageRanks, and the ten most important entries in the graph.
3. Connectivity matrix and subcliques [5 points]
The connectivity matrix for the ETH500 data set (ETH500.mat) has various small, almost entirely nonzero, submatrices that produce dense patches near the diagonal of the spy plot. You can use the zoom button to find their indices.
The first submatrix has, e.g., indices around 80. Mathematically, a graph where all nodes are connected to each other
is known as a clique. Identify the organizations within the ETH community that are responsible for these near cliques.
4. Connectivity matrix and disjoint subgraphs [10 points]
Figure 5 reports the graph of a tiny six-node subset of the Web. Unlike the example reported in Figure 1, in this case
there are two disjoint subgraphs.
1. What is the connectivity matrix G? Which are its entries?
2. What are the PageRanks if the hyperlink transition probability p assumes the default value of 0.85?
3. Describe what happens with this example to both the definition of PageRank and the computation done by
pagerank in the limit p → 1.
5. PageRanks by solving a sparse linear system [40 points]
The function pagerank(U,G) computes PageRanks by solving a sparse linear system. It then plots a bar graph and
prints the dominant URLs.
7
Numerical Computing, Fall Semester 2022
Lecturer: Prof. O. Schenk
Assistants: E. Vecchi, L. Gaedke-Merzhauser ¨
Figure 5: Another tiny Web.
1. Create pagerank1.jl by modifying pagerank.jl to use the power method instead of solving the sparse
linear system. The key statements are
G = p * G * D;
z = ((1- p) * map(!iszero, c) + map(iszero, c)) / n;
z = reshape(z, (1,500));
while termination_test
x = G * x + e * (z * x);
end
What is an appropriate test for terminating the power iteration?
2. Create pagerank2.jl by modifying pagerank.jl to use the inverse iteration. The key statements are
while termination_test
x = (alpha * I – A) \ x;
x = x / sum(x);
end
Use your functions pagerank1.jl and pagerank2.jl (set α = 0.99) to compute the PageRanks of the
six-node example presented in Figure 1. Make sure you get the same result from each of your three functions.
3. We now want to analyse the impact of α on the inverse iteration. Using the ETH500 example, set α equal
to 0.8, 0.9, 0.95 and 1. Comment on the different number of iterations the four cases take until convergence.
Analyse your results and explain what you observe.
Hint: Check your solution x for all 4 cases. Are they always the same?
4. Use your functions pagerank1.jl and pagerank2.jl (set α = 0.99) to compute the PageRanks of the
three selected graphs from exercise 2. Report on the convergence of the two methods for these subgraphs and
summarize the advantages and disadvantages of the power method implemented in pagerank1.jl against
the inverse iteration in pagerank2.jl.
8
Numerical Computing, Fall Semester 2022
Lecturer: Prof. O. Schenk
Assistants: E. Vecchi, L. Gaedke-Merzhauser ¨
Quality of the code and of the report [15 Points]
The highest possible score of each project is 85 points and up to 15 additional points can be awarded based on the
quality of your report and code (maximum possible grade: 100 points). Your report should be a coherent document,
structured according to the template provided on iCorsi. If there are theoretical questions, provide a complete and
detailed answer. All figures must have a caption and must be of sufficient quality (include them either as .eps or .pdf).
If you made a particular choice in your implementation that might be out of the ordinary, clearly explain it in the
report. The code you submit must be self-contained and executable, and must include the set-up for all the results that
you obtained and listed in your report. It has to be readable and, if particularly complicated, well-commented.
Additional notes and submission details
Summarize your results and experiments for all exercises by writing an extended LaTeX report, by using the template provided on iCorsi (https://www.icorsi.ch/course/view.php?id=14666). Submit your gzipped
archive file (tar, zip, etc.) – named project number lastname firstname.zip/tgz – on iCorsi (see [NC
2022] Project 1 – Submission PageRank Algorithm) before the deadline. Submission by email or through any other
channel will not be considered. Late submissions will not be graded and will result in a score of 0 points for that
project. You are allowed to discuss all questions with anyone you like, but: (i) your submission must list anyone you
discussed problems with and (ii) you must write up your submission independently. Please remember that plagiarism
will result in a harsh penalization for all involved parties.
In-class assistance
If you experience difficulties in solving the problems above, please join us in class either on Tuesdays 16:30-18:15 in
room C1.03 or on Wednesdays 14:30-16:15 in room D1.15.
References
[1] The power method and variants, Chapter 8: Eigenvalues and Singular Values, SIAM Book “A First Course on
Numerical Methods”, C. Greif, U. Ascher, pp. 219-229.
[2] Power iteration, http://en.wikipedia.org/wiki/Power_iteration
[3] Inverse iteration, http://en.wikipedia.org/wiki/Inverse_iteration
9
