TOWARDS UNIVERSAL CONVERGENCE OF BACKWARD ERROR IN LINEAR SYSTEM SOLVERS∗
arXiv:2604.16075v1 [math.NA] 17 Apr 2026
MICHAL DEREZIŃSKI† , YUJI NAKATSUKASA‡ , AND ELIZAVETA REBROVA§ Abstract. The quest for an algorithm that solves an n × n linear system in O(n2 ) time complexity, or O(n2 poly(1/ϵ)) when solving up to ϵ relative error, is a long-standing open problem in numerical linear algebra and theoretical computer science. There are two predominant paradigms for measuring relative error: forward error (i.e., distance from the output to the optimum solution) and backward error (i.e., distance to the nearest problem solved by the output). In most prior studies, convergence of iterative linear system solvers is measured via various notions of forward error, and as a result, depends heavily on the conditioning of the input. Yet, the numerical analysis literature has long advocated for backward error as the more practically relevant notion of approximation. In this work, we show that — surprisingly — the classical and simple Richardson iteration incurs at most 1/k (relative) backward error after k iterations on any positive semidefinite (PSD) linear system, irrespective of its condition number. This universal convergence rate implies an O(n2 /ϵ) complexity algorithm for solving a PSD linear system to ϵ backward error, and we establish similar or better complexity when using a variety of Krylov solvers beyond Richardson. Then, by directly minimizing 2 backward error over a Krylov subspace, we attain an even faster O(1/k √ ) universal rate, and we turn this into an efficient algorithm, MINBERR, with complexity O(n2 / ϵ). We extend this approach via normal equations to solving general linear systems, for which we empirically observe O(1/k) convergence. We report strong numerical performance of our algorithms on benchmark problems. Key words. Linear systems, Backward error, Krylov methods, Richardson iteration MSC codes. 65F10, 65F35, 65Y20
1. Introduction. Solving a system of linear equations Ax = b given an n × n matrix A and vector b ∈ Rn is perhaps the most important problem in computational mathematics. Direct methods, such as those based on Gaussian elimination with partial pivoting or QR factorization, remain the de facto standard for solving linear systems when n is moderate, say n < 104 in today’s desktop computers and n < 106 with supercomputers. However, the arithmetic complexity of these methods scales as O(n3 ), making them infeasible when n is much larger. An alternative approach is to use fast Strassen-type matrix multiplication algorithms, which improve the asymptotic complexity of solving linear systems to O(nω ), where currently ω ≈ 2.371 [2], but have remained largely theoretical due to extremely large constant factors. All of the aforementioned approaches tend to produce highly accurate solutions (modulo numerical precision issues), but at a relatively high computational cost. Yet, in many applications of linear system solvers we are willing to trade off accuracy for computation by relying on approximate solutions that are computed within some ϵ error relative to the optimum. This has motivated the long-standing question: Can we compute an ϵ relative error solution of an n × n linear system in O( n2 poly(1/ϵ) ) time? In response to this question, many iterative methods have been developed, from the classical iteration scheme proposed by Richardson in 1911 [40] to the modern ∗ Funding: MD was supported in part by NSF CAREER Grant CCF-233865 and a Google ML and Systems Junior Faculty Award. This work was done in part while MD was visiting the Simons Institute for the Theory of Computing. YN was supported by the EPSRC grant EP/Y030990/1. ER was partially supported by NSF DMS-2309685. † University of Michigan, Ann Arbor, MI ([email protected]). ‡ University of Oxford, Oxford, UK ([email protected]). § Princeton University, Princeton, NJ ([email protected]).
1
2
M. DEREZIŃSKI, Y. NAKATSUKASA, AND E. REBROVA
Krylov subspace methods [22]. These algorithms can quickly produce a rough approximation to a linear system by making incremental updates towards the solution. The complexity of a single iteration typically scales with the cost of applying A to a vector (denoted here by TA ), which is no more than O(n2 ) arithmetic operations for dense matrices. However, the effectiveness of iterative methods is highly dependent on the spectral properties of A. In particular, existing convergence guarantees [41] suggest that for a matrix with condition number κ(A) = ∥A∥2 ∥A−1 ∥2 , as many as O(κ(A) log(1/ϵ)) iterations may be needed to attain a relative p error approximation using popular algorithms such as LSQR [37], or at least O( κ(A) log(1/ϵ)) using Conjugate Gradient (CG, [24]) if A is positive semidefinite (PSD). While these guarantees can be finetuned given additional information about the spectrum of A [5], they suggest that the answer to our question depends on the conditioning of the input matrix. Is there no hope then of attaining universal convergence guarantees for iterative methods, i.e., relative error bounds that converge to zero at a rate that is independent of the input? Naturally, a proper answer to this question requires specifying what “relative error” refers to. The majority of existing results measure the convergence of iterative methods through the notion of a relative forward error [41]. Given an output vector x ∈ Rn , its (relative) forward error with respect to a norm ∥ · ∥, matrix A, and vector b is defined as: (forward error)
min x∗
∥x − x∗ ∥/∥x∗ ∥
s.t.
Ax∗ = b.
Here, the norm ∥ · ∥ most commonly refers to the 2-norm ∥ · ∥2 , or if A is PSD, the A-norm ∥ · ∥A . Equally common is what is known as the relative residual1 ∥Ax − b∥2 /∥b∥2 ; this can be viewed as the forward error with respect to the A⊤ Anorm ∥ · ∥A⊤ A . Thanks to the normalization, we can assert to have obtained a useful approximation as long as the error is sufficiently less than 1. Unfortunately, even attaining forward error 1/2 may take an enormous number of iterations when dealing with an ill-conditioned system, and existing lower bounds show that dependence on condition number is in general unavoidable among all algorithms that access A through matrix-vector products [12]. Despite the prevalence of forward error in the convergence analysis of iterative methods, much of the numerical analysis literature has advocated for the backward error as a more fundamental notion of approximation quality in numerical linear algebra [25]. Informally speaking, small backward error means we have the exact solution for a slightly perturbed problem. This accounts for the fact that input data to any given problem has likely already undergone some distortion before reaching the algorithm, and thus it is sufficient to solve a perturbed problem instance. Given an output vector x ∈ Rn , its (relative) backward error with respect to a norm ∥ · ∥, matrix A, and vector b is defined as: (backward error)
min
∥Ã − A∥/∥A∥
s.t. Ãx = b.
Ã
Here, the norm ∥ · ∥ typically refers to the 2-norm (which is what we will use throughout), or possibly the entrywise infinity norm, and one may also allow distortion in vector b (which we omit for simplicity). 1 The residual is often used to measure convergence in practice, because ∥x − x ∥ cannot be ∗ 2 computed without knowing x∗ . Backward error, by contrast, can be computed easily for linear systems (see [25] and Lemma 3.2) and least-squares problems [15].
TOWARDS UNIVERSAL CONVERGENCE IN LINEAR SYSTEM SOLVERS
3
Backward error analysis of numerical methods dates back over 60 years to classical works of Wilkinson [47, 48], and even before that to von Neumann and Goldstine [46]. Since then, numerical linear algebra textbooks [44, 25, 20] have argued that backward error is the most practically relevant metric of approximation quality, and popular linear system solvers such as LSQR [37] have been using it as the default stopping criterion [3, 38]. Currently, with the rise of mixed-precision computations and compressed data storage on GPU-based hardware architectures [26], as well as recent applications in training machine learning models [6, 7], backward error has an even more central role to play in the development of numerical algorithms. Fortunately, as an approximation metric for linear systems, backward error is provably less restrictive than forward error as it is insensitive to small perturbations of the input. This would suggest that it should enjoy better (ideally, universal) convergence guarantees for iterative linear system solvers. Yet, perhaps surprisingly, no such guarantees are currently known in the literature. Main results for PSD linear systems. In this work, we initiate the study of convergence in backward error for iterative linear system solvers, showing that when viewed from this perspective, our motivating question has a simple resolution for the class of PSD linear systems. We start with Richardson iteration, and prove that it attains a simple universal convergence guarantee (Theorem 3.3): After k steps of Richardson iteration on a PSD linear system, the backward error is at most 1/k. Note that, again due to normalization, we can assert to have obtained a meaningful approximation as long as backward error is sufficiently less than 1. In particular, to attain an ϵ ∈ (0, 1) backward error Richardson requires ⌈1/ϵ⌉ iterations and has O((TA +n)/ϵ) = O(n2 /ϵ) complexity. In comparison, the standard convergence analysis yields O(n2 κ(A) log(1/ϵ)) complexity which applies to both forward and backward error. Naturally, if the condition number of A is small, then the standard analysis gives an exponentially faster convergence rate, but if κ(A) is large, then our universal rate shows that Richardson makes substantial progress in backward error long before it makes any progress in forward error. Our proof of Theorem 3.3 is closely connected with the common intuition that Richardson converges faster along the top eigendirections, a phenomenon which is captured much better by the backward error than it is by the forward error. Can Richardson’s 1/k convergence rate for psd linear systems be improved by utilizing modern Krylov solvers? Interestingly, while methods such as MINRES and CG are not designed out-of-the-box to minimize backward error, we can impose a deliberate perturbation in the form of regularization, replacing matrix A with A + δI for an appropriately chosen δ, to ensure that after k steps the backward error of CG or MINRES is bounded by O(log2 (k)/k 2 ). However, this is still not the best universal rate attainable by a Krylov method. We show this by developing and analyzing an algorithm, MINBERR, that directly minimizes backward error over the Krylov subspace, and attains an O(1/k 2 ) convergence rate (Theorem 4.1). We provide an efficient implementation of MINBERR by exploiting the structure and symmetry of the Lanczos decomposition for a PSD matrix A, attaining an algorithm that takes O((TA + n)k) time to perform k iterations. This allows us to further improve on the complexity guarantee attained by Richardson iteration (Corollary 4.9): We can compute an ϵ backward error solution of an n × n PSD linear system √ in O( n2 / ϵ ) time.
4
M. DEREZIŃSKI, Y. NAKATSUKASA, AND E. REBROVA
Further results for general linear systems. Can universal convergence guarantees be obtained for general linear systems? One might expect that running the aforementioned PSD solvers on the normal equations of a general system should yield some (possibly weaker) universal rates in terms of backward error. Unfortunately, we show that Richardson on the normal equations exhibits an initial blow-up of the backward error by up to a condition number factor, and thus cannot attain any universal convergence rate (Theorem 3.4). We also observe a similar phenomenon empirically for Krylov solvers such as LSQR (Figure 3). We investigate this further by developing a variant of MINBERR which minimizes backward error over the Krylov subspace defined by the normal equations (MINBERR-NE). Remarkably, this method empirically exhibits strong convergence and avoids the blow-up of Richardson. However, we are able to construct a family of hard problem instances for which the MINBERR-NE backward error is at least log(κ(A))/k after k iterations, suggesting that no deterministic Krylov method can attain truly universal convergence for general systems. We propose to address this by introducing a random perturbation, replacing A with A + δG for a carefully chosen δ and a Gaussian matrix G, similarly to what is done in smoothed analysis [43]. This yields a condition number proportional to n while minimally affecting backward error, thereby making MINBERR-NE independent of condition number even on our hard problem instances. 2. Background and Related Work. Here, we provide additional background on backward error and convergence analysis of linear system solvers, while also highlighting an important connection to continuous optimization literature. Backward error analysis. Trefethen and Bau in their Numerical Linear Algebra textbook [44] describe backward error analysis as “a fundamental idea linking conditioning and stability, whose power has been proved in innumerable applications since the 1950s”. This idea was primarily initiated by Wilkinson, whose books on “Rounding Errors in Algebraic Processes” [49, 1963] and “The Algebraic Eigenvalue Problem” [50, 1965] have achieved the status of classics, according to Higham [25] who provides a modern treatment of the subject. Developing backward error analysis for iterative linear system solvers and Krylov subspace methods has been the subject of extensive work, including by Paige [36] and Greenbaum [21] on the stability of CG/Lanczos, and it remains an active area of research to this day [33]. These works focus on how rounding errors in finite precision arithmetic affect the backward error relative to an exact execution of the algorithm (which is different from the true backward error, since the algorithm returns an approximate solution even in exact arithmetic). Paige and Saunders [37] suggested backward error tolerance as a natural stopping criterion for LSQR,2 and later, Fong and Saunders [18] cited strong backward error performance as a highly desirable advantage of their LSMR algorithm. Kasenally and Simoncini proposed GMBACK [27, 28], which minimizes backward error over a Krylov subspace, similarly to MINBERR and MINBERR-NE, except with a much higher per-iteration cost since it does not exploit symmetry via the Lanczos decomposition. Backward error has also been recently used in machine learning to interpret model training via stochastic gradient descent as solving a nearby or modified optimization problem, providing a principled way to reason about training dynamics and implicit regularization [6, 13, 16], and to formalize stability with respect to structured perturbations of data 2 Backward error tolerance is the default stopping criterion in SciPy’s LSQR implementation [45].
TOWARDS UNIVERSAL CONVERGENCE IN LINEAR SYSTEM SOLVERS
5
and parameters [7]. Convergence analysis of iterative solvers. Even though backward error analysis has been used in the development of many iterative solvers, their convergence properties have been understood primarily in terms of forward error. This often leads to highly pessimistic guarantees where the number of iterations scales with the condition number of the problem [41, 20]. On the positive side, Axelsson and Lindskog [5] among many others have shown that by relying on polynomial approximation theory those guarantees can be improved for Krylov methods when the input matrix has clustered or isolated eigenvalues. Nevertheless, following the ideas of Nemirovsky and Yudin [34], Chou [10] showed that no deterministic algorithm can attain a universal convergence rate in forward error, since the dependence on the condition number is inevitable in the worst case. Similar hardness results apply when we restrict ourselves to PSD linear systems, and they can also be extended to all randomized algorithms [12]. Connections to continuous optimization. Our universal convergence guarantees bear some resemblance to the analysis of iterative solvers for minimizing convex objective functions, such as Gradient Descent (GD) and Nesterov’s Accelerated Gradient (NAG) method [35]. In particular, Richardson can be viewed as GD with respect to a certain convex quadratic objective (as discussed in detail in Section 3), and CG exhibits strong connections with NAG in the way it achieves acceleration [14]. In the continuous optimization literature, convergence of GD, NAG, and related methods is established in terms of the rate of decay of the so-called excess objective function value (relative to its minimum), attaining rates of the form O(L/k) for GD and O(L/k 2 ) for NAG (a.k.a. sublinear rates, see [14] for an overview), where L is a problemdependent parameter. Similar sublinear rates have been established for the forward error of CG by Axelsson and Kaporin [4], where a different problem-dependent parameter L appears. While none of these guarantees can be directly used to attain universal convergence of either backward or forward error in linear system solvers, the idea of using regularization to establish sublinear convergence of an algorithm (as we do in Theorem 3.8) is well known in continuous optimization [1]. 3. Universal Convergence Analysis. Here, we present universal convergence analysis of backward error for several classical iterative solvers, including Richardson iteration, Conjugate Gradient, and MINRES. 3.1. Preliminaries. We start with notation and preliminaries. Notation. We use ∥ · ∥2 to denote the 2-norm (spectral norm for matrices and Euclidean norm for vectors). For a positive semidefinite (PSD) matrix M, we define √ (A) to denote the condition the M-norm, ∥x∥M := x⊤ Mx. We use κ(A) := σσmax min (A) number, where σmax and σmin are the largest and smallest singular values, respectively. Backward error. We start by formulating the standard definition of backward error for general square linear systems of the form Ax = b, where for simplicity we will assume throughout that A, b ̸= 0. Here, we consider a variant of backward error where only a perturbation of A is allowed, but not of b. Naturally, if one allows perturbations of b as well, then this can only reduce the backward error, so all our convergence guarantees also implicitly apply to that variant. Definition 3.1 (Backward error). Given A ∈ Rn×n and b, x ∈ Rn , define: berrA,b (x) := min ∆A
∥∆A∥2 ∥A∥2
s.t.
(A + ∆A)x = b.
In our analysis, we rely on the classical representation of backward error for square
6
M. DEREZIŃSKI, Y. NAKATSUKASA, AND E. REBROVA
matrices in terms of the residual error normalized by the norm of the output vector. Lemma 3.2 (Lemma 1.1, [25]). For any x ̸= 0, its backward error with respect to A ∈ Rn×n and b ∈ Rn satisfies: berrA,b (x) =
∥Ax − b∥2 . ∥A∥2 ∥x∥2
Model of computation. Since our focus is on convergence rates rather than numerical stability, the theoretical analysis is performed in exact arithmetic. Nonetheless, extension of our results to finite-precision arithmetic is an important future direction. While this is beyond the scope of this work, we note that in standard floating-point arithmetic with unit roundoff u, it is in general not possible to obtain O(u) forward error [25, Ch. 7], even in terms of the relative residual [3]. By contrast, classical (and most robust) algorithms are known to be backward stable, that is, they compute a solution with O(u) backward error for any given problem. Therefore, requiring a small backward error is arguably the only sensible goal in finite precision, and hence it makes excellent sense to analyze convergence of a method in terms of the backward error. 3.2. Richardson Iteration. We start by analyzing the backward error convergence of the classical Richardson iteration for positive semidefinite linear systems. The following result shows that after k steps of Richardson iteration with step size ∥A∥−1 2 , the backward error is bounded by 1/k, thus attaining universal convergence. Theorem 3.3 (Richardson iteration). Given an n × n PSD A and b ∈ Rn , let x0 = 0,
xk+1 = xk − η(Axk − b).
1 . Then, for every k ≥ 1 the iterates xk satisfy Fix any C ≥ 1 and let η = C∥A∥ 2
C . k Proof. We rewrite the residual vector after k iterations of Richardson as follows: berrA,b (xk ) ≤
Axk − b = A(xk−1 − η(Axk−1 − b)) − b = Axk−1 − b − ηA(Axk−1 − b) = (I − ηA)(Axk−1 − b) = −(I − ηA)k b. Next, we perform a similar calculation to derive the expression for the iterate xk : xk = (I − ηA)xk−1 + ηb = η
k−1 X
(I − ηA)i b.
i=0 ∥Axk −b∥2 We can now bound the backward error berrA,b (xk ) = ∥A∥ as follows: 2 ∥xk ∥2 k−1
X −1 ∥Axk − b∥2 C∥(I − ηA)k b∥2 = Pk−1 ≤ C (I − ηA)k (I − ηA)i , ∥A∥2 ∥xk ∥2 ∥ i=0 (I − ηA)i b∥2 2 i=0 Pk−1 where in the last step we use that the matrix i=0 (I − ηA)i ⪰ (I − ηA)0 = I is invertible. Since all of the matrices under the spectral norm commute, letting λi denote the eigenvalues of A, it follows that: (1 − ηλj )k berrA,b (xk ) ≤ C max Pk−1 . i 1≤j≤n i=0 (1 − ηλj )
TOWARDS UNIVERSAL CONVERGENCE IN LINEAR SYSTEM SOLVERS
7
Thus, it suffices to bound the last expression by 1/k for each j. Here, we consider two cases. First, suppose that λj = 0. Then, (1 − ηλj )i = 1 for each i, and so the bound follows immediately. Next, suppose that λj > 0, and let x = ηλj . Using the formula for the geometric sum, we have k−1 X
(1 − x)i =
i=0
1 − (1 − x)k , x
so it suffices to show the following elementary inequality, (3.1)
x(1 − x)k ≤ 1/k, 1 − (1 − x)k
for all 0 < x ≤ 1 and k ≥ 1. This follows from Bernoulli’s inequality by observing that (1 − x)k (1 + xk) ≤ (1 − x)k (1 + x)k = (1 − x2 )k ≤ 1. Rearranging the terms, we obtain (3.1), thereby concluding the proof. The above result bears close resemblance to classical convergence guarantees for gradient descent over convex quadratic objectives [8]. In particular, when minimizing a function ϕ(x) = 12 ∥x − x∗ ∥2A for some vector x∗ ∈ Rn and a PSD matrix A, the gradient descent update xk+1 = xk −η∇ϕ(xk ) precisely coincides with the Richardson iteration ran on A and b = Ax∗ , and the textbook gradient descent analysis (e.g., [31]) yields the following bound: (3.2)
1 C ∥xk − x∗ ∥2A = ϕ(xk ) − ϕ(x∗ ) ≤ ∥A∥2 ∥x∗ ∥22 . 2 2k
We can convert this into a bound on the residual norm, obtaining a different type of convergence rate that slightly resembles our backward error guarantee: r p C ∥Axk − b∥2 ≤ ∥A∥2 ∥xk − x∗ ∥A ≤ ∥A∥2 ∥x∗ ∥2 . k However, this conversion is not without a cost. In fact, Theorem 3.3 implies a strictly sharper guarantee, since due to the monotonicity of ∥xk ∥2 we can upper bound the ∥x∗ ∥2 -normalized residual by the backward error: (3.3)
∥Axk − b∥2 ≤ berrA,b (xk ). ∥A∥2 ∥x∗ ∥2
Nevertheless, the connection to gradient descent begs the question whether universal convergence of backward error is more generally tied to convergence of the ∥x∗ ∥normalized residual, and with it, potentially to the broader family of tools from convex optimization theory. Unfortunately, those two notions of convergence can exhibit dramatically different behavior, as explained below. First, for some positive semidefinite linear systems the ∥x∗ ∥-normalized residual is not even well-defined. Concretely, when A is rank-deficient and b falls outside of its range, then such a system is inconsistent and the solution x∗ does not exist. Yet, even in these cases backward error is well-defined and converges at a 1/k rate when running Richardson (Theorem 3.3 does not require an x∗ to exist such that Ax∗ = b).
8
M. DEREZIŃSKI, Y. NAKATSUKASA, AND E. REBROVA
Second, we show that for non-symmetric linear systems, even when the solution does exist, the reverse direction of (3.3) can fail horribly, in that, the ∥x∗ ∥2 normalized residual enjoys a conditioning-free convergence but the backward error exhibits a highly conditioning-dependent behavior. The below result shows that, if we run Richardson iteration on the normal equations, then the gap between the backward error and the ∥x∗ ∥2 -normalized residual can be as large as a factor proportional to the condition number. In particular, this suggests that the rate of increase of ∥xk ∥2 plays a crucial role in the backward error convergence. Theorem 3.4 (Richardson on normal equations). Given an invertible A ∈ Rn×n and b ∈ Rn , let x0 = 0,
xk+1 = xk − ηA⊤ (Axk − b).
1 −1 b, for every k ≥ 1 the Fix any C ≥ 1 and let η = C∥A∥ 2 . Then, denoting x∗ = A 2 iterates xk satisfy r ∥Axk − b∥2 C Cκ(A) and . berrA,b (xk ) ≤ ≤ k ∥A∥2 ∥x∗ ∥2 2k
Furthermore, for any invertible A there exists b such that for every 1 ≤ k ≤ κ2 (A)−1, r C ∥Axk − b∥2 Cκ(A) ≥ κ(A) · berrA,b (xk ) ≥ . ek 4k ∥A∥2 ∥x∗ ∥2 Proof. Following the same argument as in the proof of Theorem 3.3, we can express the backward error and residual error as follows: √ p ∥AMk x∗ ∥2 Cη∥AMk x∗ ∥2 ∥Axk − b∥2 berrA,b (xk ) = , = Cη , ∥(I − Mk )x∗ ∥2 ∥A∥2 ∥x∗ ∥2 ∥x∗ ∥2 where we let Mk = (I − ηA⊤ A)k . Bounding the errors in terms of the singular values σ1 ≥ σ2 ≥ ... ≥ σn > 0 of A, we obtain: s p √ ησi2 (1 − ησi2 )k C ησi2 (1 − ησi2 )k 1 max ≤ ≤ Cκ(A) , berrA,b (xk ) ≤ C max 2 2 k 2 k i i 1 − (1 − ησi ) ησn 1 − (1 − ησi ) k where in the last step we used (3.1). On the other hand, for the residual error, we obtain: r q p √ √ ∥Axk − b∥2 C (3.4) ≤ C max ησi2 (1 − ησi2 )2k ≤ C max x exp(−2kx) ≤ , i x≥0 ∥A∥2 ∥x∗ ∥2 2k where in the last step we used the elementary inequality that x exp(−ax) ≤ 1/a for all a, x > 0. To show that the upper bound on the backward error can be sharp up to a constant factor, let b be the left singular vector of A associated with its smallest singular value. Then, x∗ is aligned with the corresponding right singular vector, so p p √ −1 √ √ ησn2 (1 − ησn2 )k ησn2 (1 − ησn2 )k Ce berrA,b (xk ) = C ≥ C ≥ p , 1 − (1 − ησn2 )k kησn2 k ησn2
TOWARDS UNIVERSAL CONVERGENCE IN LINEAR SYSTEM SOLVERS
9
where, letting x = ησn2 , we first used Bernoulli’s inequality to get 1 − (1 − x)k ≤ kx, 1 and then we used that k ≤ κ2 (A) − 1 ≤ ησ1 2 − 1 to obtain (1 − x)k ≥ (1 − x) x −1 = n
1−x
x − x (1 + 1−x ) ≥ e−1 , with the last step being the standard lower bound on the p √ constant e. Since 1/ ησn2 = Cκ(A), we get: r r Cκ(A) 2C C berrA,b (xk ) ≥ = κ(A) . ek e2 k 2k
Applying the upper bound (3.4) and using that e2 /2 < 4 concludes the proof. √ We note that the O(1/ k) upper bound on the ∥x∗ ∥2 -normalized residual can also be obtained directly from the gradient descent guarantee (3.2). Also, the regime 1 ≤ k ≤ κ2 (A) − 1, where the backward error may be much larger than the ∥x∗ ∥2 normalized residual, is natural since when k ≥ κ2 (A), the solution norm ∥xk ∥2 becomes comparable to ∥x∗ ∥2 and both metrics start to converge at the fast exponential 1 )k . rate, (1 − κ2 (A) The failure of Richardson on normal equations, illustrated by Theorem 3.4, perhaps further highlights the remarkable nature of the universal 1/k backward error rate in the positive semidefinite case (Theorem 3.3). Can this rate be further improved, for instance by employing Krylov subspace methods? 3.3. Faster Rates via Regularized Krylov Solvers. We next show that one can obtain faster-than-Richardson universal backward error rates using existing solvers such as CG and MINRES by relying on the perturbation interpretation of backward error. The key observation is that, by the definition of backward error, we do not need to solve the original problem directly, but rather, we can solve a “perturbed” problem that is chosen to have better properties for our method. The first ingredient in our analysis is the observation that a backward error bound for one linear system can be easily converted to a bound for another system based on a nearby input matrix. Lemma 3.5. Given A, Ã ∈ Rn×n , b ∈ Rn , and ϵ ≥ 0 such that ∥A − Ã∥2 ≤ ϵ∥A∥2 , for any x ∈ Rn , we have berrA,b (x) ≤ (1 + ϵ)berrÃ,b (x) + ϵ. Proof. First, observe that ∥Ax − b∥2 = ∥Ãx − b + (A − Ã)x∥2 ≤ ∥Ãx − b∥2 + ϵ∥A∥2 ∥x∥2 . Now, using Lemma 3.2 and that ∥Ã∥2 ≤ (1 + ϵ)∥A∥2 we obtain: berrA,b (x) =
∥Ax − b∥2 ∥Ãx − b∥2 ∥Ãx − b∥2 ≤ + ϵ ≤ (1 + ϵ) + ϵ, ∥A∥2 ∥x∥2 ∥A∥2 ∥x∥2 ∥Ã∥2 ∥x∥2
which completes the proof. Lemma 3.5 implies that if we are seeking a backward error bound for a PSD linear system Ax = b, then it suffices to analyze the error for a slightly regularized system (A + ϵ∥A∥2 I)x = b. Such regularization has the effect that it improves the condition number of the system sufficiently to ensure fast convergence of standard solvers in forward error. Then, it suffices to bound the backward error in terms of the forward error as follows.
10
M. DEREZIŃSKI, Y. NAKATSUKASA, AND E. REBROVA
Lemma 3.6. Given a square consistent linear system Ax∗ = b and ϵ ∈ (0, 1), suppose that x satisfies: (3.5)
∥x − x∗ ∥M ≤ ϵ, ∥x∗ ∥M
where M = (A⊤ A)α for some α ∈ [0, 1].
Then, x satisfies the following backward error guarantee: berrA,b (x) ≤
ϵ . 1−ϵ
Remark 3.7. Setting α = 0 in (3.5) recovers a 2-norm forward error guarantee, while α = 1 corresponds to the relative residual. If A is PSD then α = 1/2 corresponds to the A-norm error. Lemma 3.6 implies that each of these forward error guarantees immediately implies a nearly matching guarantee on backward error. Proof. Let A = UΣV⊤ be the singular value decomposition of A. Then, M = VΣ V⊤ , and we can bound the residual ∥Ax − b∥2 using the M-norm error as follows: 2α
∥Ax − b∥2 = ∥UΣV⊤ (x − x∗ )∥2 = ∥UΣ1−α V⊤ VΣα V⊤ (x − x∗ )∥2 ≤ ∥UΣ1−α V⊤ ∥2 ∥x − x∗ ∥M = ∥A∥1−α ∥x − x∗ ∥M . 2 Next, we seek to provide a lower bound on the 2-norm of the vector x. We do this by first observing that: ∥x∗ ∥M ≤ ∥x − x∗ ∥M + ∥x∥M ≤ ϵ∥x∗ ∥M + ∥x∥M . Rearranging the terms, we conclude that (1 − ϵ)∥x∗ ∥M ≤ ∥x∥M ≤ ∥M∥1/2 ∥x∥2 = ∥A∥α 2 ∥x∥2 . Combining these with Lemma 3.2, we have: berrA,b (x) =
∥A∥1−α ∥x − x∗ ∥M ϵ ∥Ax − b∥2 2 . ≤ ≤ 1−α ∥A∥2 ∥x∥2 1−ϵ ∥A∥2 (1 − ϵ)∥x∗ ∥M
Putting these observations together, we show that by solving a slightly regularized version of the linear system, we can attain a better universal rate with both CG and MINRES than with Richardson. We note that the strategy described below can be applied to essentially any solver with a forward error convergence guarantee. Theorem 3.8 (Regularized CG/MINRES). Given n × n positive semidefinite A and b ∈ Rn , fix k ≥ 9. If x̂ is produced from at least k steps of CG or MINRES on input A + 2( logk k )2 ∥A∥2 I and b, then log k 2 . berrA,b (x̂) ≤ 5 k Proof. Define ϵ = 2( logk k )2 . The condition number of the matrix à = A+ ϵ∥A∥2 I can be bounded by: κ(Ã) ≤
∥A∥2 + ϵ∥A∥2 = 1 + 1/ϵ. ϵ∥A∥2
TOWARDS UNIVERSAL CONVERGENCE IN LINEAR SYSTEM SOLVERS
11
So, letting x̃∗ = Ã−1 b, CG and MINRES initialized at 0 after k or more iterations return x̂ such that [44]: (3.6)
p k 1 + 1/ϵ − 1 2k ∥x̂ − x̃∗ ∥M ≤2 p ≤ 2 exp − p , ∥x̃∗ ∥M 1 + 1/ϵ + 1 1 + 1/ϵ + 1
where we use M = Ã2 for MINRES and M = Ã for CG. We next show that given our choice of ϵ and p k, the right-hand side is bounded 1 + 1/ϵ + 1) log(2/ϵ). It is easy to by ϵ. For this, it suffices to show that 2k ≥ ( p p verify that 1 + 1/ϵ + 1 ≤ 2/ϵ for any 0 < ϵ ≤ 1/8, so a sufficient condition is p p k ≥ 2/ϵ log( 2/ϵ). By our choice of ϵ we have: p
p 2/ϵ log( 2/ϵ) =
k k log ≤ k. log k log k
In order to ensure that ϵ ≤ 1/8, it suffices to assume that k ≥ 9. Thus, we have shown that the right-hand side of (3.6) can be bounded by ϵ. Applying Lemma 3.5, and then Lemma 3.6, we conclude that berrA,b (x̂) ≤ (1 + ϵ)berrÃ,b (x̂) + ϵ ≤ (1 + ϵ)
ϵ 2ϵ +ϵ= . 1−ϵ 1−ϵ
Plugging in the expression for ϵ and using that ϵ ≤ 1/8, we obtain the claim. The approach given in Theorem 3.8 has some clear drawbacks. First, one has to commit to a target backward error at the beginning in order to choose the regularization level. To obtain a more accurate estimate, one would have to restart CG or MINRES with a different regularizer. Second, the indirect nature of this procedure inevitably incurs some overhead, which is manifested by the spurious logarithmic factors. We address all of these issues in the next section by proposing a dedicated Krylov solver for backward error. 4. MINBERR: Optimal Backward Error Krylov Solver. Here, we develop and analyze MINBERR (Algorithm 4.1), a Krylov method for a PSD linear system Ax = b, that is explicitly designed to minimize the backward error. That is, letting Kk (A, b) = span{b, Ab, ..., Ak−1 b} be the Krylov subspace of rank k, consider (4.1)
xk := argmin berrA,b (x). x∈Kk (A,b)
We show that this sequence attains a faster universal convergence rate than other solvers. We also extend the method to non-symmetric systems (MINBERR-NE, Algorithm 4.2), and provide efficient implementations for both algorithms. 4.1. MINBERR: Convergence Analysis. First, we show that the MINBERR iterate sequence (4.1) attains a universal backward error convergence rate of O(1/k 2 ). Theorem 4.1 (MINBERR). For any n × n PSD matrix A, vector b ∈ Rn and k ≥ 2, the MINBERR iterate sequence (4.1) satisfies: berrA,b (xk ) ≤
3 k2 − 1
.
Remark 4.2. Numerical evidence suggests that the minimax optimal constant is 2.
12
M. DEREZIŃSKI, Y. NAKATSUKASA, AND E. REBROVA
Proof. First, define Ā = A/∥A∥2 and b̄ = b/∥A∥2 . Note that Kk (Ā, b̄) = Kk (A, b) and also berrĀ,b̄ (x) = berrA,b (x) for all x, so without loss of generality we can assume that ∥A∥2 = 1. Next, let Pk−1 denote the set of all polynomials p(x) = c0 + c1 x + · · · + ck−1 xk−1 of degree at most k − 1 with real-valued coefficients. Then, using Lemma 3.2, min x∈Kk (A,b)
berrA,b (x) = min
p∈Pk−1
∥Ap(A)b − b∥2 ≤ min ∥(Ap(A) − I)p(A)−1 ∥2 . p∈Pk−1 ∥p(A)b∥2
Since A and p(A) commute, we can express the above in terms of A’s eigenvalues 1 = λ1 ≥ ... ≥ λn ≥ 0: (Ap(A) − I)p(A)−1 2 = max i
1 λi p(λi ) − 1 ≤ max x − . p(λi ) p(x) x∈[0,1]
It remains to show that the latter term is upper-bounded by 3/(k 2 − 1). We note that this order is tight: actually, both lower and upper bounds of O(1/k 2 ) can be concluded from a more general polynomial approximation result of Levin and Saff from 1988 [32], who obtained it in a context unrelated to Krylov subspace theory: Lemma 4.3 (Theorem 2.3, [32]). For any α > 0, there exist positive constants Bα , Cα such that for any k = 1, 2, ..., the following holds: Bα 1 Cα ≤ min max xα − ≤ 2α . p∈Pk x∈[0,1] k 2α p(x) k Levin and Saff do not provide explicit values for Bα and Cα . To address this, we present a direct construction of p(x) which shows that Cα ≤ 3 when α = 1. First, for x ∈ [0, 1] and any positive integer m, define the shifted Chebyshev ∗ ∗ (x) = Tm (2x − 1), where Tm (x) = cos(m arccos x). For any (x) as Tm polynomial Tm k ≥ 2, let ℓ be an even integer ℓ ∈ {k, k + 1}, and define (4.2)
G(x) :=
1 − Tℓ∗ (x) 2ℓ2
for x ∈ [0, 1].
The following key facts about the function G(x) follow from trigonometric identities. Lemma 4.4 (Properties of shifted Chebyshev polynomials). Let G(x) be defined as in (4.2). Then, it holds that G(x) =
sin2 (ℓγ) for γ ∈ [0, π/2] such that x = sin2 (γ). ℓ2
This implies that (a) 0 ≤ G(x) ≤ x for every x ∈ [0, 1], and (b) the function G(x) − x has a root of multiplicity at least 2 at x = 0. Lemma 4.4 (proven in Appendix C) implies that (4.3)
p(x) :=
x − G(x) x2
is a well-defined polynomial of degree at most ℓ − 2 ≤ k − 1 and for any x ∈ [0, 1] (4.4)
x−
1 x2 xG(x) 1 ℓ2 sin2 (γ) sin2 (ℓγ) = x− = = 2· 2 2 , p(x) x − G(x) x − G(x) ℓ ℓ sin (γ) − sin2 (ℓγ)
where x = sin2 (γ) as in Lemma 4.4. We maximize this expression using basic calculus.
TOWARDS UNIVERSAL CONVERGENCE IN LINEAR SYSTEM SOLVERS
13
Lemma 4.5 (Explicit upper bound). For every integer ℓ ≥ 2, define Fℓ (γ) :=
ℓ2 sin2 (γ) sin2 (ℓγ)
, ℓ2 sin2 (γ) − sin2 (ℓγ)
γ ∈ [0, π/2].
Then Fℓ is well defined on [0, π/2] after continuous extension at γ = 0, and max Fℓ (γ) =
γ∈[0,π/2]
3ℓ2 ℓ2 − 1
.
The proof is in Appendix C. Since ℓ ≥ k, together with (4.4) this immediately implies: max x − x∈[0,1]
1 3 3 ≤ 2 . ≤ 2 p(x) ℓ −1 k −1
4.2. MINBERR: Implementation and Complexity. Here, we describe how to compute the MINBERR solution xk efficiently √ (and approximately), and show that MINBERR has complexity bounded by O(n2 / ϵ) to achieve ϵ backward error. First we show that finding the minimizer of the backward error in a subspace is equivalent to solving a certain generalized eigenvalue problem. Indeed, when minimizing berrA,b (x) over Kk (A, b), we seek to find a solution for Ax = b within a given k-dimensional subspace span(Qk ), parameterized by Qk ∈ Rn×k having orthonor∥Ax−b∥2 mal columns. If we write x = Qk y, then minimizing berrA,b (x) = ∥A∥ over 2 ∥x∥2 x ∈ Kk (A, b) reduces to minimizing the following generalized Rayleigh quotient: ⊤ −1 −1 y⊤ b AQk b AQk y ∥AQk y − b∥2 argmin = argmin ⊤ . ∥Q y∥ k k 0 0 −1 k 2 y∈R y∈R −1 y⊤ Ik Ik y The value of the minimum is squared in the second expression. Note that if b is in the span of AQk then MINBERR finds the exact solution to Ax = b, for which trivially berrA,b (x) = 0. Otherwise (in the usual case), b AQk has full column ⊤ b AQk is positive definite. In the latter case, minimizing rank and so b AQk such quotient is equivalent [19, Sec. 8.7.4] to finding the smallest eigenpair of the generalized eigenvalue problem ⊤ ⊤ −1 0 0 −1 b AQ b AQ (4.5) =λ , k k y Ik Ik y which is ∥b∥22 (AQk )⊤ b
(4.6)
b⊤ AQk (AQk )⊤ AQk
−1 0 =λ y
−1 . Ik y
This is a positive (semi)definite generalized eigenvalue problem, and has real eigenvalues, with one at infinity. The smallest eigenvalue λmin , which is nonnegative, is equal to the square of the backward error of the MINBERR solution, that is, λmin = (berrA,b (xk ))2 . The argument so far holds for an arbitrary subspace Qk . We note that a similar process was described by Kasenally and Simoncini [27, 28], who take Qk to be
14
M. DEREZIŃSKI, Y. NAKATSUKASA, AND E. REBROVA
the Krylov subspace for a non-symmetric A found by the Arnoldi decomposition [41, Ch. 6], and solve such generalized eigenvalue problem, which is unstructured except for the symmetry. Below we exploit the tridiagonal structure in the Lanczos decomposition to obtain a much more efficient algorithm when A is symmetric. Later, in Section 4.3 we extend this approach to non-symmetric A by working in the Krylov subspace corresponding to the normal equations, obtaining MINBERR-NE. 4.2.1. Fast implementation by exploiting symmetry. Let us assume that A is symmetric (PSD is not strictly required). In this case, Krylov subspace methods start with the Lanczos decomposition [41, Sec. 6.6]: (4.7)
AQk = Qk+1 Tk
where Qk = [q1 , . . . , qk ] is n × k, Qk+1 = [q1 , . . . , qk+1 ] is n × (k + 1), both having orthonormal columns, and Tk is (k +1)×k symmetric tridiagonal. The columns of Qk represent the orthonormal basis of the Krylov subspace at iteration k. By the standard Lanczos construction b = Qk+1 ek+1 ∥b∥2 where ek+1 = [1, 0, . . . , 0]⊤ ∈ Rk+1 , so that [b AQk ] = Qk+1 [ek+1 ∥b∥2 Tk ]. Thus the equation (4.5) simplifies to −1 0 −1 ⊤ =λ , where M = ∥b∥2 ek+1 Tk . (4.8) M M y Ik y M is a (k + 1) × (k + 1) tridiagonal and upper triangular matrix. It thus suffices to find the smallest eigenpair of the generalized eigenvalue problem (4.8). The next result shows how to do so efficiently and approximately. Lemma 4.6. Let (σmin , vmin ) be the smallest singular value and the corresponding right singular vector of the k × k matrix T̃k obtained by removing the first row from Tk in (4.7). Then, the smallest eigenpair of the generalized eigenvalue problem (4.8) 1 2 ((Tk )1,1 (vmin )1 + (Tk )1,2 (vmin )2 ). Conis given by (σmin , α1 vmin ), where α = − ∥b∥ 2 sequently, the MINBERR solution xk after k iterations satisfies: xk =
1 Qk vmin , α
berrA,b (xk ) = σmin .
1 b= b ∈ Rk , defining α ((Tk )1,1 (b v)1 + (Tk )1,2 (b v)2 ) and x Moreover, for any v b = − ∥b∥ 2 1 b b ∥2 /∥b x) = ∥T̃k v v ∥2 . α Qk v, we have berrA,b (b
Proof. Since the second matrix in (4.8) is singular, there is an eigenvalue of (4.8) at infinity λ∞ = ∞. As λ∞ is irrelevant to the MINBERR solution, we deflate, or remove, λ∞ by employing operation with respect to a nonsingular W a congruence 0 ⊤ ⊤ such that W (M M, )W is simultaneously block diagonal, as follows: We Ik choose W so that left-multiplying to M adds multiples of the first column of M to other columns so as to eliminate the other two nonzero elements in the first row of M. This is achieved by setting 1 w1 W= ∈ R(k+1)×(k+1) , where w1 = 1, −(Tk )1,1 , −(Tk )1,2 , 0, . . . , 0 . 0 Ik ∥b∥2 We then obtain 0 ⊤ ⊤ (4.9) W (M M),
Ik
W=
1
0 T̃k
⊤ 1
0 0 , T̃k
! Ik
.
TOWARDS UNIVERSAL CONVERGENCE IN LINEAR SYSTEM SOLVERS
15
Since the eigenvalues of a matrix pair are invariant under congruence (or more generally equivalence) transformations [20, Sec. 7.7], to find the MINBERR solution it suffices to examine the eigenvalues of (4.9), where the eigenvalue λ∞ is decoupled from the rest. Further, since the second matrix in (4.9) has Ik in the second block, we see that the smallest eigenvalue of the generalized eigenvalue problem (4.8) is equal to that of T̃⊤ k T̃k (i.e., a standard symmetric eigenvalue problem), which in turn can be reduced to the SVD of T̃k ∈ Rk×k ; more precisely, σmin(T̃k ) = berrA,b (xk ),and the 0 1 corresponding right singular vector vmin is such that W is parallel to ; this vmin y 0 follows from the congruence relation between (4.8) and (4.9). Now since W = v min α 0 1 ((Tk )1,1 (vmin )1 + (Tk )1,2 (vmin )2 ), we can where α = w1 = − ∥b∥ 2 vmin vmin compute3 y = α1 vmin . We obtain the MINBERR solution as xk = Qk y. b ∈Rk , letting b=Q b For the latter x ky statement, first note from (4.2) that for any y 0 0 −1 = berrA,b (b x). Note also from (4.9) that = MW . we have M b b 2 v y b T̃k v 0 1 Since α b = − ∥b∥ ((Tk )1,1 (b v)1 + (Tk )1,2 (b v)2 ) is such that the first element of αb1 W 2 b v −1 0 b = αb1 v. is −1, that is, αb1 W = 1 , we obtain the desired result by taking y b v v α b This shows that the MINBERR solution can be found via the SVD of T̃k , and that one can find an approximate minimizer of berrA,b (x) over Qk by finding a unit-norm vector v that approximately minimizes ∥T̃k v∥2 /∥v∥2 . Let us consider a common situation where a prescribed tolerance ϵ is required. The MINBERR computation therefore proceeds as follows: First, run the Krylov iterations until σmin (T̃k ) < ϵ. Then find the corresponding singular vector vmin . In the first phase, given a required tolerance ϵ for the backward error, our algorithm first tests if σmin (T̃k ) is smaller than ϵ by attempting a Cholesky factorization 2 ⊤ of the symmetric pentadiagonal matrix T̃⊤ k T̃k − ϵ Ik = R R. Cholesky succeeds if and only if σmin (T̃k ) > ϵ holds [25, Ch. 10]4 . Note that, at each iteration k, T̃k grows in size by 1, and only its final column needs computing (more specifically the bottom two elements). Similarly, the Gram matrix G̃ = T̃⊤ k T̃k and the Cholesky factor R both need to be computed only in the final columns, resulting in O(1) cost overall for the convergence check. If Cholesky breaks down, implying σmin (T̃k ) < ϵ, we terminate the Krylov iterations and proceed to compute vmin . This can be done by running the inverse −1 iteration, i.e., set v to be a random initial guess, and repeat: v := (T̃⊤ v, and k T̃k ) v := v/∥v∥2 . We then have v → vmin as we run more iterations. Each application −1 of (T̃⊤ can be done by solving two linear systems with respect to T̃⊤ k T̃k ) k and T̃k . These are tridiagonal and triangular, so each linear solve is O(k) operations. In practice, we observe that the inverse iteration typically converges to high accuracy in one 3 Unless α = 0; which happens precisely when b is orthogonal to the entire subspace AQ , as we k b see by examining (4.2). Note that when A ≻ 0 this is impossible, as b⊤ Ab > 0 and ∥b∥ is the first 2 column of Qk . When A is PSD, the same holds unless Ab = 0, that is, b lies in the null space of A. 4 The process is backward stable [25]; nonetheless, breakdown of Cholesky can occur if G̃ − ϵ2 I k is numerically rank deficient; which can happen even when σmin (T̃) > ϵ if ϵ2 = O(u), √ where u is the unit roundoff. This can only happen when we require backward error smaller than u. In this case, one is advised to go to step 7 and compute the backward error by inverse iteration.
16
M. DEREZIŃSKI, Y. NAKATSUKASA, AND E. REBROVA
or two steps for large k, which is unsurprising, since it is applied to a nearly-singular matrix with smallest singular value < ϵ. We summarize the implementation in Algorithm 4.1. Algorithm 4.1 MINBERR: Solver for symmetric positive semidefinite linear systems Input: A, b, iterations kmax , backward error tolerance ϵ, failure probability δ bk such that either k = kmax or berrA,b (xk ) < ϵ Output: x 1: Set q1 = b/∥b∥2 , q0 = 0, β1 = 1, ℓ = 0 2: for k = 1, 2, ..., kmax do 3: (Lanczos:) ŵk = Aqk −βk qk−1 , αk = ŵk⊤ qk , wk = ŵk −αk qk , βk+1 = ∥wk ∥2 , qk+1 = wk /βk+1 . This executes the Lanczos process [41, Sec. 6.6], giving AQk = Qk+1 T, where T ∈ R(k+1)×k is symmetric tridiagonal with Ti,i = αi , Ti,i+1 = βi+1 . β2 α2 β3 β3 α3 β4 .. .. .. . . . by appending the last column. 4: Update T̃ := . . .. .. βk βk αk βk+1 ⊤ Also update G̃ = T̃ T̃ by appending the last row and column. 5: (Check for ϵ convergence:) Attempt Cholesky factorization of G̃−ϵ2 Ik = R⊤ R (by computing the three nonzero entries in R’s last column); if this breaks down (i.e., σmin (T̃) = berrA,b (xk ) < ϵ), proceed to step 7. Otherwise return to step 2. 6: end for 7: Find smallest singular vector of T̃ by inverse iteration: Let v be Gaussian, and 8: while ℓ < 2.23 ln(k/δ 2 ) do 9: Solve T̃⊤ z = v for z, and T̃v = z for v. 10: Set v := v/∥v∥2 , ℓ := ℓ + 1. 11: end while 1 bk = α1 Qv where α = − ∥b∥ ((Tk )1,1 (v)1 + (Tk )1,2 (v)2 ). 12: x 2 4.2.2. Time Complexity Analysis. Next, we put together the above analysis to obtain the overall computational cost of Algorithm 4.1, and show that it does in fact approximately minimize backward error in the Krylov subspace. Theorem 4.7. Let A ≻ 0 be n × n and b ∈ Rn . Running k steps of MINBERR bk such that berrA,b (b (Algorithm 4.1) yields x xk ) ≤ (1 + 21 ) minx∈Kk (A,b) berrA,b (x) with probability at least 1 − δ, and takes O k(TA + n + log( 1δ )) operations. Proof. We look at the costs for each step of Algorithm 4.1. Step 3: one clearly needs kTA time for the k matrix-vector multiplications with A. The O(k) vector inner-products required to form the Lanczos vectors amount to O(nk) operations. Step 4 is O(1) cost, as there are only three entries to be added to T̃ and G̃. The same holds for step 5, as the only quantity that needs computing is the final column of R, which has three nonzeros. We now examine the cost of approximately computing vmin using inverse iteration, steps 7–11. To do so, we need to combine Lemma 4.6 with a result on the number of inverse iterations required to obtain a gap-independent constant-factor relative bound on the Rayleigh quotient for the smallest eigenvalue of a positive definite matrix
TOWARDS UNIVERSAL CONVERGENCE IN LINEAR SYSTEM SOLVERS
17
(here, T̃⊤ T̃). This follows analogously to the classical analysis of Kuczyński and Woźniakowski [29] for standard power iteration (see Appendix A). Lemma 4.8. Let M be n×n positive definite with smallest eigenvalue λmin , and let h 32n ln(2/δ) i ln( πδ2 ∆ ) v0 consist of n Gaussian entries. Define vℓ = (M−1 )ℓ v0 . If ℓ ≥ 21 1 + ln(1+∆/2) for ∆ ∈ (0, 1), then with probability 1 − δ we have vℓ⊤ Mvℓ /∥vℓ ∥22 ≤ (1 + ∆)λmin . Applying Lemma 4.8 to T̃⊤ T̃ shows that ℓ ≥ 2.23 ln(k/δ 2 ) steps of inverse iteration suffice to get ∆ = 1/2 relative accuracy with probability at least 1 − δ. Finally, computing xk = α1 Qk v costs O(nk) operations. Thus the total operation count is O k(TA + n + log( 1δ ) + log k) = O k(TA + n + log( 1δ )) , as k ≤ n. Let us discuss the case when A is PSD but singular, as the convergence theory in Theorem 4.1 allows it. In this case T̃k is full rank until k is equal P to the number n of nonzero distinct components ñ in the eigenvector expansion b = i=0 ci vi , and the proof remains unaffected. Once k = ñ, Qk contains a null vector of A, which can be scaled arbitrarily large to yield an arbitrarily small backward error. We note that inverse iteration would then break down in exact arithmetic, but in finite precision converges with excellent accuracy [39], often in one step. Alternatively, one can run MINBERR on A + ϵ∥A∥2 I ≻ 0 and rely on Lemma 3.5 to bound the backward error. By combining the O(1/k 2 ) convergence guarantee from Theorem 4.1 with the complexity analysis in Theorem 4.7, we immediately obtain the following. n Corollary 4.9. Let A ⪰ 0 be n × n and b ∈ R . MINBERR finds a solution 1 1 with ϵ backward error in O( √ϵ TA + n + log( δ ) time with probability at least 1 − δ.
Remark 4.10. Under the very mild assumption that TA ≥ n, this implies √ that we can find an ϵ backward error solution with probability 1 − 2−Cn in O(TA / ϵ) time. We note that other natural variants of MINBERR may be considered: p • If we fix the number of iterations upfront (e.g., to be ⌈ 2/ϵ⌉, which guarantees that backward error is bounded by ϵ), then we can bypass steps 4–5. • In our experiments, to track the convergence we ran steps 7–12 at every iteration, which incurs an additional cost of O(k 2 log(1/δ)). 4.3. MINBERR-NE: Extension to Non-Symmetric Systems. In the nonsymmetric case there are two approaches to forming the subspace Span(Qk ): GMRESstyle (standard Krylov), and LSQR-style (Krylov for the normal equations). With the standard Krylov subspace Span(b, Ab, A2 b...Ak−1 b), a well-known example, also discussed by Kasenally [27], A=
1
In−1
and
b = [0, 0, . . . , 0, 1]⊤
gives a subspace orthogonal to the solution x = [1, 0, . . . , 0, 0]⊤ until k = n, and hence backward error ≥ 1. This example shows that a universal convergence rate is not possible with the standard Krylov subspace. More generally, the convergence of GMRES is still not fully understood [22, 23], and not fully characterized by the eigenvalues or singular values. For these reasons here we opt for the Krylov subspace of the normal equations, Span(A⊤ b, (A⊤ A)A⊤ b, . . . , (A⊤ A)k−1 A⊤ b), resulting in a bidiagonalization-based decomposition (4.10)
AQk = Uk Bk ,
18
M. DEREZIŃSKI, Y. NAKATSUKASA, AND E. REBROVA
where Qk is n × k, Uk is n × (k + 1) with first column b/∥b∥2 , Bk is (k + 1) × k lower bidiagonal, and Qk , Uk both have orthonormal columns. This requires k matrixvector multiplications with A and A⊤ , as implemented for example in LSQR [37]. The orthogonalization cost is O(nk) as it needs to be done only against the current and last vectors. This also yields the Lanczos decomposition for A⊤ A via A⊤ AQk = Qk Tk + qk+1 [0, . . . , γk+1 ] where Tk = B⊤ k Bk is tridiagonal, and γk+1 = αk βk+1 . One can plug these identities into (4.6) to derive an algorithm analogous to MINBERR, which we call MINBERR-NE (Algorithm 4.2). Remarkably, the bidiagonal structure allows us to obtain an even more stable and efficient implementation. Algorithm 4.2 MINBERR-NE: Solver for non-symmetric linear systems Input: A, b, iterations kmax , backward error tolerance ϵ, failure probability δ Output: xk such that either k = kmax or berrA,b (xk ) < ϵ 1: β1 = ∥b∥2 , u1 = b/β1 , α1 q1 = A⊤ u1 s.t. ∥u1 ∥2 = ∥q1 ∥2 = 1. 2: for k = 1, 2, ..., kmax do 3: (Bidiagonalize:) Compute βk+1 uk+1 = Aqk − αk uk , and αk+1 qk+1 = A⊤ uk+1 − βk+1 qk , s.t. ∥uk+1 ∥2 = ∥qk+1 ∥2 = 1. This yields AQk = Uk Bk . 4: (One dqds step to test ϵ-convergence:) Set d := p1 − ϵ2 (only when k = 1). Let pk = (Bk )2k+1,k , ek−1 = (Bk )2k,k . Set p̂k−1 = d + ek−1 , and d := d · pk /p̂k−1 − ϵ2 . 5: If d > 0, return to step 2. Else, σmin (B̃k ) = berrA,b (xk ) < ϵ, so go to step 7. 6: end for 7: Set B̃k as the bottom k × k submatrix of Bk . Let v be Gaussian. 8: while ℓ < 2.23 ln(k/δ 2 ) do 9: Solve B̃⊤ k z = v for z, and B̃k v = z for v. 10: Set v := v/∥v∥2 , ℓ := ℓ + 1. 11: end while 1 1 12: xk = α Qv where α = − ∥b∥ (Bk )1,1 (v)1 . 2 Proposition 4.11. Let (σmin , vmin ) be the smallest singular value and the corresponding right singular vector of the k × k bidiagonal matrix B̃k obtained by remov1 ing the first row from Bk in (4.10). Then, defining α = − ∥b∥ (Bk )1,1 (vmin )1 , the 2 MINBERR-NE solution after k iterations satisfies xk := argmin berrA,b (x) = x∈Span(Qk )
1 Qk vmin , α
berrA,b (xk ) = σmin .
1 b ∈ Rn×k , defining5 α b = αb1 Qk v b , we Moreover, for any v b = − ∥b∥ (Bk )1,1 (v)1 and x 2
b ∥2 /∥b have berrA,b (b x) = ∥B̃k v v∥2 . The proof of Proposition 4.11 is similar to Lemma 4.6 and deferred to Appendix B. A very efficient way to check if σmin (B̃k ) < ϵ for a prescribed tolerance ϵ is to run a single step of the dqds algorithm [17] with shift ϵ, which is equivalent to performing the Cholesky factorization but can preserve relative high accuracy. Once convergence is detected so that σmin (B̃k ) < ϵ, we terminate the MINBERRNE iterations and compute the singular vector corresponding to the smallest singular value of B̃k . As for A ≻ 0, we can do this by inverse iteration with respect to B̃⊤ k B̃k . 5 Unlike the PSD case, α = 0 is possible here; in this case, the MINBERR solution does not exist, with the infimum of berrA,b being 1, where x approaches an ’infinite’ multiple of a vector parallel to the smallest singular vector of AQk with a residual norm ≫ ∥b∥2 , which is hardly useful.
19
TOWARDS UNIVERSAL CONVERGENCE IN LINEAR SYSTEM SOLVERS
(a) Pres Poisson data:
(b) olafu data:
(c) raefsky4 data:
PSD matrix, n = 14822 (Fluid Dynamics)
PSD matrix, n = 16146 (Aerospace Structural Model)
PSD matrix, n = 19779 (Structural Stability Analysis)
Backward error
102 100 10!2
1 k
10!4 10!6 100
1 k2
101
102
103
104
Iteration k
105
102 100
100
Backward error
Richardson LSMR LSQR MINBERR-NE
Backward error
104
10!2 1 k
10!4 10
!6
10
!8
100
Richardson LSMR LSQR MINBERR-NE 101
102
103
104
105
1 k
10!4 10!8 10
!12
10!16 100
1 k2
Richardson LSMR LSQR MINBERR-NE 101
Iteration k
102
103
104
105
Iteration k
(d) sherman3 data:
(e) bayer03 data:
(f) cyl6 data:
square non-PSD, n = 5005 (Fluid Dynamics)
square non-PSD, n = 6747 (Chemical Simulation)
square non-PSD, n = 13681 (Shell Mechanics Simulation)
Fig. 1: Empirical evaluation of MINBERR (a,b,c) and MINBERR-NE (d,e,f) on benchmark problems [11] with PSD and square non-PSD matrices, respectively.
−1 To compute (B̃⊤ v we solve B̃⊤ k B̃k ) k z = v for z, then B̃k v = z; each is a bidiagonal linear system that can be solved in O(k) operations. Finally, as before, once vmin is (approximately) obtained, one can find the MINBERR-NE solution xk = α1 Qk vmin , 1 where α = − ∥b∥ (Bk )1,1 (vmin )1 . See Algorithm 4.2 for details. 2
5. Numerical Experiments. Here, we provide numerical experiments6 supporting our analysis, comparing MINBERR and MINBERR-NE to popular methods on several test problems. We first illustrate the typical backward error convergence behavior of different methods on real-world tasks, both with PSD and general matrices. We then perform a more in-depth investigation of specific phenomena by evaluating individual methods on carefully designed synthetic hard problems. 5.1. Benchmark Problems from Real Data. We compare the performance of our proposed algorithms against that of several baselines on six problems from the SuiteSparse Matrix Collection [11]. PSD problems. In Figure 1 (a,b,c), we compare MINBERR against Richardson, CG, and MINRES on three PSD problems. For each method, we plot the backward error against iterations, and we also include lines for the two universal rates attained by our theory: 1/k (proven for Richardson) and 1/k 2 (proven for MINBERR up to a constant factor). We observe that Richardson closely follows the theoretical rate, 6 The MATLAB code is publicly available at https://github.com/nakatsukasayuji/MINBERR.
20
M. DEREZIŃSKI, Y. NAKATSUKASA, AND E. REBROVA
10 0
10 0
10 0
10 -1
10 -1
10 -1
10 -2
10 -2
10 -2
10 -3 10 0
10 2
10 -3 10 0
10 2
10 -3 10 0
10 2
Fig. 2: Comparing MINRES with MINBERR on the Small-Outlier(2000, κ, σn−1 ) synthetic task, while varying κ and σn−1 . MINBERR is denoted by solid lines, while MINRES is shown using dotted lines.
suggesting that the bound in Theorem 3.3 is often tight for typical ill-conditioned problems. On the other hand, MINBERR often outperforms its theoretical 1/k 2 rate, especially later in the convergence. The behavior of CG and MINRES varies substantially between the problems, but they consistently outperform Richardson, and are consistently outperformed by MINBERR. Non-Symmetric Problems. Next, in Figure 1 (d,e,f) we compare MINBERR-NE against Richardson-NE (Richardson on normal equations), as well as LSQR and LSMR, on three general linear system tasks. First, we observe that all three baselines exhibit some initial blow-up of backward error, in that the convergence curves start far above 1. This phenomenon aligns with our analysis of Richardson on normal equations (Theorem 3.4), and is further demonstrated below on synthetic data. MINBERR-NE does not exhibit any initial blow-up and tends to converge at least at a 1/k rate, but clearly does not match the 1/k 2 rate from the PSD setting. 5.2. Hard Problems from Synthetic Matrices. We next evaluate the algorithms on synthetic problems which are designed to induce their worst-case performance, and further illustrate our theoretical analysis. We define the following classes of problem instances: 1. Ill-Conditioned(n, κ). Let A be an n × n diagonal matrix with logarithmically spaced entries between 1 and 1/κ, and let b = [1, ..., 1, κ]⊤ ∈ Rn consist of all 1’s except for the last entry which is equal κ. 2. Small-Outlier(n, κ, σn−1 ). Let A be an n × n diagonal matrix with first n − 1 entries logarithmically spaced √ between 1 and σn−1 , and the last entry equal 1/κ. Also, we let b = [1, ..., 1, n]⊤ ∈ Rn . We can extend these constructions to non-diagonal PSD matrices by substituting A → UAU⊤ , or to non-symmetric matrices with A → UAV⊤ , and replacing b with Ub, for some orthogonal matrices U, V. Our PSD/non-symmetric experiments are invariant under those transformations (up to numerical precision). In all of the following experiments we use problem dimension n = 2000. MINRES vs MINBERR. We first demonstrate that popular Krylov solvers on PSD systems, even though often effective in practice, may exhibit a backward error rate no better than Richardson in the worst case. Here, we illustrate this by running MINRES on the Small-Outlier problem while varying κ and σn−1 . In Figure 2, we observe that by adjusting these problem parameters we can force MINRES (dotted lines) to converge nearly as slowly as 1/k for an extended number of iterations.
TOWARDS UNIVERSAL CONVERGENCE IN LINEAR SYSTEM SOLVERS
10 15
10 15
10 15
10 10
10 10
10 10
10 5
10 5
10 5
10 0
10 0
10 0
10 -5 10
10 -5 0
10
1
10
2
10
3
10
21
10 -5 0
10
1
10
2
10
3
10 0
10 1
10 2
10 3
Fig. 3: Comparing MINBERR-NE against Richardson-NE, LSMR, and LSQR on the Ill-Conditioned(2000, κ) synthetic task treated as a general non-symmetric system.
Specifically, the results suggest that for any sufficiently large k and κ, we can find a problem instance for which the backward error of MINRES after k iterations is observed to be at least min{1/k, log10 (κ)/(2k 2 )}. By increasing κ, we can force the 1/k rate to dominate, although this is limited by the mild logarithmic dependence on κ observed for the MINRES backward error. We note that MINBERR (solid lines) is not affected by this phenomenon, avoiding condition number dependence altogether. Error Blow-up on Normal Equations. Turning to non-symmetric systems, we next demonstrate that the initial backward error for Richardson-NE/LSQR/LSMR can indeed be as large as the condition number of the matrix, as suggested by Theorem 3.4. In Figure 3, we illustrate this on the Ill-Conditioned problem with a range of κ values. We observe that the backward error blow-up is nearly as severe for LSQR and LSMR as it is for Richardson, although the Krylov solvers exhibit faster convergence after that. Finally, similarly as in all of the benchmark problems, MINBERR-NE does not exhibit any backward error blow-up and its convergence curves consistently lie below the 1/k threshold. Stagnation of MINBERR-NE. Unfortunately, it turns out that MINBERR-NE does not always attain the universal 1/k rate. We show this empirically in Figure 4 (top plots), again using the Small-Outlier family of problems. Here, by adjusting the parameters we can force the backward error of MINBERR-NE to stagnate close to 1 for an extended number of iterations. Specifically, our experiments suggest that for any sufficiently large k and κ, we can find a problem instance for which the backward error of MINBERR-NE after k iterations is at least min{1, log10 (κ)/k}. As we can see in Figure 4 (dotted lines), the MINBERR-NE backward error exhibits an initial 1/k convergence until reaching a plateau around the value of σn−1 , and then eventually it resumes fast convergence. By taking the outlier singular value to zero (blowing up κ), we can extend the plateau, making backward error stagnate for as long as we want (up to the numerical limitations that arise when increasing the condition number). We note that this phenomenon also occurs if we replace the single small outlier with a large cluster of small singular values, or if we make the problem non-symmetric by applying different orthogonal matrices U and V on each side. Perturbed MINBERR-NE. The above example suggests that a universal rate is not attainable for any algorithm that relies on the Krylov subspace defined by the normal equations, due to the logarithmic dependence on the condition number of the input matrix. However, we can get around this issue by perturbing the matrix to reduce its condition number. In the psd case we could do this by simply regularizing
22
M. DEREZIŃSKI, Y. NAKATSUKASA, AND E. REBROVA
10 0
10 0
10 0
10 -1
10 -1
10 -1
10 -2
10 -2
10 -2
10 0
10 2
10 0
10 2
10 0
10 2
Fig. 4: Evaluating MINBERR-NE on the Small-Outlier(2000, κ, σn−1 ) synthetic task treated as a general non-symmetric linear system. We use dotted lines to denote MINBERR-NE running on the original task, while solid lines denote MINBERR-NE running on the matrix perturbed with Gaussian noise (in both cases, backward error is computed with respect to the original task).
with a scaled identity (see Section 3), but in the general case we must resort to a random Gaussian perturbation. Namely, in Figure 4 (solid lines), we replace A with ∥A∥2 à = A + ϵ ∥G∥ G, where G has independent Gaussian entries and we let ϵ = 10−3 . 2 We then run MINBERR-NE on à and b, while still evaluating its backward error with respect to the original matrix A. Lemma 3.5 ensures that the perturbation will not introduce more than ≈ ϵ additional error (which is negligible in our experiment), and classical results in smoothed analysis [42, 9] show that the new matrix à has condition number O(n) with high probability. This is confirmed in our experiments (Figure 4, solid lines), which show that Perturbed MINBERR-NE is insensitive to increasing the condition number of the problem. 6. Conclusions and Future Directions. We showed that, when measured in terms of backward error, the convergence of iterative linear system solvers can be made universal, i.e., independent of the conditioning of the problem. For PSD linear systems, we gave a universal 1/k backward error convergence guarantee for the classical Richardson iteration. We then developed an efficient Krylov subspace method, MINBERR, which attains an even better O(1/k 2 ) universal rate and has strong empirical performance when compared against CG and MINRES on benchmark problems. For general linear systems, we extended our algorithm via the normal equations, obtaining MINBERR-NE which empirically attains O(1/k) backward error on highly ill-conditioned benchmark problems. Our synthetic experiments suggest that even on artificially constructed hard instances, a small perturbation of the input suffices to ensure a universal convergence rate of the method. This raises the question whether a smoothed analysis of MINBERR-NE can reveal a truly universal convergence in backward error for all linear systems. Another important future direction is to provide a finite-precision backward error analysis of MINBERR and MINBERR-NE, so that universal convergence rates can be guaranteed in a numerically stable way. Acknowledgments. The authors thank Anil Damle, Ethan Epperly, Raphael Meyer, Cameron Musco, and Christopher Musco for helpful conversations. REFERENCES
TOWARDS UNIVERSAL CONVERGENCE IN LINEAR SYSTEM SOLVERS
23
[1] Z. Allen-Zhu and E. Hazan, Optimal black-box reductions between optimization objectives, Advances in Neural Information Processing Systems, 29 (2016). [2] J. Alman, R. Duan, V. V. Williams, Y. Xu, Z. Xu, and R. Zhou, More asymmetry yields faster matrix multiplication, in Proceedings of the 2025 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), SIAM, 2025, pp. 2005–2039. [3] M. Arioli, I. Duff, and D. Ruiz, Stopping criteria for iterative solvers, SIAM J. Matrix Anal. Appl., 13 (1992), pp. 138–144. [4] O. Axelsson and I. Kaporin, On the sublinear and superlinear rate of convergence of conjugate gradient methods, Numerical Algorithms, 25 (2000), pp. 1–22. [5] O. Axelsson and G. Lindskog, On the rate of convergence of the preconditioned conjugate gradient method, Numerische Mathematik, 48 (1986), pp. 499–523. [6] D. G. Barrett and B. Dherin, Implicit gradient regularization, arXiv preprint arXiv:2009.11162, (2020). [7] L. Beerens and D. J. Higham, Adversarial ink: Componentwise backward error attacks on deep learning, IMA Journal of Applied Mathematics, 89 (2024), pp. 175–196. [8] D. P. Bertsekas, Nonlinear Programming, Athena Scientific, 3rd ed., 2016, http://www. athenasc.com/nonlinbook.html. [9] P. Bürgisser and F. Cucker, Smoothed analysis of moore–penrose inversion, SIAM J. Matrix Anal. Appl., 31 (2010), pp. 2769–2783. [10] A. W. Chou, On the optimality of krylov information, Journal of Complexity, 3 (1987), pp. 26– 40. [11] T. A. Davis and Y. Hu, The University of Florida sparse matrix collection, ACM Trans. Math. Soft., 38 (2011), pp. 1–25. [12] M. Dereziński, E. N. Epperly, and R. A. Meyer, The matrix-vector complexity of Ax = b, arXiv preprint arXiv:2602.04842, (2026). [13] S. Di Giovacchino, D. J. Higham, and K. Zygalakis, Backward error analysis and the qualitative behaviour of stochastic optimization algorithms: Application to stochastic coordinate descent, arXiv preprint arXiv:2309.02082, (2023). [14] A. d’Aspremont, D. Scieur, and A. Taylor, Acceleration methods, Foundations and Trends in Optimization, 5 (2021), pp. 1–245. [15] E. N. Epperly, M. Meier, and Y. Nakatsukasa, Fast randomized least-squares solvers can be just as accurate and stable as classical direct solvers, Communications on Pure and Applied Mathematics, 79 (2026), pp. 293–339. [16] Y. Feng, T. Gao, L. Li, J.-G. Liu, and Y. Lu, Uniform-in-time weak error analysis for stochastic gradient descent algorithms via diffusion approximation, arXiv preprint arXiv:1902.00635, (2019). [17] K. V. Fernando and B. N. Parlett, Accurate singular values and differential qd-algorithms, Numer. Math., 67 (1994), pp. 191–229. [18] D. C.-L. Fong and M. Saunders, Lsmr: An iterative algorithm for sparse least-squares problems, SIAM Journal on Scientific Computing, 33 (2011), pp. 2950–2971. [19] G. H. Golub and C. F. Van Loan, Matrix Computations, The Johns Hopkins University Press, 1996. [20] G. H. Golub and C. F. Van Loan, Matrix computations, JHU press, 2013. [21] A. Greenbaum, Behavior of slightly perturbed lanczos and conjugate-gradient recurrences, Linear Algebra and its Applications, 113 (1989), pp. 7–63. [22] A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM, Philadelphia, PA, USA, 1997. [23] A. Greenbaum, V. Pták, and Z. e. k. Strakoš, Any nonincreasing convergence curve is possible for GMRES, SIAM J. Matrix Anal. Appl., 17 (1996), pp. 465–469. [24] M. R. Hestenes, E. Stiefel, et al., Methods of conjugate gradients for solving linear systems, Journal of research of the National Bureau of Standards, 49 (1952), pp. 409–436. [25] N. J. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, Philadelphia, PA, USA, second ed., 2002. [26] N. J. Higham and T. Mary, Mixed precision algorithms in numerical linear algebra, Acta Numerica, 31 (2022), pp. 347–414. [27] E. M. Kasenally, GMBACK: a generalised minimum backward error algorithm for nonsymmetric linear systems, SIAM J. Sci. Comput., 16 (1995), pp. 698–719. [28] E. M. Kasenally and V. Simoncini, Analysis of a minimum perturbation algorithm for nonsymmetric linear systems, SIAM J. Numer. Anal., 34 (1997), pp. 48–66. [29] J. Kuczyński and H. Woźniakowski, Estimating the largest eigenvalue by the power and lanczos algorithms with a random start, SIAM J. Matrix Anal. Appl., 13 (1992), pp. 1094– 1122.
24
M. DEREZIŃSKI, Y. NAKATSUKASA, AND E. REBROVA
[30] B. Laurent and P. Massart, Adaptive estimation of a quadratic functional by model selection, Ann. Stat., (2000), pp. 1302–1338. [31] C.-p. Lee and S. Wright, First-order algorithms converge faster than o(1/k) on convex problems, in International Conference on Machine Learning, PMLR, 2019, pp. 3754–3762. [32] A. Levin and E. Saff, Degree of approximation of real functions by reciprocals of real and complex polynomials, SIAM journal on mathematical analysis, 19 (1988), pp. 233–245. [33] C. Musco, C. Musco, and A. Sidford, Stability of the lanczos method for matrix function approximation, in Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, SIAM, 2018, pp. 1605–1624. [34] A. S. Nemirovsky and D. B. Yudin, Problem Complexity and Method Efficiency in Optimization, Wiley-Interscience, 1983, https://www2.isye.gatech.edu/∼nemirovs/Nemirovskii Yudin 1983.pdf. [35] Y. Nesterov, A method for solving the convex programming problem with convergence rate o (1/k2), in Dokl akad nauk Sssr, vol. 269, 1983, p. 543. [36] C. C. Paige, Error analysis of the lanczos algorithm for tridiagonalizing a symmetric matrix, IMA Journal of Applied Mathematics, 18 (1976), pp. 341–349. [37] C. C. Paige and M. A. Saunders, LSQR: An algorithm for sparse linear equations and sparse least squares, ACM Trans. Math. Soft., 8 (1982), pp. 43–71. [38] C. C. Paige and Z. Strakos, Residual and backward error bounds in minimum residual krylov subspace methods, SIAM J. Sci. Comput., 23 (2002), pp. 1898–1923. [39] G. Peters and J. H. Wilkinson, Inverse iteration, ill-conditioned equations and Newton’s method, SIAM Rev., 21 (1979), pp. 339–360. [40] L. F. Richardson, Ix. the approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam, Philosophical Transactions of the Royal Society of London. Series A, containing papers of a mathematical or physical character, 210 (1911), pp. 307–357. [41] Y. Saad, Iterative methods for sparse linear systems, SIAM, 2003. [42] A. Sankar, D. A. Spielman, and S.-H. Teng, Smoothed analysis of the condition numbers and growth factors of matrices, SIAM J. Matrix Anal. Appl., 28 (2006), pp. 446–476. [43] D. A. Spielman and S.-H. Teng, Smoothed analysis of algorithms: Why the simplex algorithm usually takes polynomial time, Journal of the ACM (JACM), 51 (2004), pp. 385–463. [44] L. N. Trefethen and D. Bau, Numerical Linear Algebra, SIAM, Philadelphia, 1997. [45] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, et al., Scipy 1.0: fundamental algorithms for scientific computing in python, Nature methods, 17 (2020), pp. 261–272. [46] J. von Neumann and H. Goldstine, Numerical inverting of matrices of high order, Bulletin of the American Mathematical Society, 53 (1947), pp. 1021–1099. [47] J. H. Wilkinson, Error analysis of floating-point computation, Numerische Mathematik, 2 (1960), pp. 319–340. [48] J. H. Wilkinson, Error analysis of direct methods of matrix inversion, Journal of the ACM (JACM), 8 (1961), pp. 281–330. [49] J. H. Wilkinson, Rounding errors in algebraic processes, Her Majesty’s Stationery Office, 1963. [50] J. H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University Press, 1965.
Appendix A. Analysis of inverse iteration: Proof of Lemma 4.8. Let λmin = λ1 ≤ λ2 · · · ≤ λn be the eigenvalues of M, and let ui be the corPn responding eigenvectors. Expressing v0 = i=1 ci ui , where by the assumption that Pn v0 is Gaussian, ci ∼ N (0, 1), we then have M−ℓ v0 = i=1 ci λ−ℓ i ui . Our goal is to find ℓ such that the Rayleigh quotient RM (vℓ ) := vℓ⊤ Mvℓ /∥vℓ ∥22 ≤ (1 + ∆)λ1 with probability 1 − δ, independent of the specific values of λi . We partition the index set I = {1, . . . , n} into two based on the target error ∆: Slow = {i : λi ≤ (1 + ∆/2)λ1 }
and
Shigh = {i : λi > (1 + ∆/2)λ1 }.
Note that index 1 is always in Slow , and the Rayleigh quotient is: ∥M−ℓ+1/2 v0 ∥22 RM (vℓ ) = = ∥M−ℓ v0 ∥22
P −(2ℓ−1) 2 −(2ℓ−1) + j∈Shigh c2j λj i∈Slow ci λi . P P 2 −2ℓ + 2 −2ℓ i∈Slow ci λi j∈Shigh cj λj
P
TOWARDS UNIVERSAL CONVERGENCE IN LINEAR SYSTEM SOLVERS
25
We seek to bound the relative error RM (vλℓ1)−λ1 . Since λi ≤ (1 + ∆/2)λ1 for all i ∈ Slow , the contribution of the first sum in the numerator is already “safe”, and the error is driven by Shigh . We simplify the quotient by considering only the c1 term in the denominator: P P −(2ℓ−1) (1 + ∆/2)λ1 i∈Slow c2i λ−2ℓ + j∈Shigh c2j λj i RM (vℓ ) ≤ P 2 −2ℓ i∈Slow ci λi P 2 −(2ℓ−1) j∈Shigh cj λj ≤ (1 + ∆/2)λ1 + . c21 λ−2ℓ 1 Focusing on the second term, for j ∈ Shigh , λj > (1 + ∆/2)λ1 . Thus: −(2ℓ−1)
λj
λ−2ℓ 1
= λ1
λ1 λj
2ℓ−1
< λ1
1 1 + ∆/2
2ℓ−1
Altogether, the relative error is bounded as follows: P 2ℓ−1 2 RM (vℓ ) − λ1 1 ∆ ∆ j∈Shigh cj + =: + Noise. ≤ 2 λ1 2 c1 1 + ∆/2 2 Now,Pusing properties of χ2 distributions we observe the following: 1. c2j follows χ2|Shigh | , so by [30, Eqn. (4.3)] we see that with probability 1−δ/2, P 2 we have cj ≤ n ln(2/δ). 2. c21 √ follows χ21 , so since the pdfpof a Gaussian random variable is bounded by 1/ 2π, we have P(c21 < η) ≤ 2η/π. To ensure that c21 ≥ η with probability 2 2 1 − δ/2, we need η ≤ πδ8 ; we thus set η = πδ8 . Combining these, with probability at least 1 − δ: P 2 16n ln(2/δ) j∈Shigh cj . (A.1) ≤ c21 πδ 2 For any j ∈ Shigh , we have λj > (1 + ∆ 2 )λ1 . The noise term is bounded by: 2 j∈Shigh cj c21
P (A.2)
Noise ≤
!
λ1 λj
2ℓ−1
≤
16n ln(2/δ) πδ 2
1 1 + ∆/2
2ℓ−1
To achieve RM (vℓ ) ≤ (1 + ∆)λ1 , we require the noise term to be ≤ ∆/2: (A.3)
1 1 + ∆/2
2ℓ−1 ≤
πδ 2 ∆ 32n ln(2/δ)
Taking the logarithm and solving for ℓ completes the proof: 32n ln(2/δ) ln 2 πδ ∆ 1 . (A.4) ℓ ≥ 1 + 2 ln(1 + ∆/2) In particular, for a 1/2-relative error (∆ = 0.5), we require ℓ ≥ 2.23 ln(n/δ 2 ).
26
M. DEREZIŃSKI, Y. NAKATSUKASA, AND E. REBROVA
Appendix B. Proof of Proposition 4.11. We start by examining (4.5), in which we note that [b, AQk ] = Uk [∥b∥2 e, Bk ] =: Uk N, where N ∈ R(k+1)×(k+1) . We then see that (4.5) reduces to ⊤ −1 0 0 −1 (B.1) N⊤ N =λ . y Ik Ik y Now N is upper bidiagonal; we can apply a congruence to zero-out transformation w1 (k+1)×(k+1) the (1, 2) entry, that is, defining the matrix W = ∈ R where 0 Ik N (Bk )1,1 N1,2 , 0, . . . , 0] (where we note that N1,2 = ∥b∥ ), that is, W is equal to w1 = [1, − N1,1 1,1 2 Ik+1 except in the (1,2) entry. Then consider ⊤ 0 ⊤ ⊤ −1 −1 ⊤ 0 −1 −1 W N NW W = λW W W , y Ik Ik y which simplifies to ⊤ ∥b∥2 ∥b∥2 0 −1 −1 (B.2) W =λ y B̃k B̃k
⊤ 0 Ik
−1 −1 W Ik y
where B̃k is the submatrix of B obtained by removing the first row. As in the PSD case, we recognize (B.2) as a block diagonal generalized eigenvalue problem (whose smallest eigenpair is desired) with the first 1 × 1 block having eigenvalue at infinity. The second block has Ik in the right-hand side, so the problem is equivalent to finding the smallest singular value and its corresponding right singular b ∥2 /∥b vector of B̃k , which is bidiagonal. The proof of berrA,b (b x) = ∥B̃k v v∥2 is as in Lemma 4.6, noting that (b x)1 = −1. Appendix C. Polynomial approximation. We use the following simple fact. Lemma C.1. If γ ∈ [0, π/2] then for any integer a ≥ 1, sin(aγ) ≤ a sin(γ). Proof. This holds by induction on a. When a = 1, the result trivially holds, assume that it holds for some a > 1. Then sin((a + 1)γ) = sin(aγ) cos(γ) + cos(aγ) sin(γ) ≤ a sin(γ) · 1 + sin(γ) · 1 = (a + 1) sin(γ). Proof of Lemma 4.4. Recall that G(x) =
1 − Tℓ∗ (x) 1 − cos(ℓ arccos(2x − 1)) = . 2ℓ2 2ℓ2
For x = sin2 (γ) we have 2x − 1 = − cos(2γ) and arccos(− cos(2γ)) = π − 2γ. Since ℓ is even by construction, cos(ℓ(π − 2γ)) = cos(2ℓγ) = 1 − 2 sin2 (ℓγ). So, G(x) =
1 − 1 + 2 sin2 (ℓγ) sin2 (ℓγ) = . 2ℓ2 ℓ2
Using Lemma C.1, we have then G(x) ≤
ℓ2 sin2 (γ) = x, ℓ2
TOWARDS UNIVERSAL CONVERGENCE IN LINEAR SYSTEM SOLVERS
27
and G(x) ≥ 0 as Tℓ∗ ∈ [−1, 1]. Finally, zero is the root of f (x) := x−G(x) multiplicity at least two because f (0) = −G(0) =
1 − cos(ℓ arccos(−1)) = 0; ℓ2
f ′ (0) = 1 − G′ (0) = 1 −
2ℓ2 (−1)ℓ =0 2ℓ2
using that ℓ is even and the standard derivative properties of Chebyshev polynomials. Proof of Lemma 4.5. First, we show that Fℓ (γ) cannot be maximized at γ ∈ π π , 2 ]. Indeed, let ( 2ℓ h πi 1 v := arcsin sin(ℓγ) ∈ 0, . ℓ 2ℓ Then, sin2 (ℓv) = sin2 (ℓγ) but v ≤ γ, so sin2 v ≤ sin2 γ since sin(·) is increasing on xy [0, π2 ]. Also, for y ∈ [0, 1], the function ϕy (x) = x−y is decreasing, so Fℓ (γ) ≤ Fℓ (v). π Now, we show that Fℓ (γ) is decreasing on γ ∈ [0, 2ℓ ]. To simplify the computation, we show that 1/Fℓ has nonnegative derivative on this interval. Define
Gℓ (u) := G′ℓ (u) =
1 1 = csc2 (ℓu) − 2 csc2 u, Fℓ (u) ℓ
2 h(u) − h(ℓu) , 2 3 ℓ u
h(x) := x3 csc2 x cot x =
x3 cos x . sin3 x
So it suffices to know that h is decreasing on (0, π/2]. Consider h′ (x) =
x2 2 3 sin x cos x − x(1 + 2 cos x) . sin4 x
Let p(x) := 3 sin x cos x − x(1 + 2 cos2 x). Then p′ (x) = 4 sin x (x cos x − sin x) ≤ 0, because x 7→ sin x/x decreases on (0, π/2], or equivalently 0≥
sin x ′ x
=
x cos x − sin x . x2
Since p(0) = 0, we get p(x) ≤ 0, hence h′ (x) ≤ 0, so h is decreasing. Therefore, for 0 < u ≤ π/(2ℓ), since ℓu ≤ π/2, we have h(u) ≥ h(ℓu), which implies that G′ℓ (u) ≥ 0. Therefore, the supremum of Fℓ is the limit as u ↓ 0. Using Taylor expansion for cosecant, we have 1 1 1 1 ℓ2 − 1 2 2 Gℓ (u) = + O(u ) − + + + O(u ) = + O(u2 ). ℓ2 u2 3 ℓ2 u2 3ℓ2 3ℓ2 Hence sup Fℓ (u) = lim Fℓ (u) = u>0
u↓0
3ℓ2 . ℓ2 − 1