Pith. sign in

REVIEW 2 major objections 5 minor 22 references

A MATLAB tool finds stable coefficient sets that cut two matrix multiplications from high-degree polynomial evaluation versus Paterson–Stockmeyer.

Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →

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

T0 review · grok-4.5

2026-07-31 12:13 UTC pith:IN6PPL4L

load-bearing objection Solid, usable MATLAB generator that turns a previously isolated 2M trick into a practical pipeline for many degrees, with the expected caveat that stable roots are discovered not guaranteed. the 2 major comments →

arxiv 2607.28286 v1 pith:IN6PPL4L submitted 2026-07-30 math.NA cs.NA

A MATLAB Tool for the Stable Generation of Matrix Polynomial Evaluation Schemes with Two-Product Savings

J. Ib\'a\~nez , J. Sastre , J.M. Alonso , E. Defez This is my paper
classification math.NA cs.NA MSC 15A1665F3065G5065Y15
keywords Matrix polynomialsPaterson–Stockmeyer methodMATLAB softwareNumerical stabilityFloating-point arithmeticMatrix functionsTwo-product savingsNonlinear systems
verification ladder T0 review T1 audit T2 compute T3 formal T4 reserved

The pith

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

Evaluating high-degree matrix polynomials is the bottleneck in many matrix-function algorithms, and the classic Paterson–Stockmeyer method has long set the multiplication count. This paper supplies a software procedure that systematically builds nested evaluation schemes saving two matrix products for degrees 18, 21, 24, 26–28 and every degree at least 30, provided the leading coefficient is nonzero. The coefficients come from reduced systems of nonlinear equations; the tool solves them in variable-precision arithmetic, then automatically filters for sets that remain accurate when rounded to single or double precision. Experiments on Taylor approximations of the exponential and logarithm, and on Euler-polynomial matrix exponentials, show that the resulting schemes match or beat the accuracy of ordinary Paterson–Stockmeyer while using two fewer multiplications. The practical payoff is a reusable MATLAB routine that lets any user test whether a given polynomial admits the cheaper evaluation path.

Core claim

For matrix polynomials of degrees m in {18, 21, 24, 26, 27, 28} and all m ≥ 30 with nonzero leading coefficient, the hybrid nested scheme Z2qs achieves computational cost exactly two matrix multiplications below optimal Paterson–Stockmeyer, and an automated stability filter can extract coefficient sets whose floating-point reconstruction error stays on the order of unit roundoff for the tested cases.

What carries the argument

The reduced s-by-s nonlinear system obtained by variable substitution from the 6s nested formulas (Y0, Y1, Y2), solved by randomized high-precision vpasolve and then filtered by the hybrid relative-absolute coefficient error ε against unit roundoff.

Load-bearing premise

A numerically stable real solution of the reduced nonlinear system must exist and be findable for the particular target coefficients; the paper itself states this is not guaranteed in advance for arbitrary polynomials.

What would settle it

Run MatrixPolEval2 on a concrete polynomial of degree ≥ 30 whose leading coefficient is nonzero; if every returned coefficient set produces reconstruction error ε many orders larger than unit roundoff in double precision, or if vpasolve returns no real solution, the claim of systematic stable 2M savings fails for that instance.

Watch this falsifier — get emailed when new claim-graph text bears on it.

If this is right

  • Any matrix-function code that evaluates Taylor or other polynomial approximants of degree ≥ 30 can attempt a two-multiplication saving by calling the supplied routine.
  • For the listed intermediate degrees the same tool also supplies alternative one-multiplication schemes when the earlier 1M generator fails.
  • Hybrid embedding of the 6s core inside a Paterson–Stockmeyer outer loop (Z2qs) can be more accurate than a pure 6s formula of slightly higher degree.
  • The reachable degree-versus-cost frontier for matrix polynomials is permanently lowered by two multiplications once a stable solution is confirmed.

Where Pith is reading between the lines

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

  • Because existence is only discovered per polynomial, libraries could cache successful coefficient sets for common series (exponential, logarithm, trigonometric) and fall back to Paterson–Stockmeyer otherwise.
  • The same reduction pipeline might be re-targeted at rational approximants or at simultaneous evaluation of several related polynomials, provided the underlying algebraic degree count still collapses.
  • Growth of single- and double-precision error for m ≥ 36 suggests that higher working precision or alternative nestings may be required before the method is routine for very high-degree approximants.

Editorial analysis

A structured set of objections, weighed in public.

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

Referee Report

2 major / 5 minor

Summary. The paper presents a MATLAB tool (MatrixPolEval2) that systematically generates evaluation schemes for matrix polynomials achieving a two-product (2M) savings over optimal Paterson–Stockmeyer for degrees m ∈ {18, 21, 24, 26, 27, 28} and all m ≥ 30 (nonzero leading coefficient). It implements the nested 6s structures of [9] via Algorithms 1–3, reduces the governing SNEs by the variable substitutions of Corollary 1, solves them with vpasolve under VPA, and filters candidate coefficient sets by a hybrid relative/absolute reconstruction error ε (Eq. 22). Theorem 2.1 establishes the hybrid cost CZ2qs(m) = CPS(m) − 2 whenever a stable solution exists and proves continuous degree coverage for m ≥ 30. Three case studies on Taylor exponential/logarithm coefficients, a 2×2 exponential example, and Euler-polynomial expm on 65 matrices from MCT/EMP show that, when stable coefficients are found, floating-point accuracy is preserved or improved relative to classic PS while the theoretical product count is realized.

Significance. If the generator reliably produces usable coefficients, the work supplies a practical, reproducible route to 2M savings that previously existed only for isolated Taylor cases. The public MATLAB implementation, explicit Algorithms 1–3, automated stability filter, and multiprecision reference comparisons (expm_mp) are concrete strengths that let other researchers test the schemes on new matrix functions. The contribution is correctly framed as a discovery-and-filter tool rather than an unconditional existence theorem; that framing is scientifically useful for the matrix-function community.

major comments (2)
  1. [Abstract, §1, Theorem 2.1, Tables 1–2, §4] The abstract and introduction claim a “generalized approach for arbitrary polynomials” and “systematically achieving the theoretical reduction of 2M.” Theorem 2.1 and the Conclusions correctly qualify this: a stable real root of the reduced s×s SNE is not guaranteed a priori and must be discovered per coefficient vector. Tables 1–2 already show that even for the well-conditioned Taylor exponential the double-precision reconstruction error jumps well above unit roundoff for m ≥ 36 (εdouble ≈ 3·10−13 and larger). The manuscript should therefore state the existence/stability caveat in the abstract and in the statement of the main claim, and should report, for a broader test suite of random or non-Taylor coefficient vectors, the fraction of instances in which vpasolve returns a coefficient set that passes the ε ≈ u filter. Without that empirical success rate the phrase “arbitrary polynomials
  2. [§3.3, Abstract] Case study 3 (the only large-scale accuracy experiment) evaluates only Euler-polynomial approximations to the matrix exponential. The abstract asserts that “numerical experiments involving various matrix functions” confirm accuracy preservation. Either additional matrix functions (logarithm, trigonometric, or general polynomials) should be added to the testbed, or the claim should be narrowed to the exponential setting actually studied. The present evidence base is too narrow to support the broader wording.
minor comments (5)
  1. [Figure 1, §2] Figure 1 caption and the surrounding text refer to MatrixPolEval1 achieving 1M savings for m = 8, 10 and all m ≥ 12; the body text earlier states that 1M is available only for a listed subset of degrees below 30. Align the two statements.
  2. [Listings 1–7] In Algorithm 2 listings the variable names y2, cs, n1, n2 are reused across nested loops without clear scoping comments; a short remark on how multiple real roots are stored would aid reproducibility.
  3. [Tables 1–2] Table 2 header uses “evpa(64)” and “esingle” while Table 1 uses “εvpa(64)” and “εsingle”; unify the notation.
  4. [§2.1–2.2] The hybrid scheme Z2qs (Eq. 23) is introduced after the pure 6s formulas; a forward reference in §2.1 would help the reader anticipate the general-degree case.
  5. [Throughout] Minor typographical inconsistencies appear (“F ramework”, “coefficients”, occasional missing spaces before citations). A final copy-edit pass is needed.

Circularity Check

0 steps flagged

No significant circularity: self-citations supply the known 6s nesting, but the tool’s claims are checked a posteriori by coefficient reconstruction and external multiprecision benchmarks.

full rationale

The paper’s load-bearing mathematical structures (the 6s nested formulas, the reduction of the full system to an s×s SNE, and the three Y1 variants) are explicitly adapted from the authors’ prior works [9], [8], [7]. That is ordinary cumulative research, not circular derivation: Proposition 1 and Corollary 1 are stated as adaptations, Theorem 2.1 is an independent cost-counting plus degree-coverage argument that holds whenever a stable coefficient set is supplied, and the stability metric ε is an a-posteriori filter that reconstructs the target polynomial coefficients and compares them to the input pk (relative or absolute). No parameter is fitted to data and then re-presented as a prediction; no uniqueness theorem is imported to forbid alternatives; experiments compare against external multiprecision references (expm_mp) and standard toolboxes. The paper itself caveats that stable real roots are not guaranteed a priori. Consequently the central software claim—generate-and-filter coefficient sets that achieve CPS(m)−2 when they exist—does not reduce to its inputs by construction. Score 1 reflects only the expected, non-load-bearing dependence on the authors’ earlier formulation of the nesting.

Axiom & Free-Parameter Ledger

2 free parameters · 4 axioms · 2 invented entities

The central claim rests on standard floating-point and matrix-polynomial arithmetic, on the algebraic reduction already proved in the authors’ earlier papers, and on the empirical hope that vpasolve will return a stable real root. No physical constants or data-fitted scales appear; the only free choices are solver precision, the three structural variants of Y1, and the stability threshold ε ~ u.

free parameters (2)
  • VPA working precision (64 decimal digits in experiments) = 64 decimal digits
    Chosen by the authors for coefficient generation; higher or lower precision can change which roots are found and how accurately they round.
  • Stability acceptance threshold ε = O(u) (u = 2^-24 or 2^-53)
    A solution is kept when reconstruction error is ‘of the order of unit roundoff’; the precise cutoff is a design choice that filters the reported coefficient sets.
axioms (4)
  • standard math Paterson–Stockmeyer cost formula CPS(m) = s + ceil(m/s) − 2 with optimal s ≈ √m
    Used as the baseline against which the 2M saving is measured (Section 2).
  • domain assumption The 6s nested evaluation formulas and the reduction of the coefficient-matching conditions to an s×s nonlinear system (Corollary 1 / [9])
    Taken as given from prior work; the present algorithms only implement the reduction.
  • domain assumption A nonzero leading coefficient p_{6s} ≠ 0 is sufficient for the variable substitutions that produce the reduced SNE
    Stated as condition (7); without it the whole construction collapses.
  • ad hoc to paper Existence of at least one numerically stable real root of the reduced SNE for the target coefficient vector
    Not proved; the tool discovers solutions when they exist. Conclusions explicitly warn that stable solutions are ‘not guaranteed a priori’.
invented entities (2)
  • MatrixPolEval2 coefficient-generation pipeline (Algorithms 1–3 + stability filter) independent evidence
    purpose: Automate construction and validation of 2M-saving evaluation schemes for arbitrary polynomials
    The concrete software artifact that turns earlier theoretical formulas into a usable generator.
  • Hybrid evaluation scheme Z2qs independent evidence
    purpose: Embed the 6s core inside a Paterson–Stockmeyer outer nesting so every degree m≥30 is covered with a 2M saving
    Introduced in equation (23) and Theorem 2.1 to remove degree gaps.

pith-pipeline@v1.2.0-daily-grok45 · 22902 in / 3160 out tokens · 50431 ms · 2026-07-31T12:13:25.628032+00:00 · methodology

0 comments
read the original abstract

Computing numerical approximations of matrix functions frequently relies on the efficient evaluation of high-degree matrix polynomials. Although computational bounds are historically governed by the Paterson--Stockmeyer (PS) method, recent theoretical developments have demonstrated the viability of evaluation schemes that eliminate two matrix products ($2M$). Existing literature documents stable instances of this $2M$ reduction only for isolated cases, such as specific degrees of Taylor approximations for the matrix exponential and the matrix logarithm. However, a generalized approach for arbitrary polynomials remains unestablished. To address this limitation, this work presents a software-driven procedure that extends these computational savings to polynomials of degrees $m \in \{18, 21, 24, 26, 27, 28\}$ and all $m \ge 30$, requiring primarily a non-zero leading coefficient. Since the underlying evaluation coefficients must be determined by solving systems of nonlinear equations (SNEs), selecting a numerically stable solution set is critical. We introduce an automated verification routine designed to filter and validate robust coefficient sets for floating-point execution. The primary contribution is a MATLAB implementation leveraging variable precision arithmetic to handle the underlying SNEs, verify stability, and project precision bounds. Numerical experiments involving various matrix functions verify that the developed implementation preserves or, in some instances, enhances the numerical accuracy of the classic PS method, while systematically achieving the theoretical reduction of $2M$.

Figures

Figures reproduced from arXiv: 2607.28286 by E. Defez, J. Ib\'a\~nez, J.M. Alonso, J. Sastre.

Figure 1
Figure 1. Figure 1: Reachable polynomial degrees m as a function of the computational cost in matrix￾matrix multiplications. The plot compares the classical Paterson–Stockmeyer (PS) base￾line, the 1M reduction scheme (MatrixPolEval1.m), the proposed 2M reduction scheme (MatrixPolEval2.m), and the specific scheme for m = 20 requiring 5 products [14]. For costs of 9 or more matrix multiplications, MatrixPolEval2.m achieves a sa… view at source ↗
Figure 2
Figure 2. Figure 2: Normwise relative errors for MATLAB functions [PITH_FULL_IMAGE:figures/full_fig_p019_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Performance profile (right) for MATLAB functions [PITH_FULL_IMAGE:figures/full_fig_p020_3.png] view at source ↗

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

22 extracted references · 3 canonical work pages

  1. [1]

    A. H. Al-Mohy, A truncated Taylor series algorithm for co mputing the action of trigonometric and hyperbolic matrix functions, S IAM Journal on Scientific Computing 40 (3) (2018) A1696–A1713

    2018

  2. [2]

    J. M. Alonso, J. Ibáñez, E. Defez, P. Alonso-Jordá, Euler polynomials for the matrix exponential approximation, Journal of Compu tational and Applied Mathematics 425 (2023) 115074

    2023

  3. [3]

    Caliari, F

    M. Caliari, F. Zivcovich, On-the-fly backward error esti mate for matrix exponential approximation by Taylor algorithm, Journal of Computational and Applied Mathematics 346 (2019) 532–548

    2019

  4. [4]

    Ibáñez, J

    J. Ibáñez, J. Sastre, P. Ruiz, J. M. Alonso, E. Defez, An im proved Taylor algorithm for computing the matrix logarithm, Mathematics 9 (17) (2021)

    2021

  5. [5]

    Bader, S

    P. Bader, S. Blanes, F. Casas, Computing the matrix expon ential with an optimized Taylor polynomial approximation, Mathematics 7 (12) (2019) 1174

    2019

  6. [6]

    M. S. Paterson, L. J. Stockmeyer, On the number of nonscal ar multiplica- tions necessary to evaluate polynomials, SIAM Journal on Co mputing 2 (1) (1973) 60–66

    1973

  7. [7]

    Sastre, Efficient evaluation of matrix polynomials, Li near Algebra and its Applications 539 (2018) 229–250

    J. Sastre, Efficient evaluation of matrix polynomials, Li near Algebra and its Applications 539 (2018) 229–250. 23

    2018

  8. [8]

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

    doi:10.48550/arxiv.2603.23143 2026

  9. [9]

    Sastre, J

    J. Sastre, J. Ibáñez, Efficient evaluation of matrix polyn omials beyond the Paterson-Stockmeyer method, Mathematics 9 (14) (2021) 160 0

    2021

  10. [10]

    Sastre, J

    J. Sastre, J. Ibáñez, E. Defez, Boosting the computatio n of the matrix exponential, Applied Mathematics and Computation 340 (201 9) 206–220

  11. [11]

    Seydaoğlu, P

    M. Seydaoğlu, P. Bader, S. Blanes, F. Casas, Computing t he matrix sine and cosine simultaneously with a reduced number of products , Applied Numerical Mathematics 163 (2021) 96–107

    2021

  12. [12]

    Bader, S

    P. Bader, S. Blanes, F. Casas, M. Seydaoğlu, An efficient a lgorithm to com- pute the exponential of skew-hermitian matrices for the tim e integration of the Schrödinger equation, Mathematics and Computers in Sim ulation 194 (2022) 383–400

    2022

  13. [13]

    Jarlebring, G

    E. Jarlebring, G. Lorentzon, The polynomial set associ ated with a fixed number of matrix-matrix multiplications (2025). arXiv:25 04.01500

    2025

  14. [14]

    Sastre, J

    J. Sastre, J. Ibánez, J. M. Alonso, E. Defez, Beyond pate rson–stockmeyer: Advancing matrix polynomial computation, WSEAS Transacti ons on Mathematics 24 (2025) 684–693. doi:10.37394/23206.2025. 24.68

    arXiv 2025

  15. [15]

    Fasi, Optimality of the paterson–stockmeyer method for evaluating ma- trix polynomials and rational matrix functions, Linear Alg ebra and its Ap- plications 574 (2019) 182–200

    M. Fasi, Optimality of the paterson–stockmeyer method for evaluating ma- trix polynomials and rational matrix functions, Linear Alg ebra and its Ap- plications 574 (2019) 182–200. doi:10.1016/j.laa.2019.0 4.001

    work page doi:10.1016/j.laa.2019.0 2019

  16. [16]

    A. H. Al-Mohy, N. J. Higham, A new scaling and squaring al gorithm for the matrix exponential, SIAM Journal on Matrix Analysis and Applications 31 (3) (2009) 970–989. doi:10.1137/09074721X

    doi:10.1137/09074721x 2009

  17. [17]

    N. J. Higham, The scaling and squaring method for the mat rix exponential revisited, SIAM Journal on Matrix Analysis and Application s 26 (4) (2005) 1179–1193

    2005

  18. [18]

    N. J. Higham, The Matrix Computation Toolbox, URL: http://www.ma.man.ac.uk/ higham/mctoolbox (2002)

    2002

  19. [19]

    T. G. Wright, Eigtool, version 2.1, URL: http://www.comlab.ox.ac.uk/pseudospectra/eigtool (20 09)

  20. [20]

    M. Fasi, N. J. Higham, Multiprecision algorithms for co mputing the ma- trix logarithm, SIAM Journal on Matrix Analysis and Applica tions 39 (1) (2018) 472–491. doi:10.1137/17M1129866

    work page doi:10.1137/17m1129866 2018

  21. [21]

    Advanpix, Multiprecision computing toolbox, URL: http://www.advanpix.com (2022). 24

    2022

  22. [2018]

    doi:10.3390/math9172018

    work page doi:10.3390/math9172018