Large-scale linear algebra: Dense matrix methods¶

The basic problem with most of the linear algebra techniques we have learned so far is that they scale badly for large matrices. Ordinary Gaussian elimination (LU factorization), Gram–Schmidt and other QR factorization algorithms, and techniques that computes all the eigenvalues and eigenvectors, all require $\sim n^3$ operations and $\sim n^2$ storage for $n \times n$ matrices.

This all assumes that you explicitly store and compute with all of the entries of the matrix, regardless of their values. Hence, they are sometimes called dense matrix algorithms (the opposite of "sparse" matrices, discussed below).

So, doubling the size of the matrix asymptotically requires about 8× more time. For any finite $n$, it is not quite a factor of 8 because computers are complicated; e.g. for larger matrices, it can use multiple processors more efficiently:

Still, if we take the $O(n^3)$ scaling as a rough guide, this would suggest that LU-factorizing (lufact) a $10^6 \times 10^6$ matrix would take $0.02\mbox{sec} \times 1000^3 \sim \mbox{months}$ and finding the eigenvectors and eigenvalues (eigfact) would take $2 \mbox{sec} \times 1000^3 \sim \mbox{decades}$.

In practice, we actually usually run out of space before we run out of time. If we have 16GB of memory, the biggest matrix we can store (each number requires 8 bytes) is $8n^2\mbox{ bytes} = 16\times 10^9 \implies 40000 \times 40000$.

Sparse Matrices¶

The saving grace is that most really large matrices are sparse = mostly zeros (or have some other special structure with similar consequences). You only have to store the nonzero entries, and you can multiply matrix × vector quickly (you can skip the zeros).

In Julia, there are many functions to work with sparse matrices by only storing the nonzero elements. The simplest one is the sparse function. Given a matrix $A$, the sparse(A) function creates a special data structure that only stores the nonzero elements:

(Of course, in practice you would want to create the sparse matrix directly, rather than first making the "dense" matrix A and then converting it to a sparse data structure.)

We've actually seen this several times in graph/network-based problems, where we often get matrices of the form:

where D is diagonal (very sparse!) and G is the incidence matrix. Since each graph node is typically only connected to a few other nodes, G is sparse and so is A.

If each node is connected to a bounded number of other nodes (say, ≤ 20), then A only has $\sim n$ (i.e. proportional to n, not equal to n) entries, and $Ax$ can be computed in $\sim n$ operations and $\sim n$ storage (unlike $\sim n^2$ for a general matrix).

So, a $10^6 \times 10^6$ sparse matrix might be stored in only a few megabytes and take only a few milliseconds to multiply by a vector.

Much of large-scale linear algebra is about devising techniques to exploit sparsity or any case where matrix × vector is faster than n².

Scalable computation of eigenvalues¶

In fact, we already learned one algorithm that works well for sparse matrices: the power method to compute eigenvalues and eigenvectors. If we just repeatedly multiply a random vector by $A$, it converges towards the eigenvector of the largest $|\lambda|$. And since this only involves matrix × vector operations, it can take advantage of sparse matrices.

Moreover, there are variants of this algorithm that work for the smallest eigenvalues as well, and it turns out that there are more sophisticated variants that converge even faster than power iterations. In Julia, these are provided by the eigs function, which lets you compute a few of the biggest or smallest eigenvalues quickly even for huge sparse matrices.

Example¶

As an example, let's consider the two-dimensional grid of masses and springs that I showed in an earlier lecture, whose eigenvectors are the vibrating modes.

This can be thought of as a discretized approximation of a vibrating drum, which is described by the partial differential equation $\nabla^2 h = \frac{\partial^2 h}{\partial t^2}$ where $h(x,y,t)$ is the height of the drum surface (= zero at the edges of the drum). (This is an example taken from the class 18.303: Linear Partial Differential Equations at MIT.)

For this class, don't worry too much about how the matrix is constructed.

The above code defines the matrix for a box-shaped drum, but for fun I will change it to define a drum over an oddly shaped domain:

This all eventually leads to the following matrix, whose eigenvalues $ \lambda = \omega^2$ are the squares of the frequencies and whose eigenvectors are the vibrating modes:

This is a $60000\times60000$ matrix, which would too big to even store on my laptop if we stored every entry. Because it is sparse, however, almost all of the entries are zero and we only need to store those.

The nnz function computes the number of nonzero entries, and we can use it to compute the fraction of nonzero entries in A:

Less than 0.01% of the entries are nonzero! It really pays to take advantage of this.

Now we'll compute a few of the smallest-|λ| eigenvectors using eigs. We'll use the Interact package to interactively decide which eigenvalue to plot.

Notice that it took less than two seconds to solve for 20 eigenvectors and eigenvalues!

This is because eigs is essentially using an algorithm like the power method, that only uses repeated multiplication by $A$.

Or, sometimes, particularly to find the smallest $|\lambda|$ eigenvectors, it might repeatedly divide by $A$, i.e. solve $Ax=b$ for $b = A^{-1} x$. But it can't actually compute the inverse matrix, and I said that LU factorization was $\sim n^3$ in general. So, what is happening?

Sparse-direct solvers for Ax=b¶

Even if $A$ is a sparse matrix, $A^{-1}$ is generally not sparse. However, if you arrange things cleverly, often the $L$ and $U$ factors are still sparse!

This leads to something called sparse-direct solvers: they solve $Ax=b$ by ordinary Gaussian elimination to find $A = LU$, but they take advantage of sparse $A$ to avoid computing with zeros. Moreover, they first re-order the rows and columns of $A$ so that elimination won't introduce too many zeros — this is a tricky problem that mostly involves graph theory, so I won't try to explain it here.

Julia (and Matlab) both use sparse-direct algorithms automatically when you do A \ b if $A$ is stored as a sparse matrix. When they work (i.e. when the L and U factors are sparse), these algorithms are great: fast, memory-efficient, reliable, and worry-free "black boxes".

Example: Helmholtz solver¶

The following code solves a scalar Helmholtz equation

This equation describes the propagation of waves u from a source f at a frequency $\omega$ is the frequency. For example, imagine water waves travelling across a shallow pond, with $u(x,y)$ being the height of the wave, where $f$ represents an vibrating disturbance that creates the wave.

We discretize this into a matrix equation $Au=f$ by discretizing space $(x,y)$ into a grid and approximating derivatives $-\nabla^2$ by differences on the grid (this is an FDFD method).

Again, don't worry too much about the details of this construction (take 18.303 to find out more). The important thing is that we will have a grid (graph!) of many unknowns, but the problem is sparse because each grid point only "talks" to its 4 nearest neighbors.

Let's set up a scattering problem with a cylindrical scatterer (a slightly slower wave speed, e.g. a different depth of water, inside a small cylindrical region, with a wavelength $2\pi/\omega$ of 1:

Again, this is a huge matrix: $230400 \times 230400$. But it is incredibly sparse, so solving it will be no problem:

For the right-hand side b, let's use a "point" source on one side. We'll use the reshape function to convert between 2d arrays and column vectors for solving with A. Note that the solution is complex, since it corresponds physically to an oscillating solution $u(x,y)e^{-i\omega t}$ and $u$ has a phase angle due to the absorbing boundary layers (which make A non-Hermitian); we'll just plot the real part.

We solved a $200000\times200000$ matrix problem in 2–3 seconds, and less than 1GB of memory. Pretty good!

Iterative solvers¶

Unfortunately, sparse-direct solvers like those we are using above have two limitations:

• They only work if the matrix is sparse. There are lots of problems where $A$ has some special structure that lets you compute $A*x$ quickly, e.g. by FFTs, and avoid storing the whole matrix, but for which $A$ is not sparse.

• They scale poorly if the sparse matrix comes from a 3d grid or mesh. For an $s$-element 1d mesh with $n=s$ degrees of freedom, they have $O(s)$ complexity. For an $s \times s$ 2d mesh with $n=s^2$ degrees of freedom, they take $O(n \log n)$ operations and require $O(n)$ storage. But for a 3d $s\times s\times s$ mesh with $n = s^3$, they take $O(n^2)$ operations and require $O(n^{4/3})$ storage (and you often run out of storage before you run out of time).

The alternative is an iterative solver, in which you supply an initial guess for the solution $x$ (often just $x=0$) and then it iteratively improves the guess, converging (hopefully) to the solution $A^{-1} b$, while using only matrix-vector operations $Ax$.

Iterative solvers are the method of choice (or, more accurately, of necessity) for the very largest problems, but they have their downsides. There are many iterative solver algorithms, and you have to know a little bit to pick the best one. They may not converge at all for non-symmetric $A$, and in any case may converge very slowly, unless you provide a "magic" matrix called a preconditioner that is specific to your problem. (It is often a research problem in itself to find a good preconditioner!)

Toy example: Steepest-descent algorithm¶

If $A$ is a real-symmetric positive-definite matrix, then solving $Ax = b$ is equivalent to minimizing the function:

Just compute $\nabla f = \cdots = 2(Ax - b)$, which equals zero at the minimum. The definiteness of $A$ means that the function $f$ is convex, so there is exactly one global minimum. Another way to see this: if $Az=b$, then $f(z+v) - f(z) = \cdots = v^T A v > 0$ for any vector $v \ne 0$, so $f(z)$ must be the minimum.

One of the simplest iterative algorithms is just to go downhill: minimize $f(x + \alpha d)$ over $\alpha$, where $d$ is the downhill direction $-\frac{1}{2}\nabla f = b - Ax = r$, where $r$ is called the residual. We can perform this line minimization analytically for this $f$, for an arbitrary $d$, to find $\alpha = d^T r / d^T A d$.

The steepest-descent algorithm simply performs this downhill line-minimization repeatedly, starting at an initial guess x (typically just x=0), e.g. stopping when the norm of the residual is less than some tolerance times the norm of b.

To see what's going on, let's try a $2\times2$ matrix where we can easily visualize the progress.

The solution "zig-zags" down the long, narrow valley defined by the quadratic function f. This is a common problem of steepest-descent algorithms: they tend to go towards the center of valleys (down the "steep" direction), rather than along the valleys towards the solution.

To fix this problem, basically we need to implement some kind of "memory": it has to "remember" that it just "zigged" in order to avoid "zagging" back where it came from.

From steepest-descent to conjugate-gradient¶

The most famous way to improve steepest descent with "memory" is the conjugate-gradient algorithm. I won't explain it here (Shewchuk's article is a good introduction to its relationship to steepest descent), but the implementation ends up being almost identical to steepest descent. However, instead of setting the line-search direction equal to the downhill direction r, the line-search direction is instead a linear combination of r with the previous search direction:

After some initial slow progress, the conjugate-gradient algorithm quickly zooms straight to the solution:

In our 2x2 problem, you can see that it actually converges in two steps:

You don't have to write your own iterative solvers¶

There are several packages out there with iterative solvers that you can use, e.g. the IterativeSolvers package:

PETSc¶

PETSc.jl provides a Julia interface to the PETSc library: large-scale iterative and sparse solvers for distributed-memory parallel systems. It is a bit harder to use, because it expects you to set up its own kind of sparse matrix that works on distributed-memory systems.

Preconditioners¶

Most iterative solvers are greatly accelerated if you can provide a preconditioner $P$: roughly, an approximate inverse of $A$ that is easy to compute. The preconditioner is applied at every step of the iteration in order to speed up convergence. (For example, you could solve $PAx = Pb$ instead of $Ax=b$; this is called a "left" preconditioner.)

For example, let's consider a problem where the matrix $A = L + S$ is a sum of the symmetric-tridiagonal discrete Laplacian $L$ (from above) and a small, sparse perturbation $S$. As our preconditioner, we'll simply use $P = L^{-1}$, since this is a good approximation for $A$ and L \ b is fast (linear time):

Our home-brewed CG function above does not accept a preconditioner, but the IterativeSolvers package cg function does, and it makes a huge difference in the convergence:

If you can find a good preconditioner, you can often speed things up by orders of magnitude. Unfortunately, finding preconditioners is hard and problem-dependent. There are some general "recipes" for things to try, but there are many problems (like the scalar-Helmholtz problem above) where good preconditioners are still an open research problem.

What makes a good preconditioner matrix $P$ for a matrix $A$? For the "ideal" preconditioner $P = A^{-1}$, you would have $PA = I$, but of course this is impractical: if you could compute $A^{-1}$ quickly, you wouldn't need an iterative solver. So, what you want is for $PA$ to be "like" $I$ in some sense. The eigenvalues of $I$ are all 1, and it turns out that a good preconditioner makes $PA$ have eigenvalues that are mostly clustered together.

Let's see how the eigenvalues of $A$ and $PA = L^{-1} A$ compare in this case.

You can see from the plot that $A$ has eigenvalues all over the map (from $10^{-3}$ to 10), but $L^{-1} A$ has eigenvalues that are mostly 1 (like $I$) with a handful of outliers.

Key points¶

• Large matrix problems ($> 10^4 \times 10^4$) almost have a matrix $A$ with a special structure that allows you to multiply Ax quickly. Most commonly, $A$ is sparse (mostly zero).

• Many specialized methods are available to take advantage of this special structure in order to solve problems much faster than the n³ time and n² memory of "dense" solvers.

  • For sparse matrices, A \ b can actually be solved directly by ordinary Gaussian elimination, cleverly re-ordering the unknowns so that the LU factors remain sparse. These "sparse direct" methods are easy to use and widely available.

  • For the very largest matrices, for eigenproblems, or non-sparse matrices with a fast Ax, the alternative is an iterative method: start with guess for the solution, and do some clever operations that rapidly improve the guess using only matrix-times-vector.

  • In many problems, you can speed up iterative solvers by orders of magnitude if you have an approximate solution for Ax=b, even a very "bad" approximation that only uses part of th matrix. Finding such preconditioners is often more of an art than a science.

• Fully understanding the methods for large sparse problems goes far beyond the scope of 18.06. Better coverage is given in 18.335. But it is important to know that such methods are out there, even if you don't understand the detail, so that you know where to begin looking if you encounter a huge matrix in practice.