Computational Science Training for Undergraduates in the Mathematical Sciences (CSUMS)

Linear systems of equations with many right hand sides in lattice Quantum Chromodynamics

Project page:   www.cs.wm.edu/~andreas/CSUMS/linsys.html

Topic description: The numerical solution of linear systems of equations of large, sparse matrices is central to many scientific and engineering applications. One of the most computationally demanding applications is lattice Quantum Chromodynamics (QCD) because not only it involves matrix sizes of more than 10 million but also requires the solution of several linear systems with the same matrix but different right hand sides. Because of size, methods that transform the matrix (e.g, Gaussian elimination) cannot be used. Iterative methods provide the only means for solving these problems.

QCD is the theory of the fundamental force known as strong interaction, which describes the interactions among quarks, one of the constituents of matter. Lattice QCD is the tool for non-perturbative numerical calculations of these interactions on a Euclidean space-time lattice (4-D).
The heart of the computations is the solution of the lattice-Dirac equation, which translates to linear system of equations. Traditionally these linear systems are solved by applying the Conjugate Gradient (CG) on each system, one by one. Unknowingly, such methods regenerate information within previously explored subspaces, thus wasting iterations in successive right hand sides.

Deflation is the process by which eigenvalues that are near zero and their eigenvectors are precomputed and then removed from the operator (simply by orthogonalizing all new CG vectors against them). This speeds up the convergence significantly.

Research opportunity: We have recently devised a technique that computes these eigenvectors during the CG process, saving all the expense of running a separate eigenvalue solver. The performance of this new technique is impressive but there are open questions as to how to apply it and tune it. Deriving theory for this technique is the focus my other eigenvalue CSUMS project . In this project we deal more with issues of tuning the linear systems. It is very easy to extend a general purpose CG code to work with our technique, but this needs to be done.

There are other ways to improve convergence for CG. One is called preconditioning. In a few words, an M = approximate_inverse(A) is computed inexpensively, and then the following system is solved.
MA x = Mb.
For lattice QCD, this approximate inverse can be much more effectively computed if one considers the deflated A operator. This idea has not been explored but is highly promising.

Prerequisities: Math 211 Linear Algebra, Matlab and C/C++ programming skills.

Suggested prerequisities: Math 413/414 or some numerical analysis background

Contact: Andreas Stathopoulos