math.NAApr 17, 2026

Towards Universal Convergence of Backward Error in Linear System Solvers

Authors: Michał DerezińskiYuji NakatsukasaElizaveta Rebrova

Organizations: University of Michigan, Ann Arbor, MI · University of Oxford, Oxford, UK · Princeton University, Princeton, NJ

Abstract

The quest for an algorithm that solves an n×nn\times n linear system in O(n2)O(n^2) time complexity, or O(n2poly(1/ε))O(n^2 \text{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/k1/k (relative) backward error after kk iterations on any positive semidefinite (PSD) linear system, irrespective of its condition number. This universal convergence rate implies an O(n2/ε)O(n^2/ε) 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 backward error over a Krylov subspace, we attain an even faster O(1/k2)O(1/k^2) universal rate, and we turn this into an efficient algorithm, MINBERR, with complexity O(n2/ε)O(n^2/\sqrtε). Finally, we extend this approach via normal equations to solving general linear systems in O(n2log(n)/ε)O(n^2\log(n)/ε) time complexity. We report strong numerical performance of our algorithms on benchmark problems.

Explore similar work

Apr 25, 2026cs.DS

Well-Conditioned Oblivious Perturbations in Linear Space

Perturbing a deterministic nn-dimensional matrix with small Gaussian noise is a cornerstone of smoothed analysis of algorithms [Spielman and Teng, JACM 2004], as it reduces the condition number of the input to O(n)O(n), and with it the complexity of many matrix algorithms. However, when deployed algorithmically, these perturbations are expensive due to the cost of generating and storing n2n^2 Gaussian random variables. We propose a perturbation that requires generating and storing O(n)O(n) random numbers in O(logn)O(\log n) bits of precision, and reduces the condition number of any deterministic matrix to O(n)O(n), matching Gaussian perturbations. Our result in particular implies a better complexity for the perturbed conjugate gradient algorithm, showing that we can solve an n×nn\times n linear system in linear space to within an arbitrarily small constant backward error using O(n)O(n) matrix-vector products. In our construction, we introduce the concept of a pattern matrix, which is a dense deterministic matrix that maps all sparse vectors into dense vectors, and we combine it with a sparse perturbation whose entries are dependent and located in a non-uniform fashion. In order to analyze this construction, we develop new techniques for lower bounding the smallest singular value of a random matrix with dependent entries.
Shabarish Chenakkod, Michał Dereziński, Xiaoyu Dong +1
Jul 23, 2026math.OC

Barzilai-Borwein Fails Superlinear Convergence on an Open Set of Quadratics for Every Dimension n4n\geq 4

Barzilai--Borwein (BB) method has shown strong practical performance in continuous optimization, yet its convergence dynamics remains poorly understood. In particular, a central unresolved question is whether BB converges superlinearly for almost every strictly convex quadratic problem and initialization. We provide a negative answer to this question. Specifically, for every finite dimension n4n\geq4, we construct a nonempty open, hence positive-Lebesgue-measure, family of strictly convex quadratic problems and initial points for which the long Barzilai--Borwein method (BB1) converges but cannot converge root-superlinearly. More precisely, with the explicit constants ρmin=106,ρmax=0.61ρ_{\min}=10^{-6},ρ_{\max}=0.61, every spectral component of the gradient is bounded above and below by the corresponding geometric sequence. Consequently, the gradient norm and the energy norm of the error satisfy two-sided geometric estimates with the same rates, while the objective gap satisfies the corresponding estimates with squared rates. In particular, all three quantities are bounded below by geometric sequences, ruling out superlinear convergence. The construction is highly nontrivial, based on a computer-assisted proof of a nonresonant, attracting seven-cycle of the projectivized BB dynamics in dimension four.
Dawei Li, Xiaotian Jiang, Mingyi Hong
Aug 12, 2026math.OC

Tight Nonasymptotic Local Convergence of Sinkhorn-Knopp

We revisit the Sinkhorn-Knopp (SK) algorithm for the matrix scaling problem. Despite extensive literature on the global convergence of SK and its variants, its local linear convergence behavior remains less understood. We address this gap by providing the first nonasymptotic local analysis of SK that matches the rate obtained from existing asymptotic Jacobian-based arguments. We show that under certain connectivity conditions, SK is a polynomial-time algorithm for doubly stochastic matrix scaling. With the developed tools, we showcase the local suboptimality of SK and provide accelerated variants. Finally, for dense matrices, we improve the complexity of existing first-order matrix scaling algorithms from O(n7/3ε2/3)O(\tfrac{n^{7/3}}{\varepsilon^{2/3}}) to O(n9/4ε)O(\tfrac{n^{9/4}}{\sqrt{\varepsilon}}).
Wenzhi Gao, Zhaonan Qu, Yinyu Ye +1