Subspace-constrained randomized coordinate descent for linear systems with good low-rank matrix approximations

1 Introduction

• The motivation is that dense KRR systems require O(n^2) storage and O(n^3) direct-solve arithmetic, while RCD can become very slow when spectral outliers worsen conditioning.Low-rank approximations are attractive because many data matrices exhibit rapid spectral decay, and RPCholesky uses O(d^2n) arithmetic, O(dn) entry evaluations, and O(dn) storage.Source: §1 · pp. 1–3

1.2 Main results

• SC-RCD restricts randomized block coordinate descent to an affine subspace defined by a rank-d Nyström approximation, making convergence depend on the residual spectrum rather than A’s largest spectral outliers.The Nyström approximation can be computed with RPCholesky, while residual-diagonal sampling selects the iteration blocks.Source: §1.2 · pp. 3–4
• SC-RCD requires only the precomputed auxiliary matrix C and sampled columns A_:,J, avoiding storage of A or its residual and using O((ℓ + d)n) arithmetic per iteration.This supports lightweight, memory-efficient iterations when d ≪ n and ℓ is at most O(√n).Source: §1.2 · pp. 4–5
• Theorem 1.3 establishes linear expected convergence, conditional on RPCholesky approximation quality, with a rate governed by eigenvalues after removing the largest spectral outliers.The result holds with high probability over the randomized approximation and bounds expected relative A-norm error during the iterative phase.Source: §1.2 · pp. 5–6
• Theorem 1.7 gives O(nd log(1/ε)) + O(n^2 · κ̄_r(A) log(1/ε)) entry evaluations and O(n(d + ℓ)) storage for approximate solution computation.The bound depends on approximation rank, block size, and the normalized tail condition number, benefiting matrices with flat spectral tails.Source: §1.2 · p. 7
• For KRR with κ̄_r(A) = O(1), r = O(√n/log n), d = O(√n), and ℓ = O(√n), SC-RCD achieves ε-relative error in O(n^2 log(1/ε)) arithmetic operations.This setting uses A = K + λI with sufficiently large regularization and rapidly decaying kernel spectrum.Source: §1 · p. 1, §1.2 · p. 7
• The general subspace-constrained sketch-and-project framework covers randomized Kaczmarz, randomized coordinate descent, and randomized Newton, with projected dynamics that effectively precondition the system.Its convergence depends on AP_B rather than A, where the projector P_B constrains iterates to the selected affine subspace.Source: §1.2 · pp. 7–8
• The framework proves that i.i.d. block sketches improve iteration complexity at least linearly with block size, including a new explicit block-size dependence for randomized block Kaczmarz.The result applies to blocks sampled with replacement according to squared row norms and is a special case of the subspace-constrained framework.Source: §1.2 · pp. 8–9

1.3 Related works

• Related work places SC-RCD alongside randomized coordinate descent, randomized Kaczmarz, accelerated sketch-and-project methods, CD++, and preconditioned conjugate gradient solvers for approximately low-rank systems.Low-rank Nyström and interpolative-decomposition approximations provide coordinate-aligned subspaces that can be enforced by updating fixed coordinate subsets.Source: §1.3 · pp. 9–10
• SC-RCD combines an explicitly learned low-rank Nyström subspace constraint with a simple iterative solver, avoiding the entire-matrix storage assumed by CD++ while retaining flexibility for related algorithmic innovations.For certain flat-tailed systems, SC-RCD has a similar complexity bound to CD++, whereas CD++ can have better normalized-tail dependence through Nesterov momentum.Source: §1.3 · pp. 9–10

1.4 Organization

• The paper develops the abstract subspace-constrained framework in Section 2, specializes its analysis to SC-RCD for positive semidefinite systems in Section 3, and evaluates it numerically in Section 4.Appendices provide omitted proofs, a least-squares extension, and additional experiments.Source: §1.4 · pp. 10–11

2.1 The subspace-constrained sketch-and-project framework

• Subspace-constrained sketch-and-project confines iterates to {x ∈ Rn : QAx = Qb} while projecting each iterate onto a sketched system.The constraint is computationally useful when the smaller system can be solved efficiently and projections onto the relevant nullspace are efficient.Source: §2.1 · pp. 11–12

2.2 Update rule and comparison with unconstrained sketch-and-project

• Its closed-form updates use projectors onto null(QAB−1/2) and the sketched subspace, while preserving an invariant error subspace.With Q = 0, the formulas recover standard sketch-and-project updates.Source: §2.2 · p. 12
• Each constrained iteration decreases error by at least as much as the corresponding unconstrained update using the same sketch.The result follows from the projector structure and Pythagorean error decomposition.Source: §2.2 · pp. 12–13
• The convergence rate depends on the eigenvalue spectrum of E[Z], under an exactness condition ensuring that all relevant error directions are covered.If the sketch distribution does not cover the entire space, unresolved components can remain in the error.Source: §2.2 · pp. 13–14
• Using the B-inner-product geometry, constrained updates parallel unconstrained sketch-and-project with the projected matrix APB replacing A.This makes the algorithmic constraint an implicit preconditioner; for SC-RCD, Q is chosen so APB approximates A well in trace norm.Source: §2.2 · p. 14

2.3 Convergence rate and block size

• The block-rate analysis decomposes a multi-row projector into rank-one updates and applies the general convergence theorem.This decomposition preserves the exactness condition across intermediate block sizes.Source: §2.3 · pp. 14–15
• For i.i.d. block sketches, increasing the block size improves iteration complexity at least linearly, so τ single-row iterations require at most τ/ℓ block iterations.The framework covers block randomized Kaczmarz and coordinate descent, including sampling schemes based on projected row norms.Source: §2.3 · pp. 15–16
• The framework specializes to subspace-constrained randomized Kaczmarz, whose rate is governed by the projected matrix AP and reduces to standard block Kaczmarz when P = I.The projected singular-value expression yields the corresponding convergence bound.Source: §2.3 · pp. 16–17

3 Analysis of subspace-constrained randomized coordinate descent

• SC-RCD restricts randomized coordinate-descent iterates to an affine subspace defined by a column Nyström approximation, while remaining well-defined for positive semidefinite matrices.The approximation is formed from a pivot set S, and the update formula remains valid even when A is rank deficient.Source: §3.1 · p. 18

3.1 Convergence rate of SC-RCD

• Diagonal sampling proportional to diag(A◦) yields a block-size-dependent convergence rate controlled by the residual matrix A◦.The same framework extends to singular systems by tracking the A-seminorm error, which is independent of the chosen solution.Source: §3.1 · pp. 18–19
• Uniform block sampling has a related rate depending on the diagonal-normalized residual spectrum, and sampling without replacement can only improve the per-iteration error decrease.Uniform sampling is therefore a simpler alternative whose effectiveness depends on residual incoherence.Source: §3.1 · p. 19, §3.2 · pp. 20, 22
• The residual-spectrum analysis improves over standard RCD when tr(A◦) is much smaller than tr(A), as can occur when RPCholesky produces a good low-rank approximation.Residual eigenvalue interlacing preserves the relevant numerator while a smaller residual trace improves the denominator; RPCholesky guarantees support this reduction.Source: §3.1 · pp. 19–20

3.2 Implementation and complexity of SC-RCD

• RPCholesky initialization costs O(d^2n) arithmetic operations, O(dn) entry evaluations, and O(dn) storage, while each SC-RCD iteration costs O(dn + ℓn + ℓ^3 + ℓ^2d) operations.The initial iterate can be computed in O(d^2) flops, and the residual vector can then be formed in O(dn) flops rather than O(n^2).Source: §3.2 · pp. 20–21
• With RPCholesky boosting, the overall solver requires O(n^2 · κ̄_r(A) log(1/ε)) entry evaluations and the stated arithmetic complexity to reach the target expected error.The guarantee combines the low-rank approximation event with the SC-RCD convergence bound.Source: §3.2 · p. 21

Synthetic psd system

• On the synthetic system, SC-RCD converges effectively because the rank-500 residual is much better conditioned than A, while standard RCD and CD++ need larger blocks to capture the leading spectral outliers.The comparison uses 8,192-dimensional systems with approximate rank r = 400 and block size ℓ = 500.Source: §4 · p. 22, Figure 4.1 · p. 23
• SC-RCD error decreases as both block size ℓ and approximation rank d increase, with marked improvement when their combined effects capture the large spectral outliers.The observed gains exceed the basic theoretical bound, consistent with an implicit low-rank approximation effect from larger blocks.Source: §4 · p. 22, Figure 4.2 · p. 23

Implementation considerations

• SC-RCD adds O(nd) work per iteration relative to RCD, but this is not dominant when d, ℓ = O(√n), and its runtime is close to RCD for equal block sizes.CD++ iterations are somewhat faster, but its theoretical guarantee requires storing and preprocessing A or assuming suitable incoherence.Source: §4 · pp. 22–23

KRR problem with fast spectral decay

• On hls4ml_lhc_jets, SC-RCD improves substantially over RCD and converges faster than PCG with the same d = 1,000 approximation, although PCG benefits more from increasing d.The experiment uses n = 100,000 and λ = 10^-9n, with kernel entries evaluated on demand rather than storing K.Source: §4 · pp. 24–25, Figure 4.3 · p. 24

KRR problem with slower spectral decay

• On sensorless, uniform SC-RCD sampling converges fastest, whereas diagonal sampling performs poorly despite helping on rapidly decaying spectra.The dataset has slower spectral decay, making a good low-rank approximation more difficult; uniform SC-RCD still improves over RCD.Source: §4 · p. 25, Figure 4.4 · p. 25

• SC-RCD combines randomized block coordinate descent with a rank-d matrix approximation, which can be computed efficiently using RPCholesky.The method is designed for solving psd linear systems.Source: §5 · p. 26
• SC-RCD uses O(nd) memory to obtain an ε-relative error solution.Source: §5 · p. 26
• The method is effective for large-scale, dense systems with rapid spectral decay, including systems arising in kernel ridge regression.Source: §5 · p. 26
• Numerical experiments support the theoretical convergence and effectiveness results.Source: §5 · p. 26
• Future work includes improving condition-number dependence, learning subspaces for different spectral regions, and extending the framework to other sketch-and-project and nonlinear methods.Source: §5 · p. 26

• The constrained update is derived by changing variables to w = B^1/2(x^k − x), imposing the sketch and affine constraints, and applying Lagrange multipliers.Source: A Subspace-constrained… · p. 30
• Projector identities imply the optimal transformed step lies in range(PB^-1/2A^TS^T), so the associated orthogonal projector fixes that step.Source: A Subspace-constrained… · p. 30
• Substituting the projected step back into the transformed variables yields the update rule and fixed-point iteration used in the subspace-constrained analysis.Source: A Subspace-constrained… · p. 30
• The proof starts from an initialization satisfying the affine constraint and shows the transformed initial error lies in the projector-defined subspace.Source: A Subspace-constrained… · p. 30
• Because every iterate continues to satisfy the affine constraint, induction keeps each transformed error in range(PB^-1/2A^T).Source: A Subspace-constrained… · p. 30

Algorithm 2 SC-RCD: least squares

• Algorithm 2 computes an approximate least-squares solution and residual using an approximation rank d and block size ℓ.Source: B Extension of SC-RCD for… · p. 31
• The implementation computes a partial QR approximation bA = QR and samples coordinate blocks J independently from a specified distribution.Source: B Extension of SC-RCD for… · p. 31

Extension to least squares

• Least squares reduces to normal equations, but SC-RCD can operate as a column-action method without explicitly forming A^TA.Source: B Extension of SC-RCD for… · p. 31

Column projection approximation

• The natural low-rank analogue is the column projection approximation ΠA,SA, whose selected-column span is represented by a partial QR factorization.Source: B Extension of SC-RCD for… · p. 31

Derivation of Algorithm 2

• The subspace constraint preserves (A_:,S)^T(Ax^k−b)=0, enabling residual updates through the residual matrix and an overdetermined least-squares solve.Source: Derivation of Algorithm 2 · pp. 31–32
• The residual update and auxiliary matrix can be implemented using sampled columns, with partial QR structure supporting efficient initialization and precomputation.Source: Derivation of Algorithm 2 · p. 32

Convergence

• The least-squares method has a convergence-rate theorem for fixed pivots, which can be combined with low-rank approximation guarantees.Source: Derivation of Algorithm 2 · p. 32

Experimental setup

• The experiments solve Gaussian-kernel ridge-regression systems with n = 20,000, λ = 10^-8n, and compare SC-RCD, RCD, CG, and PCG.SC-RCD and PCG use d = 1,000; SC-RCD and RCD use block size ℓ = 1,000.Source: C Additional numerical experiments · p. 32
• The figures report relative residual norms over 300 epochs using medians and 0.2/0.8-quantiles from 10 independent runs, together with eigenvalue spectra.Source: C Additional numerical experiments · p. 32

Rapid spectral decay

• For rapidly decaying kernel matrices, Figure C.1 presents convergence trajectories and eigenvalues for ACSIncome, Airlines_DepDelay_1M, cod-rna, and diamonds.Source: Figure C.1 · p. 33

Slower spectral decay

• For slower spectral decay, Figure C.2 presents convergence trajectories and eigenvalues for covtype.binary, creditcard, HIGGS, and SensIT Vehicle.Source: Figure C.2 · p. 34
• The slower-decay experiments specifically use uniformly sampled blocks in SC-RCD, which the caption describes as typically working quite well.Source: Figure C.2 · p. 34