pith. sign in

arxiv: 2606.24701 · v1 · pith:D7NFNIXQnew · submitted 2026-06-23 · 🧮 math.NA · cs.NA

Recursive expansion of the matrix step function using polynomials of degree eight

Pith reviewed 2026-06-25 22:49 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords matrix step functionpolynomial expansionrecursive methodsmatrix polynomialsspectral gapmatrix computationsnumerical linear algebra
0
0 comments X

The pith

Recursive expansion using degree-eight polynomials reduces matrix multiplications for the matrix step function.

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

The paper introduces a recursive method to compute the matrix step function of large dense symmetric matrices by constructing a high-degree composite polynomial from repeated applications of degree-eight component polynomials. Each component is selected to strongly amplify the spectral gap across the step while positioning the updated gap for efficient subsequent iterations. A new evaluation scheme computes any such degree-eight polynomial using only three matrix-matrix multiplications and three matrices in memory. This expands the usable class of component polynomials and, together with the selection strategy, produces a consistent reduction in total matrix multiplications compared with earlier recursive approaches.

Core claim

The authors develop a recursive polynomial expansion for the matrix step function in which composite polynomials of high degree are assembled from component polynomials of degree exactly eight. The component polynomials are chosen to achieve strong amplification of the spectral gap while favorably repositioning that gap for the next iteration. The expansion relies on a novel evaluation scheme that realizes any degree-eight matrix polynomial with only three matrix-matrix multiplications and three stored matrices, thereby making available a substantially larger family of admissible component polynomials than was previously attainable within the same multiplication budget.

What carries the argument

The three-multiplication evaluation scheme for arbitrary degree-eight matrix polynomials, which enables a wider class of component polynomials for recursive gap-amplifying expansion of the matrix step function.

If this is right

  • The total number of matrix-matrix multiplications required to reach a given accuracy is reduced compared with prior recursive expansion techniques.
  • A substantially larger set of composite polynomials becomes available because the three-multiplication budget now accommodates any degree-eight component.
  • Each iteration produces stronger amplification of the spectral gap and better positioning for the subsequent iteration.
  • The method applies directly to any large dense symmetric matrix for which the matrix step function is required.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same three-multiplication building block could be reused inside expansions for other matrix functions that admit polynomial representations.
  • Because the scheme stores only three matrices at a time, memory usage per iteration remains bounded independently of the total degree.
  • If the initial spectral gap is too small, the number of iterations required may still dominate the overall cost even after the per-iteration saving.

Load-bearing premise

Component polynomials of degree eight can always be found and selected so that each iteration amplifies the spectral gap strongly while leaving the updated gap in a position that permits equally effective amplification in the following iteration.

What would settle it

A symmetric matrix with a known spectral gap for which the new recursive expansion requires at least as many matrix-matrix multiplications as any existing recursive method while still converging to the step function.

Figures

Figures reproduced from arXiv: 2606.24701 by Elias Jarlebring, Emanuel H. Rubensson, Gustaf Lorentzon.

Figure 1
Figure 1. Figure 1: Eigenvalue spectrum vs iteration in a polynomial recursion. The polynomials [PITH_FULL_IMAGE:figures/full_fig_p003_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: Illustration of the development from the McWeeny polynomial to accelerated [PITH_FULL_IMAGE:figures/full_fig_p006_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Illustration of the equioscillatory conditions and the resulting polynomial [PITH_FULL_IMAGE:figures/full_fig_p012_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: For each value of µ, L and R we define the polynomial p(L,R,µ−ε,µ+ε)(x) as the unique polynomial satisfying the oscillatory conditions in (4.4). For each construction, we have used the narrow gap ξ = 2ε, ε = 10−5 , which is sufficiently small to represent the limiting case ξ → 0. For each value of µ, the corresponding polynomials and their derivatives are evaluated at µ. In the figure legends, the label (L… view at source ↗
Figure 5
Figure 5. Figure 5: Number of matrix–matrix multiplications as a function of the step location [PITH_FULL_IMAGE:figures/full_fig_p017_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: Number of matrix–matrix multiplications (upper panel) and maximum error [PITH_FULL_IMAGE:figures/full_fig_p018_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: Ball-and-stick representation of the Aspirin molecule used in the numerical [PITH_FULL_IMAGE:figures/full_fig_p018_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: During the conditioning phase, the condition number is steadily reduced while [PITH_FULL_IMAGE:figures/full_fig_p019_8.png] view at source ↗
read the original abstract

We consider the problem of efficiently computing the matrix step function of a large dense symmetric matrix. To this end, we introduce a recursive polynomial expansion method in which a composite polynomial of high degree is built recursively from component polynomials of degree eight. The component polynomial used in each iteration is designed to achieve strong amplification of the spectral gap across the step while favorably positioning the updated gap for subsequent iterations. A key ingredient is a novel evaluation scheme for arbitrary matrix polynomials of degree exactly eight requiring only three matrix-matrix multiplications and three matrices in memory. This scheme makes available a substantially larger class of component polynomials than previously possible within a three-multiplication budget, thereby expanding the class of composite polynomials that can be generated. Together with our polynomial selection strategy, this leads to a significant and consistent reduction in the number of matrix-matrix multiplications required to compute the matrix step function compared to existing recursive expansion methods.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

0 major / 3 minor

Summary. The manuscript introduces a recursive polynomial expansion method for computing the matrix step function of large dense symmetric matrices. Composite high-degree polynomials are built from degree-eight component polynomials, each evaluated via a novel scheme using only three matrix-matrix multiplications and three auxiliary matrices. A polynomial selection strategy is used to amplify the spectral gap at each iteration while positioning the updated gap favorably for the next step. The central claim is that this yields a significant and consistent reduction in the total number of matrix-matrix multiplications relative to prior recursive expansion methods.

Significance. If the reported multiplication-count reductions hold under the same recursive structure and test conditions as the baselines, the work supplies a practical algorithmic advance for matrix-function evaluation in numerical linear algebra. Credit is due for the explicit three-multiplication evaluation scheme, the coefficient-optimization procedure, and the direct side-by-side multiplication counts against the methods cited in the abstract; these elements make the improvement verifiable rather than asserted.

minor comments (3)
  1. [Abstract] Abstract: the phrase 'significant and consistent reduction' would be more informative if accompanied by a brief quantitative statement (e.g., typical percentage or absolute count savings) or the class of test matrices employed.
  2. The notation for the composite polynomial and the updated spectral gap after each iteration could be introduced with a small worked example to improve readability for readers unfamiliar with the recursive construction.
  3. Section describing the degree-eight evaluation scheme: a short table listing the exact sequence of multiplications and additions (or a pseudocode block) would make the three-multiplication claim immediately verifiable without reconstructing the arithmetic from the text.

Simulated Author's Rebuttal

0 responses · 0 unresolved

We thank the referee for the positive summary, significance assessment, and recommendation of minor revision. No specific major comments were raised in the report, so we have no individual points requiring rebuttal or clarification at this stage. We remain available to incorporate any editorial or minor adjustments as needed for the final version.

Circularity Check

0 steps flagged

No significant circularity; derivation is algorithmic and self-contained

full rationale

The paper constructs an explicit degree-8 polynomial evaluation scheme (three matrix-matrix multiplications) and a recursive selection strategy for composite polynomials, then directly counts multiplications against prior recursive methods. No equations reduce a claimed prediction to a fitted input by construction, no load-bearing premise rests solely on self-citation, and no uniqueness theorem or ansatz is imported from the authors' prior work. The reported reduction follows from the new scheme's ability to access a larger class of component polynomials, which is independently verifiable from the stated coefficients and multiplication budget.

Axiom & Free-Parameter Ledger

1 free parameters · 1 axioms · 0 invented entities

Abstract-only review; no explicit free parameters, axioms, or invented entities are stated. The method implicitly assumes existence of suitable degree-eight polynomials within the three-multiplication constraint and that the spectral gap can be successively amplified.

free parameters (1)
  • component polynomial coefficients
    Chosen to amplify the spectral gap; selection strategy is mentioned but not detailed in the abstract.
axioms (1)
  • domain assumption The input matrix is large, dense, and symmetric.
    Stated as the problem setting in the abstract.

pith-pipeline@v0.9.1-grok · 5684 in / 1142 out tokens · 19902 ms · 2026-06-25T22:49:23.223178+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

47 extracted references · 39 canonical work pages · 1 internal anchor

  1. [1]

    J. M. Alonso, J. Sastre, J. Ibáñez, and E. Defez , A systematic framework for stable and cost-efficient matrix polynomial evaluation , 2026, https://doi.org/10.48550/arXiv.2603.23143, https://arxiv.org/abs/2603.23143, https://arxiv.org/abs/2603.23143

  2. [2]

    Amsel, Y

    N. Amsel, Y. Baumann, P. Beckman, P. B \"u rgisser, C. Cama \ n o, T. Chen, E. Chow, A. Damle, M. Derezinski, M. Embree, E. N. Epperly, R. Falgout, M. Fornace, A. Greenbaum, C. Greif, D. Halikias, Z. Huang, E. Jarlebring, Y. Koutis, D. Kressner, R. Kyng, J. Liesen, J. Lok, R. A. Meyer, Y. Nakatsukasa, K. Pearce, R. Peng, D. Persson, E. Rebrova, R. Schneid...

  3. [3]

    The Polar Express: Optimal Matrix Sign Methods and Their Application to the Muon Algorithm

    N. Amsel, D. Persson, C. Musco, and R. M. Gower , The polar express: Optimal matrix sign methods and their application to the muon algorithm , 2025, https://doi.org/10.48550/arXiv.2505.16932, https://arxiv.org/abs/2505.16932, https://arxiv.org/abs/2505.16932

  4. [4]

    Benzi, P

    M. Benzi, P. Boito, and N. Razouk , Decay properties of spectral projectors with applications to electronic structure , SIAM Rev., 55 (2013), pp. 3--64, https://doi.org/10.1137/100814019

  5. [5]

    D. R. Bowler and T. Miyazaki , O (N) methods in electronic structure calculations , Rep. Prog. Phys., 75 (2012), p. 036503, https://doi.org/10.1088/0034-4885/75/3/036503

  6. [6]

    M. J. Cawkwell, E. J. Sanville, S. M. Mniszewski, and A. M. N. Niklasson , Computing the density matrix in electronic structure theory on graphics processing units , J. Chem. Theory Comput., 8 (2012), pp. 4094--4101, https://doi.org/10.1021/ct300442w

  7. [7]

    J. H. Cheon, D. Kim, and D. Kim , Efficient homomorphic comparison methods with optimal complexity , in Advances in Cryptology -- ASIACRYPT 2020, S. Moriai and H. Wang, eds., Cham, 2020, Springer International Publishing, pp. 221--256, https://doi.org/https://doi.org/10.1007/978-3-030-64834-3_8

  8. [8]

    E. Chow, X. Liu, M. Smelyanskiy, and J. R. Hammond , Parallel scalability of hartree–fock calculations , J. Chem. Phys., 142 (2015), p. 104103, https://doi.org/10.1063/1.4913961

  9. [9]

    Dawson and T

    W. Dawson and T. Nakajima , Massively parallel sparse matrix function calculations with ntpoly , Comput. Phys. Commun., 225 (2018), pp. 154--165, https://doi.org/10.1016/j.cpc.2017.12.010

  10. [10]

    E. D. Denman and A. N. Beavers , The matrix sign function and computations in systems , Applied Mathematics and Computation, 2 (1976), pp. 63--94, https://doi.org/10.1016/0096-3003(76)90020-5

  11. [11]

    Finkelstein, J

    J. Finkelstein, J. S. Smith, S. M. Mniszewski, K. Barros, C. F. A. Negre, E. H. Rubensson, and A. M. N. Niklasson , Mixed precision Fermi -operator expansion on tensor cores from a machine learning perspective , J. Chem. Theory Comput., 17 (2021), pp. 2256--2265, https://doi.org/10.1021/acs.jctc.1c00057

  12. [12]

    Goedecker and L

    S. Goedecker and L. Colombo , Efficient linear scaling algorithm for tight-binding molecular dynamics , Phys. Rev. Lett., 73 (1994), pp. 122--125, https://doi.org/10.1103/PhysRevLett.73.122

  13. [13]

    Goedecker and M

    S. Goedecker and M. Teter , Tight-binding electronic-structure calculations and tight-binding molecular dynamics with localized orbitals , Phys. Rev. B, 51 (1995), pp. 9455--9464, https://doi.org/10.1103/PhysRevB.51.9455

  14. [14]

    Grishina, M

    E. Grishina, M. Smirnov, and M. Rakhuba , Accelerating newton-schulz iteration for orthogonalization via chebyshev-type polynomials , 2025, https://doi.org/10.48550/arXiv.2506.10935, https://arxiv.org/abs/2506.10935, https://arxiv.org/abs/2506.10935

  15. [15]

    N. J. Higham , Functions of Matrices: Theory and Computation , Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2008

  16. [16]

    Holas , Transforms for idempotency purification of density matrices in linear-scaling electronic-structure calculations , Chem

    A. Holas , Transforms for idempotency purification of density matrices in linear-scaling electronic-structure calculations , Chem. Phys. Lett., 340 (2001), pp. 552--558, https://doi.org/10.1016/S0009-2614(01)00409-2

  17. [17]

    Huss-Lederman, A

    S. Huss-Lederman, A. Tsao, and T. Turnbull , A parallelizable eigensolver for real diagonalizable matrices with real eigenvalues , SIAM J. Sci. Comput., 18 (1997), pp. 869--885, https://doi.org/10.1137/S1064827592228833

  18. [18]

    Jarlebring and G

    E. Jarlebring and G. Lorentzon , The polynomial set associated with a fixed number of matrix-matrix multiplications , SIAM Journal on Applied Algebra and Geometry, (2026), https://doi.org/10.48550/arXiv.2504.01500, https://arxiv.org/abs/2504.01500, https://arxiv.org/abs/2504.01500. to appear

  19. [19]

    Jordan, Y

    K. Jordan, Y. Jin, V. Boza, J. You, F. Cesista, L. Newhouse, and J. Bernstein , Muon: An optimizer for hidden layers in neural networks , 2024, https://kellerjordan.github.io/posts/muon/

  20. [20]

    Katbashev, R

    A. Katbashev, R. Schade, M. Lass, M. Müller, S. Grimme, A. Hansen, and T. D. Kühne , Submatrix and gpu-accelerated implementation of density matrix tight-binding , J. Chem. Phys., 163 (2025), p. 132501, https://doi.org/10.1063/5.0271379

  21. [21]

    Khadatkar and P

    S. Khadatkar and P. Motamarri , Subspace recursive fermi-operator expansion strategies for large-scale dft eigenvalue problems on hpc architectures , J. Chem. Phys., 159 (2023), p. 031102, https://doi.org/10.1063/5.0150287

  22. [22]

    Kim and Y

    J. Kim and Y. Jung , Accelerated purification using generalized nonpurifying intermediate functions for large-scale self-consistent field calculations , J. Chem. Theory Comput., 7 (2011), pp. 3853--3858, https://doi.org/10.1021/ct200441g. PMID: 26598332

  23. [23]

    Kruchinina, E

    A. Kruchinina, E. Rudberg, and E. H. Rubensson , Parameterless stopping criteria for recursive density matrix expansions , J. Chem. Theory Comput., 12 (2016), pp. 5788--5802, https://doi.org/10.1021/acs.jctc.6b00626

  24. [24]

    Kruchinina, E

    A. Kruchinina, E. Rudberg, and E. H. Rubensson , Efficient computation of the density matrix with error control on distributed computer systems , 2019, https://doi.org/10.48550/arXiv.1909.12533, https://arxiv.org/abs/1909.12533, https://arxiv.org/abs/1909.12533

  25. [25]

    Kulichenko, R

    M. Kulichenko, R. M. Stanton, C.-H. Li, N. Lubbers, S. Tretiak, M. E. Wall, C. F. A. Negre, J. Finkelstein, and A. M. N. Niklasson , Gpu-accelerated graph-based semiempirical quantum chemistry , J. Chem. Theory Comput., 21 (2025), pp. 10780--10792, https://doi.org/10.1021/acs.jctc.5c00936

  26. [26]

    Lee, J.-W

    E. Lee, J.-W. Lee, J.-S. No, and Y.-S. Kim , Minimax approximation of sign function by composite polynomial for homomorphic comparison , IEEE Transactions on Dependable and Secure Computing, 19 (2022), pp. 3711--3727, https://doi.org/10.1109/TDSC.2021.3105111

  27. [27]

    Liang, C

    W. Liang, C. Saravanan, Y. Shao, R. Baer, A. T. Bell, and M. Head-Gordon , Improved Fermi operator expansion methods for fast electronic structure calculations , J. Chem. Phys., 119 (2003), pp. 4117--4125, https://doi.org/10.1063/1.1590632

  28. [28]

    D. A. Mazziotti , Towards idempotent reduced density matrices via particle-hole duality: Mcweeny's purification and beyond , Phys. Rev. E, 68 (2003), p. 066701, https://doi.org/10.1103/PhysRevE.68.066701

  29. [29]

    McWeeny , The density matrix in self-consistent field theory

    R. McWeeny , The density matrix in self-consistent field theory. I. Iterative construction of the density matrix , Proc. R. Soc. London Ser. A, 235 (1956), pp. 496--509

  30. [30]

    A. M. N. Niklasson , Expansion algorithm for the density matrix , Phys. Rev. B, 66 (2002), p. 155115

  31. [31]

    A. M. N. Niklasson , Density Matrix Methods in Linear Scaling Electronic Structure Theory , Springer Netherlands, Dordrecht, 2011, pp. 439--473, https://doi.org/10.1007/978-90-481-2853-2_16

  32. [32]

    A. M. N. Niklasson, S. M. Mniszewski, C. F. A. Negre, M. J. Cawkwell, P. J. Swart, J. Mohd-Yusof, T. C. Germann, M. E. Wall, N. Bock, E. H. Rubensson, and H. Djidjev , Graph-based linear scaling electronic structure theory , J. Chem. Phys., 144 (2016), p. 234101, https://doi.org/10.1063/1.4952650

  33. [33]

    A. M. N. Niklasson, C. J. Tymczak, and M. Challacombe , Trace resetting density matrix purification in O (n) self-consistent-field theory , J. Chem. Phys., 118 (2003), pp. 8611--8620, https://doi.org/10.1063/1.1559913

  34. [34]

    A. H. R. Palser and D. E. Manolopoulos , Canonical purification of the density matrix in electronic-structure theory , Phys. Rev. B, 58 (1998), pp. 12704--12711, https://doi.org/10.1103/PhysRevB.58.12704

  35. [35]

    Pederson, J

    R. Pederson, J. Kozlowski, R. Song, J. Beall, M. Ganahl, M. Hauru, A. G. M. Lewis, Y. Yao, S. B. Mallick, V. Blum, and G. Vidal , Large scale quantum chemistry with tensor processing units , J. Chem. Theory Comput., 19 (2023), pp. 25--32, https://doi.org/10.1021/acs.jctc.2c00876

  36. [36]

    E. H. Rubensson , Nonmonotonic recursive polynomial expansions for linear scaling calculation of the density matrix , J. Chem. Theory Comput., 7 (2011), pp. 1233--1236

  37. [37]

    E. H. Rubensson and A. M. N. Niklasson , Interior eigenvalues from density matrix expansions in quantum mechanical molecular dynamics , SIAM J. Sci. Comput., 36 (2014), pp. B147--B170, https://doi.org/10.1137/130911585

  38. [38]

    E. H. Rubensson and E. Rudberg , Bringing about matrix sparsity in linear-scaling electronic structure calculations , J. Comput. Chem., 32 (2011), pp. 1411--1423, https://doi.org/10.1002/jcc.21723

  39. [39]

    E. H. Rubensson, E. Rudberg, and P. Sa ek , Density matrix purification with rigorous error control , J. Chem. Phys., 128 (2008), p. 074106, https://doi.org/10.1063/1.2826343

  40. [40]

    E. H. Rubensson, E. Rudberg, and P. Salek , Methods for Hartree-Fock and Density Functional Theory Electronic Structure Calculations with Linearly Scaling Processor Time and Memory Usage , Springer Netherlands, Dordrecht, 2011, pp. 263--300, https://doi.org/10.1007/978-90-481-2853-2_12

  41. [41]

    Rudberg, E

    E. Rudberg, E. H. Rubensson, P. Sałek, and A. Kruchinina , Ergo: An open-source program for linear-scaling electronic structure calculations , SoftwareX, 7 (2018), pp. 107--111, https://doi.org/10.1016/j.softx.2018.03.005

  42. [42]

    Rudberg, E

    E. Rudberg, E. H. Rubensson, P. Sałek, and A. Kruchinina , Ergo : a quantum chemistry program for large-scale self-consistent field calculations . https://ergoscf.org, 2026. Version 3.8.2, accessed 2026-05-28

  43. [43]

    Sastre , Efficient evaluation of matrix polynomials , Linear Algebra Appl., 539 (2018), pp

    J. Sastre , Efficient evaluation of matrix polynomials , Linear Algebra Appl., 539 (2018), pp. 229--250

  44. [44]

    Schulz , Iterative berechung der reziproken matrix , Zeitschrift für Angewandte Mathematik und Mechanik (ZAMM), 13 (1933), pp

    G. Schulz , Iterative berechung der reziproken matrix , Zeitschrift für Angewandte Mathematik und Mechanik (ZAMM), 13 (1933), pp. 57--59, https://doi.org/10.1002/zamm.19330130111

  45. [45]

    Squire and G

    W. Squire and G. Trapp , Using complex variables to estimate derivatives of real functions , SIAM Rev., 40 (1998), pp. 110--112, https://doi.org/10.1137/S003614459631241X

  46. [46]

    Suryanarayana , Optimized purification for density matrix calculation , Chem

    P. Suryanarayana , Optimized purification for density matrix calculation , Chem. Phys. Lett., 555 (2013), pp. 291--295, https://doi.org/10.1016/j.cplett.2012.10.090

  47. [47]

    VandeVondele, U

    J. VandeVondele, U. Bor s tnik, and J. Hutter , Linear scaling self-consistent field calculations with millions of atoms in the condensed phase , J. Chem. Theory Comput., 8 (2012), pp. 3565--3573, https://doi.org/10.1021/ct200897x