LU Decomposition and the PageRank Algorithm

Preamble: Run the cells below to import the necessary Python packages

Revisiting the coauthorship graph¶

We previously considered the coauthorship graph among physicists based arXiv postings in astro-ph. The graph contains $N = 18772$ nodes, which correspond to unique authors observed over the period 1993 -- 2003. If Author i and Author j coauthored a paper during that period, their nodes are connected. In order to analyze this large graph, we will downsample it to a smaller graph with $N = 1000$ nodes representing the most highly-connected authors.

This dataset is from the Stanford SNAP database

Estimating degree distributions, centrality, and the PageRank algorithm¶

We previously saw that the degree distribution determines how often a random walk on the graph will visit a given node. This can be seen an an indicator of how “important” a given node is to the structure of the graph. If we wanted to rank the authors purely based on their degree distribution, we could write this ranking as

where $\mathbb{1} \in \mathbb{R}^N$ is a vector of ones, $A \in \mathbb{R}^{N\times N}$ is our adjacency matrix, and $\mathbf{R} \in \mathbb{R}^N$ is our desired ranking. Notice how the matrix product with the ones vector allows us to compute the rowwise sum. This makes sense, because we expect that this operation should take $N^2$ operations, and matrix-vector products are $N^2$.

How else might we group authors?¶

The degree ranking simply tells use who are the most well-connected authors are. They are the individuals who directly co-authored a paper with someone else at any point during the period over which the arXiv dataset was collected.

Recall that we can think of the degree distribution as the long-time limit of a random walker that “hops” among the neighboring nodes in a graph. What if we add stochasticity? What if we allow the random walker to sometimes jump to any new node, not connected to the current node? This will cause the walker to explore the graph more broadly, and discover authors who are less connected to the dominant central cluster. We will draw a loose analogy with the Ornstein-Uhlenbeck process. The spring constant pulls a random walker back to the origin, but the random noise allows the walker to explore the space more broadly. Likewise, the degree distribution pulls the random walker to the most connected nodes, but random jumps allow the walker to explore the graph more broadly.

This is the idea behind the PageRank algorithm, which is how Google originally ranked web pages. In Google’s case, the entire internet is too large to know a priori the adjacency matrix of all web pages. Instead, they use a random walk algorithm to estimate the degree distribution of the web pages. Random hops allow the algorithm to more easily find smaller subnetworks, so that it doesn’t keep revisiting dominant nodes like Wikipedia or Facebook. This is sometimes refered to as the “random surfer model” of PageRank.

Implementing the PageRank algorithm¶

We are now going to implement the PageRank algorithm on the coauthorship graph. The PageRank algorithm is a variant of the random walk algorithm that allows the random walker to sometimes jump to a new node with probability $1 - \alpha$.

If we interpret our desired ranking $\mathbf{R}$ as a probability distribution, then we can write a master equation describing the dynamics of the probability distribution of many random walkers. The vector $\mathbf{R}_t$ represents the unnormalized probability distribution of the random walkers at time $t$, and so we refer to it as the “rank” because its relative values, but not absolute values, represent the relative importance of a given node.

where $\mathbb{1} \in \mathbb{R}^N$ is a vector of ones, $A \in \mathbb{R}^{N \times N}$ is the adjacency matrix of the graph, and $\alpha$ is a parameter that controls the probability of the random surfer jumping to a random node. The first term on the right-hand side represents the process of the walker moving to a neighboring node, while the second term represents the process of jumping to a new node.

Solving this equation for the steady state distribution $\mathbf{R}_t = \mathbf{R}_{t+1}$ results in the PageRank equation

This equation can be rewritten as an equation for $\mathbf{R}$,

where $I$ is the identity matrix. We write these equations as a single matrix equation

where $\mathbf{M} \equiv (1 - \alpha) (I - \alpha A)$. If we set $\alpha = 1$ (no random jumps), then this equation becomes

Inverting a matrix¶

So far, we have taken for granted that we can easily invert a large matrix $\mathbf{A}$. Given a linear problem

with $A \in \mathbb{R}^{N \times N}$, $\mathbf{x} \in \mathbb{R}^N$, and $\mathbf{b} \in \mathbb{R}^N$, we can solve for $\mathbf{x}$ by inverting $A$:

How would we do this by hand?¶

Gauss-Jordan elimination: perform a series of operations on both the left and the right hand sides of our equation (changing both $A$ and $\mathbf{b}$), until the left side of the equation is triangular in form. Because Gauss-Jordan elimination consists of a series of row-wise addition and subtraction operations, we can think of the process as writing each row vector in terms of some linear combination of row vectors, making the whole process expressible as a matrix $M$

where we have defined the upper-triangular matrix $U \equiv MA$. This confirms our intuition that the time complexity of Gauss-Jordan elimination is $\sim\mathcal{O}(N^3)$, since that’s the time complexity of multiplying two $N \times N$ matrices.

Why do we want triangular matrices?¶

If we can reduce our matrix $A$ to upper-triangular form, then we can quickly solve for $\mathbf{x}$ using forward substitution. This should take $\sim\mathcal{O(N^2)}$ operations if $A \in \mathbb{R}^{N \times N}$

We can implement a forward substitution algorithm that takes advantage of the fact that the matrix is triangular. To simplify our notation, we will assum that $A$ is lower-triangular, so that the first row has only one non-zero element at $A_{0,0}$. We can theniterate over rows of $A$ in order to find the elements of the solution vector $\mathbf{x} \in \mathbb{R}^N$ one-by-one. For example, the first element is simplygiven by

The second element is given by

and so on. We can therefore write a general algorithm for forward substitution as

Notice how the number of operations scales with the number of non-zero elements in the row. We can apply the same approach to an upper-triangular matrix, but we would need to iterate over the rows in reverse order.

Looking at our implementation of the triangular solve, we can see that it has a single for loop over the rows of the matrix. However, within each row, it performs a dot product. In the first row, this dot product reduces to a single scalar multiplication, but in the final row, it reduces to a dot product of $N$ elements. Across the $N$ steps of the for loop, the number of product-sum operations, on average, is proportional to $N$.

Across the explicit outer loop and the vectorized dot product, we therefore expect to perform $\sim\mathcal{O}(N^2)$ operations total. Another way to see this is that the number of operations is proportional to the number of non-zero elements in the matrix, which is $N(N-1)/2$ for a lower-triangular matrix.

Let’s see if this scaling holds in practice by using the timeit module to see how the runtime of our implementation scales with the matrix size. We will compare our implementation to the built-in numpy solver.

The scaling does not match our expectation. We can see that the empirical scaling almost looks like $\sim\mathcal{O}(N)$ for our triangular solver. But we know that this is physically impossible, because we have to access $N(N-1)/2$ non-zero elements in $A$ to solve the system.

What happened? Timescale separation between our nested $\mathcal{O}(N)$ operations. The outer for loop is written in Python, which is slow, but the inner dot product is vectorized, and so it uses the highly optimized BLAS library. As a result, the inner loop is much faster than the outer loop, and so we don’t see the quadratic scaling at these matrix sizes. In fact, in order to see the quadratic scaling, we would need to access matrix sizes that are so large that we would run out of memory on a typical computer. So, in this case, we are not able to perform numerical experiments that reach the true asymptotic scaling regime.

What about numpy’s built-in solver? It performs way faster than our solver at small matrix sizes, but it appears to gets much worse as the matrix gets larger. This is because, unlike our implementation, numpy’s solver applies to any square matrix. It does not know in advance that our matrix has exploitable triangular structure, and the time cost to check for this structure before solving is substantial enough that it is better for to use a general solver instead. General matrix inversion has the same time complexity as Gauss-Jordan elimination or matrix-matrix multiplication, which is $\sim\mathcal{O}(N^3)$. However, various advanced algorithms have been developed that drive down the prefactor.

Exploiting sparsity. We can see that triangular solvers are a specific case of a common phenomenon in numerical linear algebra. If a matrix has particular structure, particularly sparsity, the time cost of many operations can be reduced to be proportional to the number of nonzero elements. Since a triangular matrix has $\mathcal{O}(N^2)$ nonzero elements, it’s worth exploiting this structure when performing a $\mathcal{O}(N^3)$ operation like matrix inversion.

To invert a full matrix, we first need reach triangular form¶

We know that we want to modify our full matrix equation

using a Gauss-Jordan elimination procedure, which we know has the form of a matrix multiplication

We define $U \equiv M A$, and we know that if can reduce $U$ to upper-triangular form, then we can solve for $\mathbf{x}$ using forward substitution in $\sim\mathcal{O}(N^2)$ operations.

But what is the form of $M$?

Our key insight is that the Gaussian elimination matrix $M$ for a full matrix inversion turns is a triangular matrix. Moreover, because the inverse of a triangular matrix is also a triangular matrix (this can be shown by writing out the algebraic terms in the matrix multiplication algebra for $U^{-1} U = I$), we therefore write a transformed version of our problem

defining $M \equiv L^{-1}$, we arrive at the celebrated LU matrix factorization

where $L$ is a lower triangular matrix, and $U$ is an upper triangular matrix.

We can therefore solve for $\mathbf{x}$ by introducing an intermediate variable $\mathbf{h} \equiv U \mathbf{x}$, and then $\mathbf{x}$

If we can find the LU form, exactly solving the linear problem consists of two steps, each of which has a time complexity of $\sim\mathcal{O}(N^2)$:

1. Solve the equation $\mathbf{h} = L^{-1} \mathbf{b}$. Since $L$ is a triangular matrix, this can be performed using substitution
2. Solve the equation $\mathbf{x} = U^{-1} \mathbf{h}$. Since $U$ is also triangular, this inverse can also be computed using substitution

The hard part, then, will be to find the LU form in the first place. This is equivalent to finding the Gaussian elimination matrix for the full matrix inversion problem.

About LU factorization¶

An iterative implementation of $LU$ factorization uses Gauss-Jordan elimination in a principled order. We first eliminate the first column of the matrix, then the second column, and so on. A more sophisticated algorithm, Crout’s method with pivoting, performs optimal order of decomposition based on the specific values present in a given row or column. However, our approach is matrix agnostic, and so we will not implement pivoting.

Decomposition into L and U is $\sim \mathcal{O}(N^3)$ for an $N \times N$ matrix, meaning that it has the same time complexity as matrix-matrix multiplication.

• Given a matrix $A$ and target $\mathbf{b}$, what is the overall runtime to find $\mathbf{x}$?

The LU decomposition algorithm¶

1. Set $L = I$ and $U = A$ (the input matrix)

2. For each column $j$ in $U$:

   1. Choose the pivot element $U_{jj}$, corresponding to the diagonal element in the $j^{th}$ column
   2. For each row $i$ below the pivot:
      1. Find the multiplier $L_{ij} = U_{ij} / U_{jj}$
      2. Subtract $L_{ij} \times$ the entire $j$ th row from the $i^{th}$ row in $U$
      3. Store $L_{ij}$ in the $i^{th}$ row of $L$

In the context of the LU decomposition algorithm, the ‘pivot’ refers to the element in the current column, denoted as $U_{jj}$, that we are processing in the outer loop. The pivot serves as a divisor to find the multipliers $L_{ij}$ which are then used to eliminate the elements below the pivot, making them zero to carve out an upper triangular matrix $U$.

Simultaneously, these multipliers are stored in the lower triangular matrix $L$. Essentially, we are utilizing the pivot to find scalar values that can help perform row operations to systematically zero out elements below the diagonal in $U$, while building up the $L$ matrix. We can see that there are three nested loops, which means that the time complexity of the LU decomposition algorithm is $\sim \mathcal{O}(N^3)$.

Applying LU Decomposition of orthogonal basis sets¶

A Hadamard matrix is a square matrix comprising a complete collection of length $N$ mutually-orthogonal vectors with all elements equal to either +1 or -1. We can convert any complete orthonormal basis to a Hadamard matrix by ordering the vectors systematically, and substituting 0 and 1 for $\pm 1$ with a recursive rule. The recursive rule to generate a $2^n \times 2^n$ Hadamard matrix is

where the base case is $H_1 = 1$. We can build this matrix bottom-up using dynamic programming, or top-down using recursion.

We can verify that the Hadamard matrix is orthogonal by computing its condition number, $\kappa(H) = \|H\| \|H^{-1}\|$.

Since all rows and columns have inner product zero with each other, the matrix is orthogonal, and so the condition number is the lowest possible value of 1. Recall that the condition number of a random $N \times N$ matrix is $\sim N$.

We can see that the Hadamard matrix has unique structure due to its self-similar construction. We can highlight this structure by decomposing the Hadamard matrix into its $L$ and $U$ factors.

We can see that the $L$ and $U$ matrices are both fractal matrices known as the Sierpinski triangle.

Applying PageRank to the coauthorship graph¶

Now that we have a method for inverting a matrix, we can now apply the PageRank algorithm to the coauthorship graph. Recall that the PageRank algorithm is a random walk algorithm that can be written as a matrix inversion problem:

where $\mathbf{M} \equiv (1 - \alpha) (I - \alpha A)$.

We can see that high values of $\alpha$ cause the PageRank algorithm to look more like the degree distribution. This emphasizes larger subnetworks and and more highly-connected nodes. Conversely, low values of $\alpha$ cause the PageRank algorithm to converge to something closer to a uniform distribution over all nodes. The effect of the damping factor $\alpha$ is to control the trade-off between the random walk limit, which samples the well-connected nodes better, and the limit in which all nodes are sampled with equal probability. The latter case helps sample rarer nodes more frequently, at the expense of accurately representing the degree distribution for higher-degree nodes.