Solving Sparse Linear Systems Faster than Matrix Multiplication

At a high level, our algorithm follows the block Krylov space method (see e.g., Chapter 6.12 of Saad¹⁶). This method is a multi-vector extension of the CG/Lanczos method, which in the single-vector setting is known to be problematic under round-off errors both in theory¹² and in practice.¹⁶ Our algorithm starts with a set of $s$ initial vectors, $B\in {\Re }^{n\times s}$, and forms a column space by multiplying these vectors by $A$ repeatedly, $m$ times. Formally, the block Krylov space matrix is

The core idea of Krylov space methods is to efficiently orthogonalize this column space. For this space to be spanning, block Krylov space methods typically choose $s$ and $m$ so that $sm=n$.

The conjugate gradient algorithm can be viewed as an efficient implementation of the case $s=1$, $m=n$, with $B$ set to $b$, the RHS of the input linear system. The block case with larger values of $s$ was studied by Eberly, Giesbrecht, Giorgi, Storjohann, and Villard⁵ over finite fields, and they gave an $O\left(n^{2.28}\right)$  time^d algorithm for computing the inverse of an $O\left(n\right)$-sparse matrix over a finite field.

Our algorithm also leverages the top-level insight of the Eberly et al. results: the Gram matrix of the Krylov space matrix, is a block Hankel matrix. Solving linear systems in this Gram matrix, ${\left(AK\right)}^T\left(AK\right)$, lead to solvers for linear systems in $A$ because

so as long as $A$ and $K$ are both invertible, composing this on the left by $K$ and on the right by $A\top K^{\top }$ gives

Eberly et al. viewed the Gram matrix as an $m$-by-$m$ matrix containing $s$-by-$s$ sized blocks, and critically leveraged the fact that the blocks along each anti-diagonal are identical:

Formally, the $s$-by-$s$ inner product matrix formed from $A^{i}B$ and $A^{j}B$ is $B^T{A}^{i+j}B$, and depends only on $i+j$. So instead of $m^2$ blocks each of size $s\times s$, we are able to represent an $n$-by-$n$ matrix with only about $m$ blocks.

Operations involving these $m$ blocks of the Hankel matrix can be handled using $\tilde{O}\left(m\right)$ block operations. This is perhaps easiest seen for computing matrix-vector products using $K$. If we use $\left\{i\right\}$ to denote the $i$th block of the Hankel matrix, and define

for a sequence of matrices $M$, we get that the $i^{\mathrm{th}}$ block of the product $Hx$ can be written in block-form as

Note this is precisely the convolution of (a sub-interval) of $M$ and $x$, with shifts indicated by $i$. Therefore, in matrix-vector multiplication (the “forward” direction), a speedup by a factor of about $m$ is possible with fast convolution algorithms. The performance gains of the Eberly et al. algorithms⁵ can be viewed as being of a similar nature, albeit in the more difficult direction of solving linear systems. Specifically, they utilize algorithms for the Padé problem of computing a polynomial from the result of its convolution.¹ Over finite fields, or under exact arithmetic, such algorithms for matrix Padé problems take $O(m log m)$ block operations,¹ for a total of $\tilde{O}\left(s^{\omega }m\right)$ operations.

The overall time complexity is based on two opposing goals:

Under these assumptions, and the requirement of $n\approx ms$, the total cost becomes about $O(n·nnz+m·s^{\omega })$, which is at most $O(n·nnz)$ as long as $m>n^{\frac{\omega -2}{\omega -1}}$. However, this runtime complexity is over finite fields, where numerical stability is not an issue. Over the reals, under round-off errors, one must contend with numerical errors without blowing up the bit complexity. This is a formidable challenge; indeed, as mentioned earlier, with exact arithmetic, the CG method takes time $O(n·nnz)$, but this is misleading since the computation is effective only when the word sizes are increased by a factor of $n$ (to about $n log \kappa$ words), which leads to an overall complexity of $O(n^2·nnz· log \kappa )$.

Our algorithm can be viewed as the numerical generalization of the algorithms from.⁵ We work with real numbers of bounded precision, instead of entries over a finite field. The core of our approach can be summarized as follows.

Doing so requires separately developing tools for two topics that have been extensively studied in mathematics:

Before we describe the difficulties and new tools needed, we first provide some intuition on why a factor $m$ increase in word lengths may be the right answer by upper-bounding the magnitudes of entries in an $m$-step Krylov space. By rescaling, we may assume that the minimum singular value of $A$ is at least $1/\kappa$, and the maximum entry in $A$ is at most $\kappa$. The maximum magnitude of (entries of) $A^{m}b$ is bounded by the maximum magnitude of $A$ to the power of $m$, times a factor corresponding to the number of summands in the matrix product:

Here the last inequality is via the assumption of $\kappa \ge n$. So by forming $K$ via horizontally concatenating $A^{i}g$ for sparse Gaussian vectors $g$, we have with high probability that the maximum magnitude of an entry of $K$, and in turn $AK$, is at most ${\kappa }^{O\left(m\right)}$. In other words, $O(m log \kappa )$ words in front of the decimal point is sufficient with high probability.

Should such a bound of $O(m log \kappa )$ hold for all numbers that arise in the algorithm, including the matrix inversion steps, and the matrix $B$ is sparse with $O\left(n\right)$ entries, the cost of computing the block-Krylov matrices becomes $O(m log \kappa ·ms·nnz)$, while the cost of the matrix inversion portion encounters an overhead of $O(m log \kappa )$, for a total of $\tilde{O}(m^2{s}^{\omega } log \kappa )$. In the sparse case of $nnz=O\left(n\right)$, and $n\approx ms$, this becomes:

Due to the gap between $n^2$ and $n^{\omega }$, setting $m$ appropriately gives improvement over $n^{\omega }$ when $log \kappa =n^{o\left(1\right)}$.

However, the magnitude of an entry in the inverse depends on the smallest magnitude, or in the matrix case, its minimum singular value. Bounding and propagating the min singular value, which intuitively corresponds to how close a matrix is to being degenerate, represents our main challenge. In exact/finite fields settings, non-degeneracies are certified via the Schwartz-Zippel lemma about polynomial roots. The numerical analog of this is more difficult — the Krylov space matrix $K$ is asymmetric, even for a symmetric matrix $A$. It is much easier for an asymmetric matrix with correlated entries to be close to singular.

Consider for example a two-banded, two-block matrix with all diagonal entries set to the same random variable $\alpha$ (see Figure 1):

In the exact case, this matrix is full rank unless $\alpha =0$, even over finite fields. On the other hand, its minimum singular value is close to 0 for all values of $\alpha$. To see this, it’s useful to first make the following observation about the min singular value of one of the blocks.

Then, in the top-left block, as long as $\left|\alpha \right|>3/2$, the top left block has minimum singular value at most ${(2/3)}^{n-1}$. On the other hand, rescaling the bottom-right block by $1/\alpha$ to get 1s on the diagonal gives $2/\alpha$ on the off-diagonal. So as long as $\left|\alpha \right|<3/2$, this value is at least $4/3$, which in turn implies a minimum singular value of at most ${(3/4)}^{n-1}$ in the bottom right block. This means no matter what value $\alpha$ is set to, this matrix will always have a singular value that’s exponentially close to 0. Furthermore, the Gram matrix of this matrix also gives such a counter example to symmetric matrices with (non-linearly) correlated entries. Previous works on analyzing condition numbers of asymmetric matrices also encounter similar difficulties; a more detailed discussion of it can be found in Section 7 of Sankar et al.¹⁸

In order to bound the bit complexity of all intermediate steps of the block Krylov algorithm by $\tilde{O}(m· log \kappa )$, we devise a more numerically stable algorithm for solving block Hankel matrices, as well as provide a new perturbation scheme to quickly generate a well-conditioned block Krylov space. Central to both of our key components is the close connection between condition number and bit complexity bounds.

First, we give a more numerically stable solver for block Hankel/Toeplitz matrices. Fast solvers for Hankel (and closely related Toeplitz) matrices have been extensively studied in numerical analysis, with several recent developments on more stable algorithms.²⁴ However, the notion of numerical stability studied in these algorithms is the variant where the number of bits of precision is fixed. Our attempts at converting these into asymptotic bounds yielded dependencies quadratic in the number of digits in the condition number, which in our setting translates to a prohibitive cost of $\tilde{O}\left(m^2\right)$ (i.e., the overall cost would be higher than $n^{\omega }$).

Instead, we combine developments in recursive block Gaussian elimination⁴ with the low displacement rank representation of Hankel/Toeplitz matrices.⁷ Such representations allow us to implicitly express both the Hankel matrix and its inverse by displaced versions of low-rank matrices. This means the intermediate size of instances arising from recursion is $O\left(s\right)$ times the dimension, for a total size of $O(n log n)$, giving a total of $\tilde{O}\left(n{s}^{\omega -1}\right)$ arithmetic operations involving words of size $\tilde{O}\left(m\right)$. We provide a rigorous analysis of the accumulation of round-off errors similar to the analysis of recursive matrix multiplication based matrix inversion from.⁴

Motivated by this close connection with the condition number of Hankel matrices, we then try to initialize with Krylov spaces of low condition number. Here we show that a sufficiently small perturbation suffices for producing a well conditioned overall matrix. In fact, the first step of our proof, showing that a small sparse random perturbation to $A$ guarantees good separations between its eigenvalues, is a direct combination of bounds on eigenvalue separation of random Gaussians¹³ as well as min eigenvalue of random sparse matrices.¹¹ This separation then ensures that the powers of $A$, $A^1,A^2,...A^m$, are sufficiently distinguishable from each other. Such considerations also come up in the smoothed analysis of numerical algorithms.¹⁸

The randomness of the Krylov matrix induced by the initial set of random vectors $B$ is more difficult to analyze: each column of $B$ affects $m$ columns of the overall Krylov space matrix. In contrast, all existing analyses of lower bounds of singular values of possibly asymmetric random matrices^18,21 rely on the randomness in the columns of matrices being independent. The dependence between columns necessitates analyzing singular values of random linear combinations of matrices, which we handle by adapting $ϵ$-net based proofs of anti-concentration bounds. Here we encounter an additional challenge in bounding the minimum singular value of the block Krylov matrix. We resolve this issue algorithmically: instead of picking a Krylov space that spans the entire ${\Re }^n$, we stop short by picking $ms=n-\tilde{O}\left(m\right)$ This set of extra columns significantly simplifies the proof of singular value lower bounds. This is similar in spirit to the analysis of the minimum singular value of a random matrix, which is easier for a non-square matrix.¹⁵ In the algorithm, the remaining columns are treated as a separate block that we handle via a Schur complement at the very end of the algorithm. Since this block is small, so is its overhead on the running time.

Our algorithm has close connections with multiple lines of research on efficient solvers for sparse linear systems. The topic of efficiently solving linear systems has been extensively studied in computer science, applied mathematics and engineering. For example, in the Society of Industrial and Applied Mathematics News’ ‘top 10 algorithms of the 20th century’, three of them (Krylov space methods, matrix decompositions, and QR factorizations) are directly related to linear systems solvers.³

At a high level, our algorithm is a hybrid linear systems solver. It combines iterative methods, namely block Krylov space methods, with direct methods that factorize the resulting Gram matrix of the Krylov space. Hybrid methods have their origins in the incomplete Cholesky method for speeding up elimination/factorization based direct solvers. A main goal of these methods is to reduce the $\Omega \left(n^2\right)$ space needed to represent matrix factorizations/inverses. This high space requirement is often even more problematic than time requirements when handling large sparse matrices. Such reductions can occur in two ways: either by directly dropping entries from the (intermediate) matrices, or by providing more succinct representations of these matrices using additional structure.

The main structure of our algorithm is based on the latter line of work on solvers for structured matrices. Such systems arise from physical processes where the interactions between objects have invariances (e.g., either by time or space differences). Examples of such structure include circulant matrices, Hankel/Toeplitz matrices and distances from $n$-body simulations.⁷ Many such algorithms require exact preservation of the structure in intermediate steps. As a result, many of these works develop algorithms over finite fields.

More recently, there has been work on developing numerically stable variants of these algorithms for structured matrices, or more generally, matrices that are numerically close to being structured.²⁴ However, these results only explicitly discussed in the entry-wise Hankel/Toeplitz case (which corresponds to block size $s=1$). Furthermore, because they rely on domain-decomposition techniques similar to fast multiple methods, they produce one bit of precision for each outer iteration loop. As the Krylov space matrix has condition number $exp(\Omega (m\left)\right)$, such methods would lead to another factor of $m$ in the solve cost if invoked directly.

Instead, our techniques for handling and bounding numerical errors are more closely related to recent developments in provably efficient sparse Cholesky factorizations.⁹ These methods generate efficient preconditioners using only the condition that intermediate steps of Gaussian elimination, known as Schur complements, have small representations. They avoid the explicit generation of the dense representations of Schur complements by treating them as operators, and apply randomized tools to directly sample/sketch the final succinct representations, which have much smaller size and algorithmic cost.

On the other hand, previous works on sparse Cholesky factorizations required the input matrix to be decomposable into a sum of simple elements, often through additional combinatorial structure of the matrices. In particular, this line of work on combinatorial preconditioning was initiated through a focus on graph Laplacians, which are built from 2-by-2 matrix blocks corresponding to edges of undirected graphs.¹⁹ Since then, there have been substantial generalizations to the structures amenable to such approaches, notably to finite element matrices, and directed graphs/irreversible Markov chains. However, recent works have also shown that many classes of structures involving more than two variables are complete for general linear systems.²⁵ Nonetheless, the prevalence of approximation errors in such algorithms led to the development of new ways to bound numerical round-off errors in algorithms, which will be critical to our elimination routine for block-Hankel matrices.

Key to recent developments in combinatorial preconditioning is matrix concentration.²² Such bounds provide guarantees for (relative) eigenvalues of random sums of matrices. For generating preconditioners, such randomness arises from whether each element is kept or not, and a small condition number (which in turn implies a small number of outer iterations using the preconditioner) corresponds to a small deviation between the original and sampled matrices. In contrast, we introduce randomness in order to obtain block Krylov spaces whose minimum eigenvalue is large. As a result, the matrix tool we need is anti-concentration, which somewhat surprisingly is far less studied. Previous works on it are related to similar problems of numerical precision^18,21 and address situations where the entries in the resulting matrix are independent. Our bound on the min singular value of the random Krylov space also yields a crude bound for a sum of rectangular random matrices.

Nie¹⁴ gave a more general and tighter version of matrix anti-concentration that also works for square matrices, answering an open question we posed. For an $m$-step Krylov space instantiated using $\frac{n}{m}$ vectors, Nie’s bound reduces the middle term in our analysis from $O\left(n^2{m}^3\right)$ to $O\left(n^{2}m\right)$, thus leading to a running time of $nnz·n·m+n^{\omega }·m^{2-\omega }$, which matches the bound for finite fields. Moreovver, it does so without the padding step at the end. We elect to keep for this article our epsilon-net based analyses, and the padded algorithm required for it, both for completeness, and due to it being a more elementary approach toward the problem with a simpler proof.

Faster matrix multiplication is an active area of work with recent progresses. Due to the dependence on fast multiplication in our algorithm, such improvements also lead to improvements in the running time of solving sparse systems as well.

Our main result for solving sparse linear systems has also been extended to solving sparse regression, with faster than matrix multiplication bounds for sufficiently sparse matrices.⁶ The complexity of sparse linear programming remains an interesting open problem.