Convergence of AdaRBFGS: adaptive randomized BFGS matrix inversion with iterate-dependent sketching

The sketch-and-project framework of Gower and Richtárik yields randomized quasi-Newton methods that invert a symmetric positive definite matrix by repeatedly projecting the current estimate onto randomly sketched inverse equations; when the sketching matrix is sampled i.i.d. from a fixed distribution, these methods converge globally and linearly in expectation. The adaptive variant AdaRBFGS instead builds each sketch from the current iterate itself ($S=L_k\tilde S$ where $X_k=L_kL_k^\top$), so the sampling distribution co-evolves with the approximation and, if the method converges, asymptotically approaches the rate-optimal distribution; numerically it dramatically outperforms non-adaptive randomized BFGS (Broyden–Fletcher–Goldfarb–Shanno) and classical fixed-point iterations for inversion such as Newton–Schulz, apparently converging superlinearly. Yet no convergence proof is known — only a conditional one-step bound — and the source work states it as an explicit open question whether AdaRBFGS converges linearly or superlinearly. The problem is to establish a global convergence theory (an explicit linear or superlinear rate) for AdaRBFGS and related adaptive randomized matrix-inversion methods under implementable assumptions, or to refute such guarantees by an explicit counterexample. Later literature has settled the analogous adaptive-sampling question for the linear-systems (vector) version of sketch-and-project but not for matrix inversion.

Background

Randomized iterative methods for linear systems — such as the randomized Kaczmarz method, which projects the current iterate onto a randomly chosen constraint hyperplane $\{x: a_i^\top x = b_i\}$, and randomized coordinate descent — were unified by Gower and Richtárik into a single sketch-and-project framework (Gower and Richtárik, SIAM J. Matrix Anal. Appl. 36(4), 2015). The idea is to replace a large deterministic problem by a randomly sketched simpler one that retains all its solutions: to solve $Ax=b$, draw a random sketching matrix $S$ and project the current iterate, in a weighted inner product, onto the solution set of the sketched system $S^\top A x = S^\top b$. The same principle applies to matrix inversion: to compute $A^{-1}$, one iteratively projects a matrix estimate $X_k$ onto the solution set of the sketched inverse equations $S^\top A X = S^\top$, i.e. a random subset of the equations $AX=I$. Since $S$ typically has only $q \ll n$ columns, each step requires only a $q\times q$ factorization, making the method cheap for large $n$.

For symmetric positive definite $A$, choosing the projection weight $B=A$ and enforcing symmetry of the iterates yields the randomized block BFGS (RBFGS) update — a randomized block variant of the classical Broyden–Fletcher–Goldfarb–Shanno (BFGS) quasi-Newton update:

$X_{k+1} = S(S^\top AS)^{-1}S^\top + \left(I - S(S^\top AS)^{-1}S^\top A\right)X_k\left(I - AS(S^\top AS)^{-1}S^\top\right).$

If $X_0$ is symmetric positive definite, so are all iterates, which is why such methods are attractive for building preconditioners (approximate inverses used to accelerate the iterative solution of linear systems) and variable-metric (quasi-Newton) optimization methods. When $S$ is sampled i.i.d. from a fixed distribution $\mathcal{D}$ at every iteration, Gower and Richtárik proved global linear convergence in expectation of the weighted Frobenius error, where $\|M\|^2_{F(B)} = \mathrm{Tr}(M^\top B M B)$:

$\mathbb{E}\left[\|X_k - A^{-1}\|^2_{F(B)}\right] \le \rho^k \|X_0 - A^{-1}\|^2_{F(B)}, \qquad \rho = 1 - \lambda_{\min}\left(B^{-1/2}\mathbb{E}[Z]\,B^{-1/2}\right),$

with $Z = A^\top S(S^\top AB^{-1}A^\top S)^{-1}S^\top A$ and $\lambda_{\min}$ the smallest eigenvalue; for the RBFGS specialization $B=A$ this is convergence in $\|\cdot\|_{F(A)}$ with $\|M\|^2_{F(A)} = \mathrm{Tr}(M^\top A M A)$. In the complete discrete sampling case — a fixed list of sketch blocks $S_1,\dots,S_r$ with $S=[S_1,\dots,S_r]$ and sampling probabilities $p_i=\mathrm{Tr}(S_i^\top AS_i)/\mathrm{Tr}(S^\top AS)$ — this rate becomes $\rho = 1-\lambda_{\min}(S^\top AS)/\mathrm{Tr}(S^\top AS)$, which is minimized when $S^\top AS \propto I$, i.e. when $S$ is a column selection from $A^{-1/2}$ — precisely the matrix being computed.

This observation motivated the adaptive method AdaRBFGS (Gower's thesis, 2016; Gower and Richtárik, SIAM J. Matrix Anal. Appl. 38(4), 2017). Maintain the iterate in factored form $X_k = L_kL_k^\top$ and at iteration $k$ use the sketch $S = L_k\tilde S$, where $\tilde S$ is drawn i.i.d. from a fixed distribution $\tilde{\mathcal{D}}$ (e.g. a random block of $q$ coordinate vectors, or a Gaussian matrix). If the method succeeds then $L_k \to A^{-1/2}$, so the sketch distribution automatically approaches the rate-optimal one: computing the inverse and optimizing the sampling distribution are, in the authors' words, in synchrony. Numerically, AdaRBFGS is dramatically faster than non-adaptive randomized BFGS and than the Newton–Schulz iteration (the classical fixed-point recursion $X_{k+1}=X_k(2I-AX_k)$ for the inverse) and the self-conditioned minimal-residual method (a minimal-residual inversion scheme that reuses its own iterates as preconditioners) on large positive definite matrices, with apparently superlinear residual decay.

However, the convergence proofs break down: because the distribution of $S$ now depends on the random iterate $X_k$, the fixed $\mathbb{E}[Z]$ contraction argument does not apply. The only known result is a conditional one-step bound (Theorem 9.1 of the SIMAX paper; Theorem 41 of the thesis):

$\mathbb{E}\left[\|X_{k+1}-A^{-1}\|^2_{F(A)} \mid X_k\right] \le \left(1 - \frac{\lambda_{\min}(AX_k)}{\mathrm{Tr}(AX_k)}\right)\|X_k - A^{-1}\|^2_{F(A)}.$

This bound cannot be telescoped into a global rate: the factor $1-\lambda_{\min}(AX_k)/\mathrm{Tr}(AX_k)$ approaches $1$ whenever $AX_k$ is ill-conditioned, and no a priori control of the conditioning of $AX_k$ along the random trajectory is available. It only shows that if $X_k \to A^{-1}$, then the asymptotic one-step factor tends to the dimension-optimal value $1-1/n$ — consistent with, but not proving, the observed superlinear behavior. The thesis states: "we have still to prove that the AdaRBFGS methods are even linearly convergent. Thus it is now an open question whether or not the AdaRBFGS convergent linearly or superlinearly", and the published paper leaves "a deeper analysis of this and other adaptive variants" to future work.

For context on what has since been settled: for the linear-systems (vector) version of sketch-and-project, Gower, Molitor, Moorman and Needell (SIAM J. Matrix Anal. Appl. 42(2), 2021) proved global linear convergence of adaptive sampling rules whose probabilities depend on the current iterate (max-distance, proportional, and capped rules based on sketched residuals), and adaptive sketch-and-project methods with proven linear convergence have been extended to tensor equations (Tang, Zhang and Li, 2026). None of these analyses covers the matrix-inversion setting of AdaRBFGS, where the sketch is built from a factor of the iterate rather than from residual magnitudes; the iterate-dependent sketching of AdaRBFGS remains without a convergence theory.

Problem Statement

Let $A\in\mathbb{R}^{n\times n}$ be symmetric positive definite and let $X_0 = L_0L_0^\top$ be symmetric positive definite with $L_0$ invertible. Fix a distribution $\tilde{\mathcal{D}}$ over $n$-row matrices of full column rank (for instance: $\tilde S = I_{:C}$, a uniformly chosen block of $q$ coordinate vectors, where $I_{:C}$ denotes the $n\times q$ matrix formed by the columns of the identity matrix indexed by the set $C$; or a standard Gaussian $n\times q$ matrix). A closely related adaptive variant is the complete discrete sampling case: blocks $C_1,\dots,C_r$ partitioning $\{1,\dots,n\}$ chosen with probabilities $p_i=\mathrm{Tr}(S_i^\top AS_i)/\mathrm{Tr}(S^\top AS)$, where $S_i=L_kI_{:C_i}$ and $S=[S_1,\dots,S_r]=L_k$, so that both the blocks' probabilities and the sketches themselves depend on the iteration index; since these probabilities depend on the current factor $L_k$, this scheme no longer fits the fixed-$\tilde{\mathcal{D}}$ framework and is included in the target as an adaptive variant in its own right. The AdaRBFGS method (Gower and Richtárik, SIAM J. Matrix Anal. Appl. 38(4), 2017, Algorithm 4) maintains $X_k = L_kL_k^\top$ and at each iteration draws an independent copy $\tilde S_k\sim\tilde{\mathcal{D}}$, forms the iterate-dependent sketch $S = L_k\tilde S_k$, and performs the randomized block BFGS update

$X_{k+1} = S(S^\top AS)^{-1}S^\top + \left(I - S(S^\top AS)^{-1}S^\top A\right)X_k\left(I - AS(S^\top AS)^{-1}S^\top\right),$

updating the factor $L_k$ so that $X_{k+1}=L_{k+1}L_{k+1}^\top$. For the non-adaptive case ($S$ i.i.d. from a fixed distribution independent of the iterates) global linear convergence in $\mathbb{E}\|X_k-A^{-1}\|^2_{F(A)}$ is known, and for AdaRBFGS the only known guarantee is the conditional one-step bound

$\mathbb{E}\left[\|X_{k+1}-A^{-1}\|^2_{F(A)}\mid X_k\right] \le \left(1-\frac{\lambda_{\min}(AX_k)}{\mathrm{Tr}(AX_k)}\right)\|X_k - A^{-1}\|^2_{F(A)},$

which cannot be telescoped because the conditioning of $AX_k$ is uncontrolled along the trajectory.

Question: establish or refute the global convergence of the AdaRBFGS methods — the general fixed-$\tilde{\mathcal{D}}$ scheme above and its complete-discrete-sampling adaptive variant.

(1) Positive resolution — prove that $\mathbb{E}\|X_k - A^{-1}\|^2_{F(A)} \to 0$ (convergence in $\|\cdot\|_{F(A)}$ is equivalent to convergence in the standard Frobenius norm because $A$ is fixed; almost-sure convergence of $X_k$ to $A^{-1}$ is also acceptable) for all symmetric positive definite $A$ and symmetric positive definite $X_0$, or under precisely stated additional assumptions that are implementable, i.e. verifiable without knowledge of $A^{-1}$ or of other quantities the method is trying to compute; the proof must yield an explicit nonasymptotic convergence rate (a linear rate, or the superlinear rate the numerical evidence suggests), and must genuinely handle the defining feature that the sampling distribution of $S$ depends on the random iterate $X_k$.

(2) Negative resolution — exhibit an explicit symmetric positive definite $A$, an explicit initial factor $L_0$ with $X_0 = L_0L_0^\top$ positive definite, and either an explicit sketch distribution $\tilde{\mathcal{D}}$ (for the general fixed-$\tilde{\mathcal{D}}$ scheme) or an explicit partition $\{C_1,\dots,C_r\}$ of $\{1,\dots,n\}$ (for the complete-discrete-sampling adaptive variant), for which the corresponding AdaRBFGS iterates provably fail to converge in expectation to $A^{-1}$ (or provably admit no linear convergence rate), together with a rigorous demonstration of that failure.

References

1. Robert M. Gower, Sketch and Project: Randomized Iterative Methods for Linear Systems and Inverting Matrices, PhD thesis, University of Edinburgh, 2016. arXiv:1612.06013, https://arxiv.org/abs/1612.06013
2. Robert M. Gower and Peter Richtárik, Randomized Quasi-Newton Updates Are Linearly Convergent Matrix Inversion Algorithms, SIAM Journal on Matrix Analysis and Applications 38(4):1380–1409, 2017. DOI 10.1137/16M1062053, arXiv:1602.01768, https://arxiv.org/abs/1602.01768
3. Robert M. Gower and Peter Richtárik, Randomized Iterative Methods for Linear Systems, SIAM Journal on Matrix Analysis and Applications 36(4):1660–1690, 2015. DOI 10.1137/15M1025487, arXiv:1506.03296, https://arxiv.org/abs/1506.03296
4. Robert M. Gower, Denali Molitor, Jacob Moorman, and Deanna Needell, Adaptive Sketch-and-Project Methods for Solving Linear Systems, SIAM Journal on Matrix Analysis and Applications 42(2):954–989, 2021 (journal version retitled “On Adaptive Sketch-and-Project for Solving Linear Systems”). DOI 10.1137/19M1285846, arXiv:1909.03604, https://arxiv.org/abs/1909.03604
5. Robert M. Gower, Filip Hanzely, Peter Richtárik, and Sebastian Stich, Accelerated Stochastic Matrix Inversion: General Theory and Speeding up BFGS Rules for Faster Second-Order Optimization, Advances in Neural Information Processing Systems 31 (NeurIPS 2018). arXiv:1802.04079, https://arxiv.org/abs/1802.04079
6. Robert M. Gower and Peter Richtárik, Linearly Convergent Randomized Iterative Methods for Computing the Pseudoinverse, arXiv:1612.06255, 2016, https://arxiv.org/abs/1612.06255
7. Nicolas Loizou and Peter Richtárik, Momentum and Stochastic Momentum for Stochastic Gradient, Newton, Proximal Point and Subspace Descent Methods, Computational Optimization and Applications 77:653–710, 2020. DOI 10.1007/s10589-020-00220-z, arXiv:1712.09677, https://arxiv.org/abs/1712.09677
8. Peter Richtárik and Martin Takáč, Stochastic Reformulations of Linear Systems: Algorithms and Convergence Theory, SIAM Journal on Matrix Analysis and Applications 41(2):487–524, 2020. DOI 10.1137/18M1179249, arXiv:1706.01108, https://arxiv.org/abs/1706.01108
9. Ling Tang, Yanjun Zhang, and Hanyu Li, On Sketch-and-Project Methods for Solving Tensor Equations, Applied Mathematics and Computation 511:129735, 2026. DOI 10.1016/j.amc.2025.129735, arXiv:2210.08241, https://arxiv.org/abs/2210.08241
10. Michał Dereziński and Elizaveta Rebrova, Sharp Analysis of Sketch-and-Project Methods via a Connection to Randomized Singular Value Decomposition, SIAM Journal on Mathematics of Data Science 6(1):127–153, 2024. DOI 10.1137/23M1545537, arXiv:2208.09585, https://arxiv.org/abs/2208.09585
11. Dmitry Kovalev, Robert M. Gower, Peter Richtárik, and Alexander Rogozin, Fast Linear Convergence of Randomized BFGS, arXiv:2002.11337, 2020, https://arxiv.org/abs/2002.11337
12. Ning Zheng and Tianwen Li, A Unified Convergence Analysis of Random Sketch Methods for Rank Deficient Linear Systems, Numerical Algorithms 101(1):703–739, 2026 (published online December 2025). DOI 10.1007/s11075-025-02297-1, https://doi.org/10.1007/s11075-025-02297-1
13. Julianne Chung, Matthias Chung, J. Tanner Slagel, and Luis Tenorio, Sampled Limited Memory Methods for Massive Linear Inverse Problems, Inverse Problems 36(5):054001, 2020. DOI 10.1088/1361-6420/ab77da, arXiv:1912.07962, https://arxiv.org/abs/1912.07962
14. Peter Richtárik and Martin Takáč, On Optimal Probabilities in Stochastic Coordinate Descent Methods, Optimization Letters 10(6):1233–1243, 2016. DOI 10.1007/s11590-015-0916-1, arXiv:1310.3438, https://arxiv.org/abs/1310.3438
15. Yansheng Su, Deren Han, Yun Zeng, and Jiaxin Xie, On the Convergence Analysis of the Greedy Randomized Kaczmarz Method, Numerical Algorithms, published online 2026. DOI 10.1007/s11075-026-02352-5, arXiv:2307.01988, https://arxiv.org/abs/2307.01988
16. Nicolas Loizou and Peter Richtárik, Revisiting Randomized Gossip Algorithms: General Framework, Convergence Rates and Novel Block and Accelerated Protocols, IEEE Transactions on Information Theory 67(12):8300–8324, 2021. DOI 10.1109/TIT.2021.3113285, arXiv:1905.08645, https://arxiv.org/abs/1905.08645