Pith. sign in

REVIEW 3 major objections 6 minor 16 references

A Portable and Versatile Limited-Memory BFGS Implementation in PETSc/TAO

T0 review · 3 major / 6 minor · reviewed 2026-08-04 · deepseek-v4-flash

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

Pith's one-line read A new 'intermediate dense' formulation of limited-memory BFGS computes each quasi-Newton step with a single application of the base Hessian inverse and dense BLAS operations, avoiding the extra passes and recomputations that slow compact de

desk verdict The intermediate dense L-BFGS factorization is sound and genuinely useful; the 'most robust' performance claim outruns the evidence. read the letter →

arxiv 2608.00196 v1 pith:YOT6F4IC submitted 2026-07-31 cs.DC

Hansol Suh , Tobin Isaac , Alp Dener , Todd Munson , Hong Zhang , Richard Tran Mills This is my paper · ORCID
classification cs.DC MSC 65K0590C3090C5365Y0565Y10
keywords limited-memoryBFGSquasi-Newtonmethodsvariable-metricdenseBLASGPUaccelerationperformanceportabilityPETSc/TAO
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 claims that the standard ways of implementing limited-memory BFGS—the recursive two-loop recursion and the compact dense representation—have complementary weaknesses, and that a third representation, called the intermediate dense form, preserves the dense BLAS efficiency and low synchronization of the compact form without its penalty when the base inverse Hessian H0 changes between iterations. The representation factors the L-BFGS inverse Hessian into a product that applies H0 exactly once and does all other work with small dense matrices and matrix-vector products, using only the triangular matrix R = upper-triangle(S^T Y) and a cached S^T g product. If this claim holds, quasi-Newton optimizers with variable metrics (e.g., diagonal scaling) can enjoy dense-level performance on CPUs and GPUs without paying m applications of H0 per step, which the paper's microbenchmarks and one MNIST training run support. The paper also proves that no analogous intermediate dense form exists for the forward Hessian application B.mult, so that direction remains tied to the compact dense representation.

What carries the argument

The key machinery is the intermediate dense identity (6.17), obtained by multiplying out the blocks of the compact dense representation and noticing that the H0 term appears only once, as an oblique projection (I - S R^{-T} Y^T) H0 (I - Y R^{-1} S^T) plus a small rank correction S R^{-T} D R^{-1} S^T. All history-dependent data are packed into S, Y, and the m x m matrix R = htriu(S^T Y); R is stored in cyclic order matching the history vectors, and products with R are triangular solves via a three-step 'trsv cyclic' routine. This eliminates the need for the Q = H0 Y basis and Z = Y^T H0 Y matrix that the compact dense form must recompute whenever H0 changes.

What would settle it

Run the three implementations inside a complete line-search L-BFGS loop (including objective and gradient evaluation) with a diagonal variable-metric H0 at n=10^7, m=50 on a 4-GPU node; if the recursive two-loop recursion is consistently faster in wall-clock than the intermediate dense form, then the paper's central performance claim that the intermediate dense approach is 'most robust' would fail. The missing piece is any comparison with line search enabled.

Watch

Extended reading notes

Core claim

The central discovery is a factored representation of the L-BFGS inverse Hessian, equation (6.17), which expresses H_k^m as the product of two rank-m 'wings' around a block diagonal containing H0 and D = diag(d_i). Applied to a vector, it requires one solve with H0, four dense matrix-vector products with S and Y, and two triangular solves with the small m x m matrix R = htriu(S^T Y) (the history-order upper triangle of the s-y inner products), plus caching of S^T g to avoid recomputation in the line-search step. This yields O(mn) BLAS work and O(1) synchronization points per quasi-Newton step, with only one H0 application, whereas the compact dense form needs m applications and O(m^2n) BLAS

Load-bearing premise

The paper's performance conclusions rest on a stripped-down loop that skips line search and restricts the base matrix to two diagonal choices; if full optimizations change the timing balance, its recommendation could shift.

Editorial extensions

If this is right

  • Under a changing base inverse Hessian H0, a variable-metric L-BFGS step costs one H0 solve instead of m, so expensive base preconditioners become affordable in quasi-Newton loops.
  • The intermediate dense form matches the recursive method's BLAS work (roughly (8m+2)n flops) but with Level 2/3 BLAS and O(1) collectives, giving the best effective bandwidth at large n on both CPUs and GPUs.
  • When H0 is constant and cheap, the compact dense form remains the right choice; the two dense forms are complementary, not interchangeable, for B.solve.
  • No analogous intermediate dense form exists for B.mult: applying the forward approximate Hessian still requires the compact dense representation or recursive O(m^2n) updates, as proven by an orthonormal-set argument.
  • The same dense patterns have been extended to Broyden, bad Broyden, symmetric Broyden, and SR1 limited-memory matrices in the library, so the variable-metric benefit is not BFGS-specific.

Reading between the lines

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

  • Because the MNIST end-to-end test used only a diagonal H0, the speedup there (up to 4.7%) is a lower bound for applications whose H0 is an expensive preconditioner; the per-step savings should grow with the cost of each H0 solve.
  • The paper's own measurements show the fixed cases where recursive still wins (n=10^5, m=50 and n=10^6 with 4 CPU ranks) coincide with their single-rank storage of R; a variant that duplicates R on all ranks (allreduce instead of reduce+bcast) could plausibly close that gap.
  • The negative result for B.mult suggests that an inexact forward-Hessian application—for example, using fewer history pairs in the mult operator than in the solve operator—could be a useful compromise for trust-region subproblems, a direction the paper does not explore.
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.

Referee Report

3 major / 6 minor

Summary. The paper introduces an 'intermediate dense' representation of the limited-memory BFGS inverse Hessian, Eq. (6.17), and implements it in PETSc/TAO alongside the standard recursive two-loop recursion and the compact dense representation of Byrd, Nocedal, and Schnabel. The key algebraic claim is that the intermediate dense form applies H_k^m v with one B0^{-1} application and O(mn) BLAS work, avoiding the recomputation of Q and Z matrices that makes compact dense expensive in variable-metric settings. The authors provide a complexity analysis, a non-existence argument for an analogous intermediate dense B.mult, single-node microbenchmarks on Polaris and Frontier (CPU and GPU), and one MNIST training benchmark.

Significance. The algebraic derivation from (6.8) to (6.16)-(6.17) is correct, and Algorithm 9 implements it with the stated BLAS counts; this is a genuine and useful algorithmic contribution for PETSc/TAO. The paper also gives a clean complexity comparison and an interesting negative result for B.mult in Section 6.3.2. If the performance claims hold, the intermediate dense method would be a practical default for variable-metric L-BFGS. However, the strongest performance conclusion in Section 11 — that the intermediate dense approach is 'the most robust across architectures and base-matrix choices' — is not supported by the measurements as presented, because the microbenchmarks omit line search and objective evaluation, use only diagonal variable-metric H0, lack error bars, and the end-to-end MNIST result compares only one 'dense' variant against recursive BFGS.

major comments (3)
  1. [Section 7.2 and Section 11] The claim that the intermediate dense approach 'provides the most robust performance across architectures and base-matrix choices' is stronger than the evidence. The Tstep loop in Section 7.2 is explicitly only Update+Solve, with no line search and no objective/gradient evaluation, as the paper itself acknowledges. Figures 2-12 show single point estimates with no error bars or repeated runs. Line search adds vector operations and synchronization that can change the relative ranking of the algorithms, especially for the recursive method's sequential rank-1 reductions. Please either add end-to-end line-search benchmarks for all three implementations (or a defensible proxy that includes line-search costs) and report variability, or restrict the conclusion to the microbenchmark setting.
  2. [Section 7.2 and Section 3] The only variable-metric H0 tested is a diagonal matrix built with full-memory restricted Broyden scaling. This does not exercise the setting that motivates the intermediate dense algorithm: a general, expensive variable base matrix such as the multigrid preconditioner mentioned in Section 3. With a diagonal H0, B0.solve is a cheap vector scaling, so the compact dense method's additional B0.solve applications are not costly; the observed difference is dominated by BLAS flops rather than by the number of base-solve applications. The complexity argument for one B0.solve per step instead of m is clear, but the practical advantage in expensive variable-metric cases remains unquantified. I request at least one experiment with a non-diagonal variable base operator, or a cost model that combines the measured T(B0.solve) of a realistic preconditioner with the Table 6 complexity counts.
  3. [Section 8, Table 8] The MNIST result compares 'dense' with recursive BFGS, but it is not stated whether 'dense' is the intermediate dense or compact dense implementation; the text and table caption should identify it. More importantly, this is a single training problem, and the reported improvement (0.9-4.7%) is in solver time excluding forward/backward propagation. It does not substantiate the conclusion that the intermediate dense approach is most robust, nor does it provide an intermediate-dense-versus-compact-dense comparison. If the end-to-end claim is load-bearing, please add comparisons of all three variants on more than one workload, or soften the conclusion accordingly.
minor comments (6)
  1. [Section 7.2] After the sentence 'included in the PETSc test suite. 4' there is an orphan fragment '1ksp ksp utils lmvm tests-solve performance 0.' It appears to be a leftover command; remove or format as code.
  2. [Figures 2-12] The meaning of the legend labels '1/2 BS' and '1/4 BS' is not defined. The text refers to 1 and 4 MPI ranks/GPUs, but the figure labels are cryptic and should be explained in the captions or the text.
  3. [Table 6] The abbreviation 'allr.' is not defined; use 'allreduce' or a footnote. Also, the row label 'B0 = B0,k var.' conflicts with the notation m_update used elsewhere; align the notation.
  4. [Section 6.3.2] The 'generally impossible' conclusion should state the class of algorithms being excluded (for example, algorithms restricted to BLAS-level matvecs/solves and linear operations) and provide a formal proof or a reference. The reduction to computing (U^T B0 U)^{-1} x is convincing, but the impossibility claim is currently too informal for a journal.
  5. [Algorithm 15] In Algorithm 15, line 3 defines l = max{m-1, k-1} but l is never used; remove it.
  6. [Section 9.1] The inline command-line strings (e.g., '-mat lmvm J0 mat type diagonal ...') are embedded in prose and footnotes; set them off as code and explain their role in the text for readers not familiar with PETSc option syntax.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity: the intermediate dense form is an algebraic rearrangement of the independent compact-dense representation; performance claims are empirical.

full rationale

The paper's central derivation (Section 6.3) is not circular. Equation (6.17) is obtained by multiplying out blocks of the compact dense representation (6.8), which is attributed to the independent external result of Byrd, Nocedal, and Schnabel [6]. The claim that H_k^m v needs one B0^{-1} application and O(mn) BLAS work follows by counting operations in Algorithm 9, not by assuming the conclusion. The cached S^T g product reused across update/solve is an implementation optimization already present in the compact-dense literature; it is not fitting a parameter and renaming it as a prediction. Section 6.3.2's non-existence argument for an intermediate dense B.mult is a self-contained reduction to the impossibility of computing (U^T B0 U)^{-1}x with O(1) applications of B0, not an appeal to the paper's own prior claims. The benchmarks in Section 7.2 are empirical microbenchmarks; they involve no fitted quantities, and the paper explicitly discloses that the Tstep loop omits line search and objective evaluation, so the limited scope is a generality limitation rather than circular reasoning. Self-citations ([2], [7], [10], [14], [16]) are to the PETSc/TAO codebase and prior library work; none is used as the sole justification of the algebraic or performance conclusions. Accordingly, there are no circular steps to report.

Assumptions & free parameters 0 free parameters · 5 assumptions · 0 invented entities

The central derivation rests on the standard compact L-BFGS representation of Byrd/Nocedal/Schnabel and the BFGS recurrence; no free parameters are fit and no new entities are invented. The performance comparisons additionally assume the reported cost model and that the test H0 strategies represent the variable-metric setting.

assumptions (5)
  • standard math The compact dense representation (Eq. 6.8) of L-BFGS is taken as given from Byrd, Nocedal, and Schnabel [6].
    Section 6.2 builds the intermediate dense form by algebraic manipulation of Eq. (6.8); no re-derivation of the compact representation is given.
  • standard math The BFGS recurrence (Eqs. 2.1-2.5) is taken as background.
    Section 2 states the recurrence and the rank-2m expression without proof, as standard prior results.
  • domain assumption The cost model counting BLAS levels, memory traffic, and MPI collectives is assumed to capture performance.
    Table 6 and Section 7.1 compare algorithms on the basis of Level 1/2/3 BLAS operations, memops, and MPI synchronizations; kernel-launch overlap and other overheads are not modeled.
  • domain assumption The diagonal matrix with diagonal scaling is treated as representative of the variable-metric H0 case.
    Section 7.2 uses a full-memory restricted Broyden diagonal update [11] as the only variable-metric test; general variable-metric SPD B0 is not benchmarked.
  • domain assumption The non-existence proof in Section 6.3.2 assumes an 'analog' intermediate dense B.mult must use O(1) B0.mult/B0.solve applications for every SPD B0 and orthonormal U, and that exact computation of (U^T B0 U)^-1 x is impossible under that budget.
    This defines the scope of the negative result; it does not exclude dense B.mult variants that use O(m) base-matrix applications or relaxed exactness.

how reviews work

0 comments
Cite this review

Pith. "Pith review of A Portable and Versatile Limited-Memory BFGS Implementation in PETSc/TAO." pith.science (2026). https://pith.science/paper/YOT6F4IC

@misc{pith2026260800196,
  author       = {Pith},
  title        = {Pith review of: A Portable and Versatile Limited-Memory BFGS Implementation in PETSc/TAO},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/YOT6F4IC}},
  note         = {Machine review of arXiv:2608.00196}
}
abstract

The limited-memory BFGS (L-BFGS) Hessian update scheme is the critical kernel in many quasi-Newton optimization algorithms. The most common approach to implementing L-BFGS uses $2m$ sequential rank-1 updates as part of solving a linear system when there are $m$ history steps. The performance of this approach suffers when the latency of synchronization is significant, and its poor temporal locality increases the memory traffic when vectors do not fit in cache. The compact dense representation of L-BFGS results in an approach that has minimal synchronization latency and better temporal locality, but it requires an additional pass over the basis vectors and an additional basis that must be recomputed when the $B_0$ matrix changes as in variable-metric methods. In the Portable Extensible Toolkit for Scientific Computation and the Toolkit for Advanced Optimization (PETSc/TAO), we have implemented an intermediate dense formulation of BFGS that retains most of the good characteristics of both the recursive and compact dense approaches. We report single-node performance tests of these implementations on the U.S. Department of Energy's Polaris and Frontier machines, testing both GPU-based and CPU-based computations.

Figures

Figures reproduced from arXiv: 2608.00196 by the authors.

Figure 1
Figure 1. Cyclic ordering. Left: cyclic ordering of the history vectors [PITH_FULL_IMAGE:figures/full_fig_p006_1.png] view at source ↗
Figure 2
Figure 2. BFGS Solve comparison, constant diagonal H0 without scaling. Be is the effective bandwidth, as defined in (7.1). Polaris (EPYC 7543P), 1 MPI CPU process In [PITH_FULL_IMAGE:figures/full_fig_p017_2.png] view at source ↗
Figure 3
Figure 3. BFGS Solve comparison, constant diagonal H0 without scaling. Polaris (EPYC 7543P), 4 MPI CPU processes 10−5 10−4 10−3 10−2 10−1 100 0 5 10 15 20 25 102 103 n = 104 105 106 107 1 2BS 1 4BS Tstep (s) Be (GB/s) method recursive comp. dense dense history size m 5 10 20 50 [PITH_FULL_IMAGE:figures/full_fig_p018_3.png] view at source ↗
Figures from the paper (9 more)
Figure 4
Figure 4. Figure 4: BFGS Solve comparison, diagonal H0 with diagonal scaling. Polaris (EPYC 7543P), 1 MPI CPU process sizes, for both 1 and 4 MPI processes. This behavior is expected: as summarized in [PITH_FULL_IMAGE:figures/full_fig_p018_4.png]
Figure 5
Figure 5. Figure 5: BFGS Solve comparison, diagonal H0 with diagonal scaling. Polaris (EPYC 7543P), 4 MPI CPU processes introduces overhead that slows the intermediate dense method in these cases. Even so, for the largest problem size, n = 107 , the intermediate dense formulation remains …
Figure 6
Figure 6. Figure 6: BFGS Mult comparison, diagonal H0 without scaling, Polaris (EPYC 7543P), 4 MPI CPU processes 10−3 10−2 0 200 400 600 105 106 n = 107 1 2BS 1 4BS 1 2BS 1 4BS Tstep (s) Be (GB/s) method recursive comp. dense dense history size m 5 10 20 50 [PITH_FULL_IMAGE:figures/full_…
Figure 7
Figure 7. Figure 7: BFGS Solve comparison, constant diagonal H0 without scaling. Polaris A100, 1 GPU Similar to the CPU-based results, we observe the compact dense approach exhibit￾ing much worse performance when H0 is not constant. This is not surprising, because different architecture d…
Figure 8
Figure 8. Figure 8: BFGS Solve comparison, constant diagonal H0 without scaling. Frontier MI250X, 1 GPU 10−3 10−2 0 100 200 300 400 500 105 106 n = 107 1 2BS 1 4BS 1 2BS 1 4BS Tstep (s) Be (GB/s) method recursive comp. dense dense history size m 5 10 20 50 [PITH_FULL_IMAGE:figures/full_f…
Figure 9
Figure 9. Figure 9: BFGS Solve comparison, diagonal H0 with diagonal scaling. Polaris A100, 1 GPU relative to the recursive algorithm. We find that for variable H0 with diagonal scaling, the highest speedup factor of 4.73 is achieved, with n = 105 with history size m = 50, and the speedup…
Figure 10
Figure 10. Figure 10: BFGS Solve comparison, constant diagonal H0 without scaling. Polaris A100, 4 GPUs 10−3 10−2 10−1 0 500 1,000 1,500 2,000 106 n = 107 108 1 2BS 1 4BS Tstep (s) Be (GB/s) method recursive comp. dense dense history size m 5 10 20 50 [PITH_FULL_IMAGE:figures/full_fig_p02…
Figure 11
Figure 11. Figure 11: BFGS Solve comparison, diagonal H0 with diagonal scaling. Polaris A100, 4 GPUs algorithm performs slightly worse than the recursive one. For the Mult() timings in [PITH_FULL_IMAGE:figures/full_fig_p022_11.png]
Figure 12
Figure 12. Figure 12: BFGS Mult comparison, diagonal H0 without scaling. Polaris A100, 4 GPUs convolutional neural network on the MNIST image classification benchmark. We integrated PETSc/TAO’s BFGS solvers with PyTorch, enabling their use for training a neural network (NN). The MNIST trai…

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

16 extracted references · 2 canonical work pages

  1. [1]

    H. Anzt, Y. M. Tsai, A. Abdelfattah, T. Cojean, and J. Dongarra,Evaluating the performance of NVIDIA’s A100 Ampere GPU for sparse and batched computations, in 2020 IEEE/ACM Performance Modeling, Benchmarking and Simulation of High Performance Computer Systems (PMBS), IEEE, Nov. 2020, https://doi.org/10.1109/pmbs51919.2020. 00009

    arXiv 2020

  2. [2]

    Balay, S

    S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. M. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, J. F aibussowitsch, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan,...

    2023

  3. [3]

    Balay, S

    S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, E. M. Constantinescu, A. Dener, J. F aibussowitsch, W. D. Gropp, et al.,PETSc/TAO users manual revision 3.22, tech. report, Argonne National Laboratory (ANL), Argonne, IL (United States), 09 2024, https://doi.org/10.2172/2476320, https://www.osti.gov/biblio/ 2476320

    doi:10.2172/2476320 2024

  4. [4]

    Brown, V

    J. Brown, V. Barra, N. Beams, L. Ghaffari, M. Knepley, W. Moses, R. Shakeri, K. Stengel, J. L. Thompson, and J. Zhang,Performance portable solid mechanics via matrix-freep-multigrid, 2022, https://arxiv.org/abs/2204.01722

    arXiv 2022

  5. [5]

    Brown and P

    J. Brown and P. Brune,Low-rank quasi-Newton updates for robust Jacobian lagging in Newton-type methods, in International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering, 2013, pp. 2554–2565

    2013

  6. [6]

    R. H. Byrd, J. Nocedal, and R. B. Schnabel,Representations of quasi-newton matrices and their use in limited memory methods, Mathematical Programming, 63 (1994), pp. 129–156, https://doi.org/10.1007/bf01582063

    doi:10.1007/bf01582063 1994

  7. [7]

    Dener and T

    A. Dener and T. Munson,Accelerating limited-memory quasi-newton convergence for large- scale optimization, in Lecture Notes in Computer Science, Springer International Publishing, 2019, pp. 495–507, https://doi.org/10.1007/978-3-030-22744-9 39

    work page doi:10.1007/978-3-030-22744-9 2019

  8. [8]

    J. B. Erway, V. Jain, and R. F. Marcia,Shifted limited-memory DFP systems, in 2013 Asilomar Conference on Signals, Systems and Computers, IEEE, 2013, pp. 1033–1037

    2013

Show all 16 references
  1. [9]

    J. B. Erway and R. F. Marcia,On solving large-scale limited-memory quasi-newton equations, Linear Algebra and its Applications, 515 (2017), pp. 196–225, https://doi.org/10.1016/j.laa. 2016.11.003

    2017 doi

  2. [10]

    F aibussowitsch, M

    J. F aibussowitsch, M. F. Adams, R. T. Mills, S. Zampini, and J. Zhang,Safe, seamless, and scalable integration of asynchronous GPU streams in PETSc, 2023, https://arxiv.org/ abs/2306.17801

    2023 arXiv

  3. [11]

    J. C. Gilbert and C. Lemar ´echal,Some numerical experiments with variable-storage quasi- Newton algorithms, Mathematical Programming, 45 (1989), pp. 407–435. 26H. SUH, T. ISAAC, A. DENER, T. MUNSON, H. ZHANG, R. T. MILLS

    1989

  4. [12]

    Larsson Traff, W

    J. Larsson Traff, W. D. Gropp, and R. Thakur,Self-consistent MPI performance guidelines, IEEE Transactions on Parallel and Distributed Systems, 21 (2010), pp. 698–709, https: //doi.org/10.1109/tpds.2009.120, http://dx.doi.org/10.1109/TPDS.2009.120

    2010 doi

  5. [13]

    D. C. Liu and J. Nocedal,On the limited memory BFGS method for large scale optimization, Mathematical Programming, 45 (1989), pp. 503–528, https://doi.org/10.1007/bf01589116

    1989 doi

  6. [14]

    R. T. Mills, M. F. Adams, S. Balay, J. Brown, A. Dener, M. Knepley, S. E. Kruger, H. Morgan, T. Munson, K. Rupp, B. F. Smith, S. Zampini, H. Zhang, and J. Zhang, Toward performance-portable PETSc for GPU-based exascale systems, Parallel Computing, 108 (2021), p. 102831, https:...

    2021

  7. [15]

    Nocedal and S

    J. Nocedal and S. J. Wright,Numerical Optimization, Springer New York, 2006, https: //doi.org/10.1007/978-0-387-40065-5

    2006 doi

  8. [16]

    Zhang, J

    J. Zhang, J. Brown, S. Balay, J. F aibussowitsch, M. Knepley, O. Marin, R. T. Mills, T. Munson, B. F. Smith, and S. Zampini,The PetscSF scalable communication layer, IEEE Transactions on Parallel and Distributed Systems, 33 (2022), pp. 842–853, https: //doi.org/10.1109/tpds.20...

    2022

Pith tools

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