Nightmare Matrices: When Sparse Solvers Fail on GPUs

I introduce a family of difficult sparse matrices (“Nightmare matrices”) \(A\) designed to frustrate several common sparse direct and iterative strategies for solving \(Ax = b\). The construction is symmetric positive semidefinite, and is positive definite when the underlying random factor has full column rank, as it did here. I then analyze the performance of Krylov methods for these matrices on GPUs and show how the expander structure leads to uncoalesced memory access in cusparse::csrmv, causing warp stalls and poor GPU utilization per iteration. Next I demonstrate how deflation can improve GPU utilization by adding work the hardware handles efficiently while still accelerating convergence. Finally, I outline how to compute a deflation subspace that uses the GPU effectively and support these claims with performance counter and profiler results.

Nightmare Sparse Matrices

These matrices are created in three phases. First, nonzeros are distributed randomly across a matrix \(B\) (each with a value of 1) to match a prescribed number of nonzeros per row. The identity matrix is then added to ensure a nonzero diagonal, and the result is symmetrized. Random regular graph models are expanders with high probability. My construction is slightly different, but has the same expander-like separator structure that frustrates fill-reducing orderings. In the second phase the nonzero values (currently all 1s) are replaced with numbers uniformly distributed in [-1,1], centering the eigenvalue distribution. Finally, I set \(A = B^T B\). This increases the number of nonzeros and squares the singular values of \(B\), often making the smallest modes much smaller relative to the largest. The result is always positive semidefinite and is positive definite exactly when \(B\) has full column rank. Because the nonzero values are sampled randomly, the eigenvalue distribution of the intermediate symmetric factor tends to center near zero and taper off; squaring maps those modes to a nonnegative, poorly conditioned spectrum. This is a very poor (but not the absolute worst) distribution for unpreconditioned Krylov solvers. The worst case would be chebyshev points, but it’s not trivial to match an arbitrary spectrum to a fixed sparsity pattern.

I will investigate unpreconditioned Krylov methods in more detail later, but just to illustrate the poor performance I used Conjugate Gradients on a nightmare matrix with m=100000 rows on an NVIDIA T4 GPU. It reduced the error very slowly; the run below stopped with a relative solution error around 5e-3 rather than demonstrating mathematical nonconvergence. I use relative error here instead of relative residual because the exact manufactured solution is known.

it: 175593 relative error = 0.005438385018617882

real    3m52.127s
user    3m44.893s
sys     0m8.932s

Since unpreconditioned Krylov methods fail, we might consider preconditioners or sparse direct solvers, but these approaches also fall short

Failure of sparse direct solvers on Nightmare matrices

Expander graphs ensure that any two points in the underlying graph are connected by a short path. The graph can be quite sparse, yet in many important ways it behaves like a dense matrix. In particular, sparse direct methods produce significant fill-in regardless of the fill-reducing ordering, as illustrated below with a banded matrix as reference

For the tested sizes and orderings, these figures show factor storage approaching dense scaling even though \(A\) is extremely sparse. For matrices with 100000 rows this is impossible on a GPU with limited memory. We therefore turn to iterative methods and preconditioning, but those options mostly fail as well.

Failure of Sparsifying Preconditioners

“Sparsifying” preconditioners here means approaches that do not use the full matrix—for example incomplete LU, algebraic multigrid, or additive Schwarz/Block Jacobi. For Nightmare matrices these approaches fail to approximate the underlying system in important ways. One could design expander matrices where they succeed; for instance, Block Jacobi could work if off-diagonal blocks have low rank, even when the underlying graph is an expander. Nightmare matrices are specifically crafted to avoid such properties. I omit numerical results here to keep the post manageable, though I may explore them later.

But in a nutshell here is what ends up happening in each case:

1. Incomplete LU
   • Drops so many nonzeros to meet the threshold that the approximation collapses
2. Algebraic multigrid
   • Coarse spaces become dense quickly, defeating the point.
3. Additive Schwarz/Block Jacobi
   • Discards too many off-diagonal components

These techniques may offer some benefit, but not enough to outweigh the additional memory and computational burden they impose

Solving Nightmare Systems on GPU

Given the failure of sparse direct solvers, we turn to iterative methods. Most preconditioners offer little help, leaving unpreconditioned Krylov methods such as Conjugate Gradients. Here we encounter another awful property of Nightmare matrices: evaluating \(Ax\) for any \(x\) resembles pure random access, the worst possible memory-access pattern on the GPU.

Plain Conjugate Gradient

Each step of Conjugate Gradient involves a single evaluation of \(Ax\), and this dominates the runtime. For a Nightmare matrix that work makes poor use of the GPU. I profiled a run with m=100000 and ~60 nonzeros per row using cupy which mostly calls into cusparse in this case. The table below from an nsys profile confirms this, with 70% of the time spent in csrmv which is the cusparse sparse matrix-vector product for CSR matrices

The problem with csrmv is that the matrix comes from an expander graph, so access is nearly random. This leads to significant uncoalesced memory traffic and warp stalls on the GPU, as shown by the performance counters below

The profile shows poor DRAM bandwidth (21% utilization) and few warps in flight because uncoalesced reads cause stalls. In the next section I show how to recover GPU utilization by doing more work per iteration, but choosing work the GPU executes efficiently.

Deflated Conjugate Gradient

Deflated methods seek an orthonormal matrix \(V\) that projects out components slowing Krylov convergence. The challenge is computing such a \(V\) and integrating it with a Krylov method, each of which handles \(V\) slightly differently. The scheme below follows Kahl and Rittich.

For ( A x = b ), let ( V \in \mathbb{R}^{n \times k} ) have orthonormal columns and define

\[ M = V^\top A V. \]

The deflated system is

\[ (I - A V M^{-1} V^\top)A x_h = (I - A V M^{-1} V^\top)b. \]

After solving for \( x_h \), the reconstruction is

\[ x = V M^{-1} V^\top b + \bigl( x_h - V M^{-1} V^\top A x_h \bigr). \]

I choose \(V\) to be the matrix of eigenvectors associated with the smallest-magnitude eigenvalues of \(A\) and postpone discussion of how to compute it until the next section. This adds work per iteration on top of evaluating \(Av\), but the extra work is productive for the following reasons:

1. It involves level-2 BLAS ( \(Vx\), \(V^Tx \) with dense \(V\) ), which coalesces well and achieves higher GPU utilization
2. \( V \) projects out the extremely small eigenvalue components which caused our observed CG to stagnate

and indeed we see point (1) plainly in the trace, with aggregates below. Almost 80% of the time goes to *gemv-style level-2 BLAS, and the share spent on csrmv drops from about 70% to roughly 9%. The csrmv call is unchanged; each iteration now includes work that uses the GPU efficiently, so far fewer iterations are needed to reach the same error (quantified later, after the GPU efficiency discussion).

Our utilization metrics improve markedly. DRAM bandwidth nearly triples, and many more warps are active because memory access patterns become largely coalesced.

And just for ease of comparison I include both runs - with and without deflation - below

and finally the payoff that we actually care about: achieving much lower error in significantly less solve time once the deflation space is available:

The deflation basis is not free. For \(m=100000\), \(k=1024\), and FP64 storage, \(V\) alone occupies about 0.8 GB, and the table excludes the cost of computing it. This setup makes the most sense when it can be reused across many right-hand sides.

Computing the deflation subspace

The results from deflated conjugate gradient are impressive but of course require computing a suitable \( V \). This takes us back to the beginning in a way because now we need an eigensolver capable of producing the smallest-magnitude eigenvalues for \( A \), which is a nightmare matrix and many straightforward techniques are unattractive. Standard Krylov approaches include Arnoldi or Lanczos iterations, but their single-vector recurrences introduce a dependency between matrix-vector products and require maintaining orthogonality with previous iterates. I have already established that evaluating \(Ax\) for single vectors uses the GPU poorly. The solution is to attempt a block method that jointly advances multiple eigenvectors simultaneously. This brings us back to coalesced accesses because evaluating \(AX\) for a matrix \(X\) allows each loaded nonzero of \( A \) to be reused across several columns of \(X\). The exact coalescing behavior depends on the sparse-kernel layout and how the dense block is stored.

To accomplish this I make use of the fact that symmetric matrices can be diagonalized with orthonormal eigenvectors \(W \)

\[ A = W \Lambda W^T \]

and this decomposition shows us that all polynomials \( P \) satisfy

\[ P(A) = W P(\Lambda) W^T \]

This means that if we choose \( P \) to approximate \( P(x) \approx \frac{1}{x} \) then we can magnify the small eigenvalues. Now this is basically restating the underlying principle of a Krylov solver, but we can statically compute a polynomial for this as well. This is the principle underlying Chebyshev iteration.

For an SPD \(A\) with \(\operatorname{spec}(A)\subset[\lambda_{\min},\lambda_{\max}]\), set \[ \mu = \tfrac{\lambda_{\max}+\lambda_{\min}}{2}, \qquad \nu = \tfrac{\lambda_{\max}-\lambda_{\min}}{2}. \]

Initialize \(x^{(0)}\) and \(r^{(0)} = b - A x^{(0)}\). Let \(\alpha_0 = \tfrac{1}{\mu},\; \beta_0 = 0\) and form \[ x^{(1)} = x^{(0)} + \alpha_0 r^{(0)}. \]

For \(k \ge 1\), \[ \beta_k = \left(\frac{\nu \,\alpha_{k-1}}{2}\right)^2, \qquad \alpha_k = \frac{1}{\mu - \dfrac{\beta_k}{\alpha_{k-1}}}, \] \[ x^{(k+1)} = x^{(k)} + \alpha_k r^{(k)} + \beta_k \big(x^{(k)} - x^{(k-1)}\big), \qquad r^{(k)} = b - A x^{(k)}. \]

The advantage of Chebyshev iteration in this context compared to a separately adapted Krylov solve for every column is that it is simple to advance multiple vectors simultaneously because the polynomial coefficients they get are all exactly the same (as opposed to say Conjugate Gradient, which computes the coefficients through dot products of previous iterates, which will result in different coefficients for each vector).

Here is the pseudocode for computing \( V \) with all of the above observations in mind:

# Chebyshev–Rayleigh–Ritz for k smallest eigenpairs
# Assume a routine:
#   Y = ChebyshevSolve(A, RHS, inner, λ_min, λ_max)
# which does 'inner' Chebyshev iterations to approximately solve A·Y = RHS.

Inputs:
  A            # symmetric SPD n×n operator (matvec/matmat supported)
  k            # number of eigenpairs
  outer        # outer RR iterations
  inner        # inner Chebyshev iterations
  λ_min, λ_max # spectral bounds

Procedure:
  V ← RandomOrthonormal(n, k)      # e.g., QR of random matrix
  Θ ← undefined

  for t = 1..outer:
      # 1) Block Chebyshev solve: approximate A^{-1}·V
      W ← ChebyshevSolve(A, V, inner, λ_min, λ_max)   # n×k

      # 2) Orthonormalize subspace
      Q, _ ← qr(W)                                    # n×k with QᵀQ = I

      # 3) Rayleigh–Ritz on span(Q)
      T ← Qᵀ A Q                                      # k×k
      S, diag(Θ) ← eig(T)                             # Θ ascending
      V ← Q S                                         # n×k (orthonormal)

      # (optional) stop if max_j ||A v_j − Θ_j v_j|| ≤ tol
  end for

Outputs:
  V  # Ritz vectors (≈ k smallest eigenvectors)
  Θ  # Ritz values (≈ k smallest eigenvalues)

I also provide code in the appendices.

From the trace we can see significant time in cusparse::spmm (note the asterisk: I renamed this from the symbol in cusparse which had a lot of template magic in it) and subsequent calls to DGEMM which is level-3 blas, absolutely the best place to be on a GPU.

Here we can see extremely high utilization of the GPU with significant number of warps active per cycle, fully loaded GR and SM, high memory bandwidth.

For the deflated case with m=100000 I was able to compute 512 highly accurate eigenvectors (relative residuals around 1e-12); accuracy then decayed for the remaining requested vectors. This was the practical limit of my run rather than an inherent cutoff of the block method.

Appendices

Machine

I used a GCP instance for these measurements:

Code

I have started collecting many related functions into a github repository that has examples and other useful functions not showcased here as well.

References

1. Friedman, Joel. “A Proof of Alon’s Second Eigenvalue Conjecture and Related Problems.” arXiv:cs/0405020. Preprint, arXiv, May 5, 2004. https://doi.org/10.48550/arXiv.cs/0405020.
2. Kahl, K., and H. Rittich. “The Deflated Conjugate Gradient Method: Convergence, Perturbation and Accuracy.” Linear Algebra and Its Applications 515 (February 2017): 111–29. https://doi.org/10.1016/j.laa.2016.10.027.
3. Zhaojun, Bai, Demmel James, Dongarra Jack, Ruhe Axel, and van der Vorst Henk, eds. Templates for the Solution of Algebraic Eigenvalue Problems. Software, Environments and Tools. Society for Industrial and Applied Mathematics, 2000. https://doi.org/10.1137/1.9780898719581.