CS 5220

Applications of Parallel Computers

Intro to Sparse LA

Please click the play button below.

World of Linear Algebra

Dense methods (last week)
Sparse direct methods (Thurs)
Iterative methods (today and Thurs)

Dense methods

% Dense (LAPACK)
[L,U] = lu(A);
x = U\(L\b);

Direct representation of matrices with simple data structures (no need for indexing data structure)
Mostly \(O(n^3)\) factorization algorithms

Before we sketch the new, though, let’s sketch the old.

Dense linear algebra is all about dealing with matrices where we have to store almost every entry explicitly, and we use simple data structures to do so. The typical algorithms in this space involve factoring matrices into products of simpler matrices. For example, we saw last week that we can write the Gaussian elimination method for solving linear systems as LU factorization, or factorization of the coefficient matrices into a permutation, a lower triangular matrix, and an upper triangular matrix. Then we can solve linear systems via forward and backward substitution (also known as lower and upper triangular system solves). Gaussian elimination is an O(n^3) algorithm, and this is true for most of the algorithms in this area – cubic work and quadratic data, which is pretty good for caches.

Software Strategies: Dense Case

Assuming you want to use (vs develop) dense LA code:

Learn enough to identify right algorithm
(e.g. is it symmetric? definite? banded? etc)
Learn high-level organizational ideas
Make sure you have a good BLAS
Call LAPACK/ScaLAPACK!
For \(n\) large: wait a while

If you are going to use dense linear algebra, you need to know a little bit about common structural properties so you can choose the right algorithm. And it helps to have some understanding of the high-level organization of the dense codes, so that you know about things like blocking, and can do some simple blocked method design by hand when your problem has special structure that merits it. This is the type of thing that I teach in CS 4220 and in CS 6210. But you can get a long way just knowing the names of a few standard factorizations.

Dense linear algebra is relatively easy in this way. There are good algorithms for these types of factorizations in the LAPACK library, which is used by languages like MATLAB, R, Julia, and NumPy/SciPy. LAPACK runs fast because we can make level-3 BLAS routines – such as matrix-matrix multiply – run fast. If we want to scale beyond the performance we can get from LAPACK with maybe some shared-memory parallelism in the BLAS, there are packages like MAGMA, PLASMA, and ScaLAPACK out there. But we will probably mostly use one of a small number of common libraries in this space.

Sparse direct methods

% Sparse direct (UMFPACK + COLAMD)
[L,U,P,Q] = lu(A);
x = Q*(U\(L\(P*b)));

Direct representation, keep only the nonzeros
Factorization costs depend on problem structure (1D cheap; 2D reasonable; 3D gets expensive; not easy to give a general rule, and NP hard to order for optimal sparsity)
Robust, but hard to scale to large 3D problems

In sparse matrices, we have relatively few nonzeros, and so we use a data structure that only explicitly keeps track of those nonzeros. We’ve already mentioned one such data structure, the compressed sparse row (or compressed sparse column) format. As you might imagine, there are others as well.

I want to for the moment distinguish between this type of explicit sparse representation and what we call “data sparse” matrices (which may have a lot of nonzeros, but can be written in terms of fewer than O(n^2) parameters). This is because explicit sparse representations are a natural setting for sparse direct methods. Sparse direct methods for solving linear systems involve LU factorization, just as in the dense case. But the sparse case is way more complicated. The factorization costs depend on the sparsity pattern of the matrix in a fairly subtle way. We usually talk about this in terms of the “graph” associated with a sparse matrix; we’ll talk about this later. Anyhow, two matrices with the same dimension and the same number of nonzero entries may be very different from the perspective of Gaussian elimination, with one factorization being much more expensive to compute than the other. For problems that come from 2D spatial things, like PDEs on a two-dimensional domain, a sparse direct factorization can often be the best tool to use. For 3D PDE problems, matters get much more complicated.

Sparse Direct Strategies

Assuming you want to use (vs develop) sparse LA code

Identify right algorithm (mainly Cholesky vs LU)
Get a good solver (often from list)
- You don’t want to roll your own!
Order your unknowns for sparsity
- Again, good to use someone else’s software!
For \(n\) large, 3D: get lots of memory and wait

Sparse iterations

% Sparse iterative (PCG + incomplete Cholesky)
tol = 1e-6;
maxit = 500;
R = cholinc(A,'0');
x = pcg(A,b,tol,maxit,R',R);

Only need \(y = Ax\) (maybe \(y = A^T x\))
Produce successively better (?) approximations
Good convergence depends on preconditioning
Best preconditioners are often hard to parallelize

Alas, for 3D PDE discretizations, problems from social networks, and a variety of other problems, we might not be able to afford a sparse factorization method. Or we might not have an explicit representation for the matrix! In this case, it is often natural to turn to iterative methods.

Iterative methods turn everything we’ve said about direct factorization on its head. It is not insane to go to a textbook and implement most of the standard sparse iterations out there. They often involve only a few lines of code, and while the analysis is subtle, the implementation mostly isn’t (with a few exceptions). Even if you use a standard implementation, you may find that you want to use an application-specific preconditioner, a transformation to the problem that converges acceleration of the iterative method. The choice of preconditioners is still as much art as science for most problems, and depends a lot on details of the application. Moreover, the convergence of these methods is often all over the place; they sometimes need a certain amount of hand-holding to get to the right solution. Worse, a preconditioner that is great for convergence may often be lousy for parallelism.

The good news about all of this is that there is a lot of room for creativity and there is still a lot of application-specific “low hanging fruit” when it comes to preconditioning. The bad news, if you just want to solve the problem, is that these methods often require a certain amount of hand-holding, particularly for problems that are “difficult” in some way.

LA Software: the Wider World

Dense: LAPACK, ScaLAPACK, PLAPACK, PLASMA, MAGMA
Sparse direct: UMFPACK, TAUCS, SuperLU, MUMPS, Pardiso, SPOOLES, ...
Sparse iterative: too many!
- PETSc (Argonne, object-oriented C)
- Trilinos (Sandia, C++)

If you are going to solve dense linear systems of equations, you’ll probably use LAPACK or some successor. Maybe you’ll use it via some other interface, like the Armadillo package in C++ or the solver interfaces in Julia and R and SciPy and MATLAB; but it’s LAPACK under the hood.

If you are going to do sparse direct solves, there’s a wider playing field. But it’s still a not-too-big set, mostly written by a handful of experts in the field. I’ve hyperlinked a survey article by Xiaoye Li, updated as of June 2020 – Xiaoye is the author of the SuperLU sparse solver (and is also my academic big sister).

For sparse iterative methods, there are too many to easily enumerate, and they all seem to have a lot of knobs. If you’re going to do things with these types of methods, I recommend getting to know either PETSc or Trilinos, two of the big parallel numerical framework libraries out there. They include implementations of many of these methods and preconditioners, and if you do something clever in this space, you probably want to use one of these frameworks to be able to compare with lots of competitors, and also to make sure that other people can benefity from your code.

Sparse Iterative Software

Assuming you want to use (vs develop) sparse LA software...

Identify a good algorithm (GMRES? CG?)
Pick a good preconditioner
- Often helps to know the application
- ... and to know how the solvers work!
Play with parameters, preconditioner variants, etc...
Swear until you get acceptable convergence?
Repeat for the next variation on the problem

What strategies should you use as a user of sparse iterative solvers? You probably don’t want to develop a new basic iteration: use a standard one like GMRES or CG. On the other hand, you will want to spend some time picking a good preconditioner. There are lots of default preconditioners that you can try, but if your problem is hard, you may need to go beyond the defaults. Or, at the very least, you’ll need to play with the parameters. Swearing may be involved. And passive voice may be used. Because I, of course, have never wanted to swear at an iteration that refused to converge.

(One of my grad school flatmates used to have a verbal running commentary when he was watching his iterations converge, or fail to do so. It was like listening to a sports commentator. It helped that he was doing earthquake simulations, so his “aah, oh no!” could either mean a convergence failure or a successful simulation that showed a catastrophic collapse of some structure. Fun times.)

Stacking Solvers

(Typical) example from a bone modeling package:

Outer load stepping loop
Newton method corrector for each load step
Preconditioned CG for linear system
Multigrid preconditioner
Sparse direct solver for coarse-grid solve (UMFPACK)
LAPACK/BLAS under that

I’ve been talking about these different types of solvers so far as if any given problem only requires one or the other. Of course, the picture is not that simple.

Let me tell you about a nonlinear finite element code that I wrote some years ago as a consulting gig. This was a code that simulated a hysteretic material, and so it mattered how it was loaded. So the outer loop was something that incrementally applied ever greater loads, usually up to a failure event. At each load step, we ran a Newton iteration to solve the force balance equation; and at each Newton step, we had a linear system that we solved with an iterative method (preconditioned conjugate gradients). To get reasonable convergence, we used a geometric multigrid precondiitioner that I wrote. For the “top level” coarse subproblem, I used a sparse direct solver, UMFPACK. And UMFPACK gets a lot of its performance by organizing as much as it can around dense sub-problems treated with LAPACK and level 3 BLAS subroutines. So in one code, you have the whole list of solvers that we’ve talked about so far.

By the way, I wrote the first three levels of the solver stack – the load stepper, the Newton iteration, and the CG solver – in the Lua scripting language. The multigrid code was in C++, and UMFPACK in C.

Overview of the week

Up next: Stationary iterations
Thurs: Krylov and sparse direct

Reading: Templates book