Sequential Preconditioned Conjugate Gradient Method for Linear Statistical Models

Dong-Yue Xie; Guan-Yu Chen; Xi Yang; Zun-Hao Zheng

REVIEW 2 major objections 3 minor 15 references

Sequential Preconditioned Conjugate Gradient Method for Linear Statistical Models

T0 review · 2 major / 3 minor · reviewed 2026-08-01 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read Sequential sketching reaches full-data least-squares accuracy in fewer iterations

desk verdict The sequential warm-started sketched PCG setup is a genuinely new combination with mostly sound theory, but the per-subproblem iteration bound rests on an unstated variance condition, so the efficiency claims are not fully proven; still, the method deserves a serious referee. read the letter →

arxiv 2607.25272 v1 pith:DL6UKUWO submitted 2026-07-28 math.NA cs.NAstat.ML

Guan-Yu Chen , Dong-Yue Xie , Xi Yang , Zun-Hao Zheng This is my paper · ORCID

classification math.NAcs.NAstat.ML MSC 65F2065F1062J05

keywords sequentialpreconditionedconjugategradientleast-squaresestimationrandomizedsketchingsubspaceembeddingwarm-startingordinaryleastsquarescomputationalcomplexity

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper proposes SPCG, a two-stage solver for large least-squares problems that solves a chain of increasingly large sketched subproblems with preconditioned conjugate gradients, using each solution to warm-start the next, and then applies a short full-data refinement. Its central claim is that this chain can reach the prediction accuracy of the ordinary least-squares estimator with far less work than solving on the full data from scratch. The main theorem bounds the final prediction error by an exponentially decaying initialization term plus the intrinsic error of the last sketch plus the OLS error, so as full-data iterations grow the method matches OLS. The paper also proves that only a few PCG iterations are needed per subproblem, and shows numerically that SPCG reaches the target with fewer iterations and less CPU time than full-data PCG and a competitor baseline.

What carries the argument

The load-bearing object is a nested, warm-started sequence of sketched least-squares systems with a common randomized preconditioner. The paper forms a fixed sketched Hessian Ĥ = (ŜX)^⊤(ŜX) once, then runs PCG on each subproblem (S_i X, S_i Y); the CG contraction on the preconditioned normal equations gives a per-iteration factor θ=1/4, while the subspace embedding gives the intrinsic accuracy level δ_i. Because subproblem i stops at δ_i instead of numerical zero, its output is a high-quality starting point for subproblem i+1; a final full-data PCG pass converts the accumulated error into OLS-level error. The iteration bounds rely on the ratio r_{i−1,i}=√((m_i−d)/(m_{i−1}−d)) and the doubli

What would settle it

Compute, for a fixed design and a doubling sketch schedule, Monte Carlo estimates of Var‖X(β̃_i−β)‖ and Var(c_i‖X(β̂−β)‖); if the former is smaller than the latter at any stage, the δ_i ≤ r_{i−1,i}δ_i derivation collapses and the per-subproblem iteration bound no longer follows. A direct experimental probe would run SPCG with a slowly increasing sketch schedule and check whether the required iterations per subproblem exceed the theoretical α.

Watch

Extended reading notes

Core claim

Under a subspace-embedding condition for a randomized Hadamard sketch, SPCG's global prediction error obeys E‖X(β_T−β)‖ ≤ θ^T E‖X(β_0−β)‖ + Σ_{i=1}^K (1+θ^{a_i}) θ^{T−τ_i} δ_i + (1+θ^s)δ_OLS, with θ=1/4. If each sketch level's intrinsic error δ_i stays within a constant C of the OLS error, the final bound becomes θ^T times the initial error plus (1+θ^s(1+5C/3))δ_OLS; as the full-data refinement length s grows, the expected prediction loss therefore tends to δ_OLS, the error of the OLS estimator itself. Each subproblem need only be solved to its own intrinsic accuracy, roughly O(log(1/ω)) iterations, because its output is a warm start for the next larger subproblem. The total cost is dominate

Load-bearing premise

The iteration-count theorem hinges on an unproved variance comparison in Appendix A.4: the sketched estimator's prediction error must have at least as much variance as a rescaled copy of the OLS error, a condition used to convert second-moment bounds into the first-moment bounds δ_i ≤ r_{i−1,i}δ_i but not stated as a hypothesis of Theorem 3.3.

Editorial extensions

If this is right

If the bounds hold, SPCG is statistically equivalent to OLS in prediction loss: the error decays to (1+θ^s(1+5C/3))δ_OLS, and letting the full-data refinement length grow gives δ_OLS.
Most PCG iterations run on systems of size m_i ≪ N, so the solver cost is dominated by the one-time transform O(Nd log N); with a sparse sketch the initialization cost can be reduced to O(Nd).
The per-subproblem iteration count is bounded by a constant α independent of N whenever sketch sizes double, so the number of solver iterations does not grow with sample size.
In the paper's experiments, SPCG reaches the OLS target with 1.25–1.76× fewer iterations and 2.10–2.83× less CPU time than full-data PCG, and it remains faster even on ill-conditioned problems with condition number 10^8.
Because the final full-data pass starts from an already accurate warm start, only a few full-data iterations are needed, making the method suitable for very large N where full-data Krylov methods are costly.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

If the variance condition in Appendix A.4 can be replaced by a direct finite-sample bound, the iteration-count guarantees would extend to designs far from the Gaussian/asymptotic regime where the second-moment prediction-efficiency formula is currently used.
The framework suggests a natural adaptive stopping rule: increase sketch size while an inexpensive error estimate says the warm-started iterate is still above target; the paper lists automatic sketch-schedule selection as an open problem.
Because the stagewise accuracy condition is generic, the same warm-start chain could be applied to other iterative solvers or regularized objectives, potentially carrying the same speedups to ridge or generalized least squares.
The numerical success of the sparse-sketch variant indicates that the conservative theoretical sketch-size requirements may be much larger than what is needed in practice, so the practical cost of the method may be lower than the stated O(Nd log N) bound.

Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

New here is the laddered warm-start scheme: solve a sequence of sketched least-squares problems with doubling sketch sizes, each PCG solve warm-started from the previous, then finish on the full data using the same sketched Hessian preconditioner. The individual pieces are known—SRHT/CountSketch embeddings, SAP preconditioning, warm starts—but the sequential multi-sketch structure with a unified preconditioner is, as far as I can tell, new. The nested SRHT construction that produces all sketches from one Hadamard transform is a nice practical detail. Theorems 3.1 and 3.2 are coherent: the CG contraction, the embedding event, and the triangle inequality chain line up, and the prose around (5) is honest that the simplified bound is conditional.

The soft spot is exactly where the stress-test puts it. Theorem 3.3, which is what guarantees "only a few PCG iterations per subproblem," is proven in A.4 using two assumptions stated in the appendix but not in the theorem hypothesis: ξ_i > m_i/N > d/N > γ and Var(Z_i) ≥ Var(c_i Z_OLS). The second-moment PE bound from Lemma A.7 gives E Z_i^2 ≤ c_i^2 E Z_OLS^2. To conclude E Z_i ≤ c_i E Z_OLS you need the variance condition, and it is simply asserted, no proof or reference. If it fails, the stagewise accuracy bound has no derivation, and Theorem 3.4's O(Nd log N) complexity estimate inherits the problem. That is a real gap in the efficiency argument, not a cosmetic one. The main statistical convergence claim (Theorem 3.2) does not depend on this, so the method's asymptotic OLS accuracy is on solid ground.

Two smaller points. The display √(5/9) < 1/4 in the proof of Theorem 3.1 is wrong as written (√(5/9) ≈ 0.745 > 0.25), but the preceding line has the correct inequality 2ε√((1+ε)/(1−ε)) < 1/4 for ε < 1/9, so it is a harmless typo. And the experiments set a_i = 3 based on a runtime plot; the theoretical lower bounds in Figure 1 are consistent with that, but there is no systematic sensitivity study, and the tuned schedule in Section 4.3 is explicitly empirical. The comparisons against IDS and PCG are honest, and the speedups are real but modest—1.4–1.9× versus IDS, 2.1–2.8× versus PCG.

I think this deserves a serious referee. The gap is concrete and fixable: either prove the variance condition under stated assumptions (it may follow from normality or from Marchenko–Pastur asymptotics with more careful finite-sample control) or make it an explicit hypothesis of Theorem 3.3. I would send it out with the request to address that condition and the sensitivity of a_i. For someone working on randomized least-squares solvers, this is worth reading and citing.

Referee Report

2 major / 3 minor

Summary. The paper proposes SPCG, a two-stage randomized iterative method for ordinary least-squares estimation in large-scale linear statistical models. The method solves a sequence of sketched least-squares subproblems with increasing sketch sizes, using PCG with a shared sketched preconditioner, warm-starting each subproblem from the previous solution, and then performs a final refinement on the full problem. The authors prove a per-subproblem contraction bound (Theorem 3.1), a global convergence bound showing that the final iterate attains OLS-level prediction accuracy as the number of full-data iterations grows (Theorem 3.2), per-subproblem iteration bounds (Theorem 3.3), and a complexity estimate with dominant cost O(Nd log N) (Theorem 3.4). Numerical experiments compare SPCG with Iterative Double Sketching and full-data PCG, reporting fewer iterations and less CPU time.

Significance. If the per-subproblem iteration bounds are valid, SPCG is an attractive combination of sketch-and-solve and sketch-and-precondition ideas, and the numerical results support its practical efficiency. Theorems 3.1 and 3.2 are coherent: the CG contraction, embedding event, and triangle inequality steps align, and the global OLS-accuracy claim follows from Theorem 3.2. The paper also gives clear experimental comparisons and a careful cost breakdown. However, the key 'few iterations per subproblem' claim rests on Theorem 3.3, whose proof relies on an unproved variance inequality and on a finite-sample use of an asymptotic prediction-efficiency lemma. This is a load-bearing gap, but it appears fixable by adding explicit hypotheses or supplying the missing argument.

major comments (2)

[Appendix A.4, proof of Theorem 3.3] The transition δ_{i-1} ≤ r_{i-1,i}δ_i is derived from the second-moment bound E Z_{i-1}^2 ≤ r_{i-1,i}^2 E Z_i^2 together with the variance inequality Var(Z_{i-1}) ≥ r_{i-1,i}^2 Var(Z_i). The variance inequality is asserted as a 'condition' but is never proved, cited, or included as a hypothesis of Theorem 3.3. It is not a consequence of the SRHT embedding properties used elsewhere in the paper. Since Theorem 3.3 is the only result that yields the logarithmic per-subproblem iteration counts, and Theorem 3.4 inherits that conclusion, the central efficiency claim is currently unsupported. The variance condition must either be proved for the assumed noise model or stated explicitly as an assumption in Theorem 3.3 and propagated to Theorem 3.4.
[Appendix A.4, Lemma A.7 and finite-sample use] Lemma A.7 is an asymptotic result: it assumes d/N→γ and m_i/N→ξ_i as N→∞. In the proof of Theorem 3.3, this asymptotic statement is converted into the finite-sample inequality PE(i) ≤ (N−d)/(m_i−d) with no error term or non-asymptotic derivation. The displayed condition 'ξ_i > m_i/N > d/N > γ' also mixes the limiting quantity ξ_i with finite ratios. Because Theorems 3.3 and 3.4 are finite-N statements, a finite-sample version of the prediction-efficiency bound is needed. Alternatively, Theorem 3.3 should be formulated as conditional on the finite-sample PE bound, and the complexity theorem should be adjusted accordingly.

minor comments (3)

[Appendix A.2, proof of Theorem 3.1] The displayed inequality √(5/9) < 1/4 is numerically false. The intended argument is 2ε√((1+ε)/(1−ε)) < 1/4 for ε<1/9, which does hold. Since the text immediately uses the product form with ε^{a_i−1}, the conclusion of Theorem 3.1 is unaffected.
[Theorem 3.2, Eq. (5)] The simplified bound (5) is conditional on δ_i ≤ C δ_OLS, and the paragraph after the theorem notes that one may take C = max_i δ_i/δ_OLS. With that choice, (5) is essentially a rearrangement rather than a closed-form predictive bound. If the intent is to state an explicit OLS-accuracy guarantee, a concrete upper bound on δ_i/δ_OLS in terms of the sketch sizes is needed; otherwise (4) is the substantive convergence result.
[Appendix A.4, first displayed condition] The condition 'ξ_i > m_i/N > d/N > γ' mixes the asymptotic limit ξ_i with finite-N quantities. Please clarify, for example by replacing ξ_i with m_i/N in the finite-sample condition or by writing the asymptotic limits separately.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity: proofs rely on external lemmas; the unproved variance condition in Appendix A.4 is a proof gap, not a circular step.

full rationale

SPCG's derivation chain is not circular. Theorem 3.1's per-subproblem contraction follows from the SRHT subspace embedding guarantee (Lemma 2.2, [14]) and the standard PCG bound (Lemma A.2, [10]); the intrinsic errors δ_i = E||X(β̃_i − β)|| are defined from the independent sketched estimators, not from the iterates being bounded. Theorem 3.2 is a recursive triangle-inequality chaining of Lemmas A.5–A.6; the C appearing in (5) is introduced as C = max_i δ_i/δ_OLS for a simplified display and is not a fitted parameter. Theorem 3.3's iteration bound uses the external asymptotic prediction-efficiency result Lemma A.7 ([2]) together with two conditions stated in Appendix A.4. The second condition, Var(Z_i) ≥ Var(c_i Z_OLS), is asserted without proof or reference and is not listed as a hypothesis of Theorem 3.3; this is a genuine rigor gap in the 'few iterations' claim, and it should be flagged as missing support. But it is not circular: the condition is not equivalent to the stagewise-accuracy conclusion, no self-citation is used to establish it, and the main OLS-accuracy bound (Theorem 3.2) does not depend on it. There are no self-citations of the present authors in the reference list, no ansatz smuggled in by citation, and no renaming of an externally known result. Numerical comparisons are against external baselines (IDS, full-data PCG). Hence no significant circularity.

Assumptions & free parameters 3 free parameters · 6 assumptions · 0 invented entities

The only invented items are algorithm parameters and an unproved variance/spectral condition used to derive stagewise iteration bounds. No new mathematical or physical entities are introduced.

free parameters (3)

inner iteration count a_i = 3 (experiments); lower bound from Theorem 3.3
In Section 4.1, a_i=NoI is chosen uniformly for all sketched subproblems by minimizing runtime on the N=2^20, d=2^6, κ=10^4 test problem (Fig.1b). The theoretical bounds give only a lower bound, not an optimal value.
sketch size schedule m_i = m_1=8d, m_{i+1}=2m_i, m_K=N/2 (experiments); m_K=N/2 and K=O(1) in Theorem 3.4
Doubling schedule is prescribed by the analysis; the starting size m_1=8d and terminal size N/2 are hand-chosen; Section 4.3 tunes the schedule empirically, showing the choice affects runtime.
preconditioner sketch size r = 6d (experiments); r=O(d log d) in theory
The sketched Hessian rank r=6d is a hand-set implementation parameter; the theory only requires r ≥ c ε^{-2}(d+log(N/δ))log(ed/δ).

assumptions (6)

standard math SRHT subspace embedding guarantee (Lemma 2.2, from [14])
Basis of all subsequent conditioning bounds.
standard math CG convergence estimate (Lemma A.2, from [10])
Used in every stage to bound PCG error contraction.
domain assumption Linear model Y=Xβ+ζ with zero-mean noise, covariance σ²I_N, and X full column rank (Eq.1)
Defines the OLS target and δ_OLS.
domain assumption N and all m_i are powers of two (SRHT construction)
SRHT requires N power of two; Theorems 3.1/3.4 assume it.
ad hoc to paper Finite-sample PE bound PE(i) ≤ (N−d)/(m_i−d) inferred from asymptotic Lemma A.7, plus Var(Z_i) ≥ Var(c_i Z_OLS) and ξ_i > m_i/N > d/N > γ
Used in A.4 to get δ_i ≤ c_i δ_OLS and interstage bounds; not proven and not listed in Theorem 3.3.
ad hoc to paper δ_i ≤ C δ_OLS for a constant independent of i (Theorem 3.2 bound (5))
For a finite sequence one can define C=max_i δ_i/δ_OLS, so (5) is conditional on a problem-dependent C.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Sequential Preconditioned Conjugate Gradient Method for Linear Statistical Models." pith.science (2026). https://pith.science/paper/DL6UKUWO

@misc{pith2026260725272,
  author       = {Pith},
  title        = {Pith review of: Sequential Preconditioned Conjugate Gradient Method for Linear Statistical Models},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/DL6UKUWO}},
  note         = {Machine review of arXiv:2607.25272}
}

read the original abstract

We propose a randomized iterative method for the ordinary least-squares estimation problem in large-scale linear statistical models, namely the Sequential Preconditioned Conjugate Gradient Method (SPCG). SPCG constructs a sequence of sketched least-squares subproblems with increasing sketch sizes, applies PCG as the inner solver, and warm-starts each subproblem from the previous solution. A final refinement stage is then performed on the full-scale problem. Since most iterations are carried out on smaller subproblems, the overall computational cost is significantly reduced. We establish the convergence theory, prove that SPCG attains OLS prediction accuracy, and derive per-subproblem iteration bounds and complexity estimates. Numerical experiments show that SPCG reaches the target prediction accuracy with fewer iterations and less CPU time than full-data PCG and Iterative Double Sketching (IDS).

Figures

Figures reproduced from arXiv: 2607.25272 by the authors.

**Figure 2.** FIG. 2 [PITH_FULL_IMAGE:figures/full_fig_p007_2.png] view at source ↗

**Figure 3.** FIG. 3 [PITH_FULL_IMAGE:figures/full_fig_p007_3.png] view at source ↗

Figures from the paper (2 more)

**Figure 4.** Figure 4: FIG. 4 [PITH_FULL_IMAGE:figures/full_fig_p009_4.png]

**Figure 5.** Figure 5: FIG. 5 [PITH_FULL_IMAGE:figures/full_fig_p009_5.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

15 extracted references · 2 linked inside Pith

[1]

Blendenpik: supercharging LAPACK’s least-squares solver

Avron, H., Maymounkov, P., Toledo, S., 2010. Blendenpik: supercharging LAPACK’s least-squares solver. SIAM Journal on Scientiﬁc Computing 32, 1217–1236

2010
[2]

Asymptotics for sketching in least squares regression, in: Advances in Neural Information Processing Systems, pp

Dobriban, E., Liu, S., 2019. Asymptotics for sketching in least squares regression, in: Advances in Neural Information Processing Systems, pp. 3670–3680

2019
[3]

RandNLA: randomized numerical linear algebra

Drineas, P., Mahoney, M.W., 2016. RandNLA: randomized numerical linear algebra. Communications of the ACM 59, 80–90

2016
[4]

Matrix Computations

Golub, G.H., Van Loan, C.F., 2013. Matrix Computations. Johns Hopkins University Press, Baltimore

2013
[5]

Optimal randomized ﬁrst-order methods for least-squares problems, in: International Conference on Machine Learning, PMLR

Lacotte, J., Pilanci, M., 2020. Optimal randomized ﬁrst-order methods for least-squares problems, in: International Conference on Machine Learning, PMLR. pp. 5587–5597

2020
[6]

Fast convex quadratic optimization solvers with adaptive sketching-based preconditioners

Lacotte, J., Pilanci, M., 2021. Fast convex quadratic optimization solvers with adaptive sketching-based preconditioners. arXiv preprint arXiv:2104.14101

arXiv 2021
[7]

Randomized numerical linear algebra: foundations and algorithms

Martinsson, P.G., Tropp, J.A., 2020. Randomized numerical linear algebra: foundations and algorithms. Acta Numerica 29, 403–572

2020
[8]

LSRN: a parallel iterative solver for strongly over- or underdetermined systems

Meng, X., Saunders, M.A., Mahoney, M.W., 2014. LSRN: a parallel iterative solver for strongly over- or underdetermined systems. SIAM Journal on Scientiﬁc Computing 36, C95–C118

2014

Show all 15 references

[9]

Randomized numerical linear algebra: a perspective on the ﬁeld with an eye to software

Murray, R., Demmel, J., Mahoney, M.W., Erichson, N.B., Melnichenko, M., Malik, O.A., Grigori, L., Luszczek, P., Derezinski, M., Lopes, M.E., Liang, T., Luo, H., Dongarra, J., 2023. Randomized numerical linear algebra: a perspective on the ﬁeld with an eye to software. arXiv pr...

2023 arXiv
[10]

Numerical Optimization

Nocedal, J., Wright, S.J., 2006. Numerical Optimization. Springer, New York

2006
[11]

Iterative Hessian sketch: fast and accurate solution approximation for constrained least-squares

Pilanci, M., Wainwright, M.J., 2016. Iterative Hessian sketch: fast and accurate solution approximation for constrained least-squares. Journal of Machine Learning Research 17, 1–38. Page 13 of 14 Sequential preconditioned conjugate gradient method

2016
[12]

A statistical perspective on randomized sketching for ordinary least-squares

Raskutti, G., Mahoney, M.W., 2016. A statistical perspective on randomized sketching for ordinary least-squares. Journal of Machine Learning Research 17, 1–31

2016
[13]

A First Course in Probability

Ross, S.M., 2020. A First Course in Probability. Pearson, Harlow

2020
[14]

Iterative double sketching for faster least-squares optimization, in: International Conference on Machine Learning, PMLR

Wang, R., Ouyang, Y., Xu, W., 2022. Iterative double sketching for faster least-squares optimization, in: International Conference on Machine Learning, PMLR. pp. 22935–22963

2022
[15]

Sketching as a tool for numerical linear algebra

Woodruﬀ, D.P., 2014. Sketching as a tool for numerical linear algebra. Foundations and Trends in Theoretical Computer Science 10, 1–157. Page 14 of 14

2014

Pith tools

Reviewed August 1, 2026 · model on record in the stance chip above.

[1] [1]

Blendenpik: supercharging LAPACK’s least-squares solver

Avron, H., Maymounkov, P., Toledo, S., 2010. Blendenpik: supercharging LAPACK’s least-squares solver. SIAM Journal on Scientiﬁc Computing 32, 1217–1236

2010

[2] [2]

Asymptotics for sketching in least squares regression, in: Advances in Neural Information Processing Systems, pp

Dobriban, E., Liu, S., 2019. Asymptotics for sketching in least squares regression, in: Advances in Neural Information Processing Systems, pp. 3670–3680

2019

[3] [3]

RandNLA: randomized numerical linear algebra

Drineas, P., Mahoney, M.W., 2016. RandNLA: randomized numerical linear algebra. Communications of the ACM 59, 80–90

2016

[4] [4]

Matrix Computations

Golub, G.H., Van Loan, C.F., 2013. Matrix Computations. Johns Hopkins University Press, Baltimore

2013

[5] [5]

Optimal randomized ﬁrst-order methods for least-squares problems, in: International Conference on Machine Learning, PMLR

Lacotte, J., Pilanci, M., 2020. Optimal randomized ﬁrst-order methods for least-squares problems, in: International Conference on Machine Learning, PMLR. pp. 5587–5597

2020

[6] [6]

Fast convex quadratic optimization solvers with adaptive sketching-based preconditioners

Lacotte, J., Pilanci, M., 2021. Fast convex quadratic optimization solvers with adaptive sketching-based preconditioners. arXiv preprint arXiv:2104.14101

arXiv 2021

[7] [7]

Randomized numerical linear algebra: foundations and algorithms

Martinsson, P.G., Tropp, J.A., 2020. Randomized numerical linear algebra: foundations and algorithms. Acta Numerica 29, 403–572

2020

[8] [8]

LSRN: a parallel iterative solver for strongly over- or underdetermined systems

Meng, X., Saunders, M.A., Mahoney, M.W., 2014. LSRN: a parallel iterative solver for strongly over- or underdetermined systems. SIAM Journal on Scientiﬁc Computing 36, C95–C118

2014