Solving Dense Linear Systems Faster Than via Preconditioning
A DPP-analyzed stochastic block solver achieves ˜O((n²+nk^{ω−1})log(1/ε)) for dense systems.
Key Findings
Methodology
The paper develops a randomized block coordinate descent framework combining Kaczmarz, coordinate descent, and sketch-and-project ideas. k-DPP sampling and elementary symmetric polynomials provide the ideal convergence analysis; a randomized Hadamard transform and Markov-chain coupling replace expensive DPP sampling; sketch-preconditioned conjugate gradient approximately solves each block subproblem.
Key Results
- For dense n×n systems, the algorithm returns ‖Aẋ−b‖≤ε‖b‖ in ˜O((n²+nk^{ω−1})log(1/ε)), where k counts singular values exceeding a constant multiple of the smallest positive singular value.
- With ω<2.372, k=O(n^0.729) yields ˜O(n²), improving over direct O(n^ω) solution and avoiding the ˜O(n²k^{ω−2}) power-iteration cost often associated with low-rank preconditioning.
- Extensions give ˜O((nnz(A)+d²+dk^{ω−1})log(1/ε)) for least squares and ˜O((nnz*(A)+nk^{ω−1})log(1/ε)) for PSD systems. No real datasets or numerical benchmark tables are reported.
Significance
The work turns a flat spectral tail—previously useful mainly for analyzing CG—into a direct algorithmic resource. It targets regularized, noisy, covariance, and Hessian-related systems in which only a few singular values are large and most are comparable. Instead of paying upfront for an accurate global low-rank preconditioner, the method distributes spectral exploitation across randomized iterations, retaining near-input-scale n² complexity for a substantially wider range of k.
Technical Contribution
The technical contribution has three layers. First, Schur-concavity of ratios of elementary symmetric polynomials, combined with majorization, sharpens k-DPP block-coordinate convergence: blocks of size ˜O(k) require only ˜O(n/k) iterations. Second, Markov-chain coupling shows that cheaper sampling after randomized Hadamard preprocessing preserves the guarantee up to logarithmic factors. Third, sketch-preconditioned PCG solves each highly underdetermined block approximately in O(log(n/k)) inner iterations, reducing an update to ˜O(nk+k^ω).
Novelty
The novelty is not any isolated use of DPPs, PCG, or sketching. It is their end-to-end integration into a dense solver that exploits spectral tails without explicitly constructing a low-rank preconditioner. Relative to power-iteration or low-rank-SVD approaches, the algorithm avoids globally approximating all leading singular directions while retaining a k-sensitive complexity guarantee.
Limitations
- The theory is stated in the exact-arithmetic Real-RAM model. The authors argue that sparse subspace embeddings, SVD, and PCG should support numerical stability, but a complete finite-precision theorem and empirical error study are absent.
- The paper reports complexity theorems rather than datasets, wall-clock measurements, or numerical comparisons with CG and preconditioned CG; practical constants, sampling overhead, and memory behavior therefore remain uncertain.
- If k is close to n, the nk^{ω−1} term approaches direct-solver cost; if the spectral tail is not flat, the central structural advantage largely disappears.
Future Work
Important next steps include Word-RAM and floating-point stability proofs, large-scale synthetic and application benchmarks, and optimized implementations of Hadamard preprocessing, sampling, and inner PCG. Further directions include adaptive estimation of k, rectangular matrix-multiplication exponents, distributed and GPU implementations, dynamic systems, and extensions to nonsymmetric operators.
AI Executive Summary
Dense linear systems Ax=b are classically solved in O(n^ω) time, while conjugate gradient can still require roughly ˜O(n²k) work when only k singular values are large. Low-rank preconditioning reduces iterations but often pays a substantial power-iteration or spectral-approximation cost first. Dereziński and Yang ask whether a flat spectral tail can be exploited without building a global preconditioner.
Their answer is a randomized block coordinate solver. The analysis uses k-DPP sampling, majorization, and elementary symmetric polynomials to show that ˜O(k)-sized blocks need only ˜O(n/k) outer updates. A randomized Hadamard transform and Markov-chain coupling replace expensive DPP sampling with a cheaper procedure. Matrix sketching and preconditioned conjugate gradient approximately solve each block subproblem, giving ˜O(nk+k^ω) work per update and overall ˜O((n²+nk^{ω−1})log(1/ε)).
Thus k=O(n^0.729) yields ˜O(n²), while extensions cover sparse PSD systems and least squares. The result applies to regularized systems such as (K+λI)x=b, yet the paper supplies theoretical bounds rather than dataset experiments. Open issues include floating-point stability, engineering constants, parallel implementation, and performance when k is large. The broader contribution is a principled alternative to explicit preconditioning for flat-tailed spectra.
Deep Analysis
Background
Dense systems have long been addressed by elimination, Cholesky, or fast matrix multiplication at O(n^ω). CG is effective for favorable conditioning, but systems with k large singular values can still incur ˜O(n²k). Regularization A=B+λI and noise A=B+δN often create flat tails. Existing preconditioners based on power iteration or low-rank SVD exploit this structure but may cost nearly as much as direct solution.
Core Problem
Given A∈R^{n×n}, compute ẋ with ‖Aẋ−b‖≤ε‖b‖ while scaling with the number k of dominant singular values rather than the full dimension or condition number. The difficulty is simultaneous: sampling must be spectrally informative and cheap; block updates solve underdetermined systems; and sampling and inexact-solve errors must preserve global convergence.
Innovation
- ��Use k-DPPs and elementary symmetric polynomials for a sharp spectral convergence bound.
- ��Apply majorization and Schur-concavity of ESP ratios to obtain a k-factor improvement.
- ��Use Markov-chain coupling to transfer the guarantee to cheaper sampling after randomized Hadamard preprocessing.
- ��Use sketch-preconditioned PCG for approximate block updates, eliminating an explicit global low-rank preconditioner.
Methodology
- ��Input: A, b, and ε; characterize the tail through the averaged condition number of A−Ak.
- ��Sampling: analyze ideal k-DPP subsets, whose probabilities are determinant-proportional; implement a cheaper post-Hadamard sampler.
- ��Outer update: form a small row or column block and perform a sketch-and-project-style correction.
- ��Inner solve: construct a sparse subspace embedding and run PCG on the block subproblem for O(log(n/k)) iterations.
- ��Convergence: ˜O(n/k) outer updates, each costing ˜O(nk+k^ω), plus preprocessing, producing the main theorem.
Experiments
This is primarily a theory and algorithm paper. It reports no ImageNet, LIBSVM, or other application datasets, and no train/test protocol. The implicit baselines are direct solvers, CG, low-rank-preconditioned CG, and sketching methods. The evaluated quantities are asymptotic runtime, residual accuracy, and dependence on k and the matrix-multiplication exponent ω, rather than measured speedups.
Results
The dense guarantee is ˜O((n²+nk^{ω−1})log(1/ε)); because 1/(ω−1)≈0.729, k=O(n^0.729) gives ˜O(n²). Least squares requires ˜O((nnz(A)+d²+dk^{ω−1})log(1/ε)), while PSD systems require ˜O((nnz*(A)+nk^{ω−1})log(1/ε)). These are proven complexity bounds, not experimentally measured percentage improvements.
Applications
Relevant settings include kernel ridge regression, covariance and Hessian-related PSD systems, l2-regularized least squares, and scientific problems whose spectra are flattened by noise or compression. Sparse PSD instances can exploit nnz*(A)=n·max_i nnz(A_i,:). Practical deployment requires fast matrix operations, sketching, stable PCG, and residual-based stopping.
Limitations & Outlook
The exact-arithmetic assumption leaves finite-precision behavior incompletely justified. The absence of datasets and wall-clock studies leaves constants, parallelism, and memory traffic unknown. When k is large or the tail is not flat, nk^{ω−1} can dominate and the advantage over direct methods vanishes. Future work should address floating-point proofs, real benchmarks, GPU/distributed execution, and adaptive block selection.
Plain Language Accessible to non-experts
Imagine a factory trying to coordinate thousands of machines so that every order is completed correctly. The slow traditional strategy is to build one enormous master schedule before production starts. That resembles directly solving the whole system or constructing a detailed global preconditioner; the preparation can be almost as expensive as the task itself.
This paper instead checks a small group of machines at a time. The ideal lottery, called DPP sampling, prefers groups whose information is different rather than repetitive. The authors prove that these diverse groups reduce the factory’s mistakes quickly, then design a cheaper lottery that behaves similarly. Each local correction is also allowed to be approximate, provided it is accurate enough to keep the whole process moving toward the target.
The method is especially useful when only a few machines are unusually difficult while most behave similarly. In that case, the expensive directions are handled through blocks and the easy tail is processed cheaply. The promised cost is ˜O((n²+nk^{ω−1})log(1/ε)); when the number k of difficult directions is at most about n^0.729, it is near n². The paper proves this mathematically but does not yet demonstrate the advantage on real factory-like datasets.
ELI14 Explained like you're 14
Suppose you are solving a giant school timetable puzzle. Checking every possible class and room at once would take forever. A smarter plan is to choose a small group of classes, fix the biggest clashes, and repeat. But choosing random classes can be bad: you might keep checking the same kind of clash and learn almost nothing.
The researchers study a clever way to pick groups that are different from one another. Their ideal picker is called a DPP, but using it directly is too slow, so they create a cheaper picker that behaves almost the same after a random shuffle. Then, for each group, they solve a smaller puzzle instead of the entire timetable. They do not need the smaller puzzle to be perfectly solved; “good enough” is enough if every round improves the timetable.
The cool part is that the speed depends on how many directions are genuinely difficult. Call that number k. If k is not too large—roughly no more than n^0.729—the total work is close to n², about the cost of even reading a dense problem. This can help when noise or regularization makes most directions similarly easy.
There is an important catch: this is mainly a mathematical result, not a flashy app demonstration. The paper does not give real datasets or stopwatch comparisons. It also assumes exact arithmetic in its main proofs. So the next challenge is to test the idea on real systems, computers, and noisy data!
Glossary
Singular value
A singular value measures how strongly a matrix stretches a particular direction. The collection of singular values controls conditioning and reveals whether most directions form a flat tail.
The paper defines k by counting singular values much larger than the smallest positive one.
k-DPP
A determinantal point process conditioned to select exactly k items, with subset probabilities proportional to determinants. It favors diverse, nonredundant selections.
It is the ideal block-sampling distribution used in the convergence analysis.
Majorization
A partial order describing whether one vector is more evenly distributed than another. With Schur-concavity, it yields bounds for symmetric functions.
It sharpens the elementary-symmetric-polynomial convergence analysis.
Matrix sketching
A randomized compression that approximately preserves the geometry needed by a computation. It reduces the cost of constructing and solving auxiliary problems.
It creates fast approximate block updates and preconditioners.
Preconditioned conjugate gradient
CG enhanced by an easily solvable transformation that improves convergence. Here it solves each underdetermined block only approximately.
The inner solver needs O(log(n/k)) iterations.
Flat-tailed spectrum
A spectrum with a few large leading values and many tail values of comparable size. Such structure separates a small difficult part from an easy bulk.
It is the structural reason the proposed runtime can depend on k.
Open Questions Unanswered questions from this research
- 1 A complete finite-precision theorem is still missing. It must show that sampling, sketching, and PCG errors preserve the claimed convergence on real hardware.
- 2 There are no application datasets or wall-clock studies, so practical constants, communication costs, cache behavior, and comparison with optimized CG remain unknown.
- 3 The advantage weakens when k is near n or the tail is not flat; adaptive online estimation of spectral structure and block size remains open.
Applications
Immediate Applications
Kernel ridge regression and regularized systems
For (K+λI)x=b and related PSD systems, the solver may reduce work when only a few eigenvalues lie well above the regularization scale. Users need matrix products or sparse access, a stable PCG implementation, and residual-based accuracy control.
Noisy covariance and scientific models
Measurement noise, compression, or privacy perturbations can flatten spectral tails. Teams could apply the method to covariance, Hessian, or implicit dense operators, provided randomized transforms, sketches, and reliable stopping tests are available.
Long-term Vision
Parallel flat-spectrum linear algebra
GPU or distributed implementations could process sampled blocks concurrently and approach input-scale cost on favorable spectra. Major obstacles are synchronization, communication, finite-precision stability, and efficient rectangular matrix multiplication.
Abstract
We give a stochastic optimization algorithm that solves a dense $n\times n$ real-valued linear system $Ax=b$, returning $\tilde x$ such that $\|A\tilde x-b\|\leq ε\|b\|$ in time: $$\tilde O((n^2+nk^{ω-1})\log1/ε),$$ where $k$ is the number of singular values of $A$ larger than $O(1)$ times its smallest positive singular value, $ω< 2.372$ is the matrix multiplication exponent, and $\tilde O$ hides a poly-logarithmic in $n$ factor. When $k=O(n^{1-θ})$ (namely, $A$ has a flat-tailed spectrum, e.g., due to noisy data or regularization), this improves on both the cost of solving the system directly, as well as on the cost of preconditioning an iterative method such as conjugate gradient. In particular, our algorithm has an $\tilde O(n^2)$ runtime when $k=O(n^{0.729})$. We further adapt this result to sparse positive semidefinite matrices and least squares regression. Our main algorithm can be viewed as a randomized block coordinate descent method, where the key challenge is simultaneously ensuring good convergence and fast per-iteration time. In our analysis, we use theory of majorization for elementary symmetric polynomials to establish a sharp convergence guarantee when coordinate blocks are sampled using a determinantal point process. We then use a Markov chain coupling argument to show that similar convergence can be attained with a cheaper sampling scheme, and accelerate the block coordinate descent update via matrix sketching.