Pith. sign in

REVIEW 3 major objections 5 minor 21 references

This paper claims that any matrix polynomial of degrees 18, 21, 24, 26, 27, 28, or at least 30 with a non-zero leading coefficient can be evaluated with two fewer matrix–matrix multiplications than the optimal Paterson–Stockmeyer method, wh

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 · deepseek-v4-flash

2026-08-03 01:34 UTC pith:IN6PPL4L

load-bearing objection Useful MATLAB tool with honest experiments, but the abstract's blanket 2M-savings claim is contradicted by the paper's own Table 1 for Taylor exponential coefficients. the 3 major comments →

arxiv 2607.28286 v2 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 methodtwo-product savingsmatrix functionsnumerical stabilitynonlinear systemsvariable precision arithmeticfloating-point arithmetic
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.

This paper sets out to show that a computational saving previously confined to a few special cases—evaluating a matrix polynomial with two fewer matrix–matrix multiplications than the optimal Paterson–Stockmeyer method—can be made systematic. Its central theorem claims the saving is guaranteed for degrees 18, 21, 24, 26, 27, 28, and for every degree at least 30, provided the nonlinear system that defines the evaluation coefficients has a numerically stable solution. The paper's concrete contribution is a MATLAB tool that solves that reduced system in variable precision arithmetic, filters candidate coefficient sets by a roundoff-error criterion, and outputs a scheme ready for floating-point use. This matters because matrix products dominate the cost of evaluating high-degree matrix polynomials, so two fewer products can translate into roughly a 20% reduction in practice; in the paper's experiments on 128-by-128 matrices the generated scheme used 455 matrix products versus 585 for Paterson–Stockmeyer while matching or improving accuracy. The claim is conditional: the paper states the stable-solution assumption explicitly, and its own tables show it failing for exponential Taylor coefficients at degrees 36 and 42 in double precision.

Core claim

The paper's central claim (Theorem 2.1) is that for matrix polynomials of degrees 18, 21, 24, 26, 27, 28, and every degree m≥30, there is an evaluation scheme whose cost is exactly two matrix–matrix multiplications less than the optimal Paterson–Stockmeyer scheme: C_{Z_{2q}^s}(m)=C_PS(m)−2. The reduction is achieved by a nested 6s-degree core (polynomials Y0, Y1, Y2 from formulas (10)–(14)) wrapped in a generalized PS hybrid that handles the residual q=m−6s degrees, giving total cost s+2+⌈q/s⌉. Equating the nested form to the target polynomial reduces, under the non-zero leading coefficient condition p_{6s}≠0, to an s×s nonlinear system; the coefficients are then whatever solution set is cho

What carries the argument

The load-bearing object is the 6s nested evaluation scheme built from three polynomials: Y0=A^s∑_{i=1}^s c_{5s+i}A^i, Y1 as one of three product structures (formulas (11)–(13)), and Y2=Y1(Y0+∑_{i=1}^s c_{s+i}A^i)+∑_{i=0}^s c_iA^i. When Y2 is forced to equal the target polynomial P_m of degree 6s, Corollary 1 reduces the 6s+1 equations to an s×s nonlinear system whose solutions determine all coefficients; the hybrid formula Z_{2q}^s embeds this 6s core in a PS-style loop to cover degrees m=6s+q. The stability filter is the error epsilon from (22): max relative (or absolute, if p_k=0) discrepancy between the target coefficients and the coefficients reconstructed from candidate coefficients rou

Load-bearing premise

The load-bearing premise, stated as an assumption in Theorem 2.1 and exposed by the paper's own Table 1, is that for every degree in the claimed range and every coefficient vector with non-zero leading coefficient, the reduced nonlinear system has at least one solution set that remains accurate when rounded to single or double precision; for Taylor exponential coefficients at degrees 36 and 42 this premise fails in double precision.

What would settle it

Take the degree-42 Taylor polynomial of the matrix exponential, call the tool to generate its evaluation coefficients, and compute the stability error epsilon from (22) in double precision; the paper's Table 1 reports epsilon=3.12e-12, thousands of times the unit roundoff, so this is a concrete case where no stable 2M scheme was obtained. A reader could verify that for this case the tool's own filter rejects or degrades the scheme, or, conversely, find a different solution branch that passes the filter—either outcome would sharpen or refute the claimed universality.

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

If this is right

  • For every degree in the claimed set, the cost formula s+2+⌈q/s⌉ applies, so a degree-30 polynomial can be evaluated in 7 products instead of PS's 9, and the gap remains exactly 2 for all higher degrees.
  • The prior isolated 2M examples (specific Taylor degrees for the exponential and logarithm) become instances of one general recipe; the same routine can be pointed at any coefficient vector with p_6s≠0.
  • For degrees 12–29 outside the 2M set, the tool still provides a 1M saving for nine degrees and no saving for 15, 16, 17, giving users a fallback when the 1M tool from the earlier work fails its accuracy check.
  • On 128-by-128 matrices, using the generated scheme for Euler-polynomial exponential codes reduced total matrix products from 585 to 455 versus the PS version, a 22.2% drop, while producing smaller relative errors than PS in 90.77% of the test matrices and than a scaling-and-squaring Padé code in 83.08%.

Where Pith is reading between the lines

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

  • Extension (not claimed by the paper): the same coefficient generator should work for Padé, Chebyshev, or other polynomial families with non-zero leading coefficient; a quick way to test is to feed those coefficient vectors into the tool and check that the stability filter passes.
  • Extension: because the solver is randomized, the degree-by-degree existence statement currently rests on search, not proof; a deterministic root-finding or continuation step would turn 'a stable solution was found' into 'a stable solution exists.'
  • Extension: the stability filter only measures coefficient reconstruction error, so in an application one would still need to combine the 2M-saving evaluation with backward error estimation for the specific matrix function; the tool is best seen as a coefficient oracle inside a larger algorithm.
  • Extension: the measured saving is asymptotic in matrix size; for small matrices the extra scalar preprocessing and coefficient generation cost may offset the two-product gain, so the practical payoff is concentrated in large-scale computations.

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

3 major / 5 minor

Summary. The paper presents a MATLAB tool, MatrixPolEval2, that generates evaluation coefficients for matrix polynomials of degree m=6s+q using a nested 6s structure, aiming to achieve a two-matrix-product (2M) saving over the Paterson–Stockmeyer (PS) method. The main theoretical result, Theorem 2.1, is a cost-counting proof that, if the associated system of nonlinear equations admits a numerically stable solution, then evaluation with the proposed scheme costs C_Z = s+2+ceil(q/s), and for m in {18,21,24,26,27,28} and all m≥30 this is two products fewer than optimal PS. The paper also introduces an automated stability filter based on reproducing target scalar coefficients in finite precision, and reports numerical experiments on Taylor approximations of the exponential and logarithm and on Euler-polynomial approximations of the matrix exponential.

Significance. If the unqualified claim were true, this would be a practically useful contribution: a public MATLAB tool that systematically reduces matrix polynomial evaluation cost by 2M relative to PS while preserving accuracy. The paper documents the algorithms and code, and the cost-counting part of Theorem 2.1 is straightforward and likely correct. However, the paper's own data and statements show that the existence of numerically stable coefficient sets is not guaranteed, and the abstract overstates the result. The significance is therefore conditional: the tool is a valuable exploratory instrument, but the claimed universal 2M saving for arbitrary polynomials of the listed degrees is not established.

major comments (3)
  1. [Abstract; Theorem 2.1 (Section 2.2); Table 1; Conclusion] The abstract claims that the savings extend to degrees m∈{18,21,24,26,27,28} and all m≥30 'requiring primarily a non-zero leading coefficient.' This is contradicted by the theorem itself, which explicitly assumes 'the algebraic system defining the 6s evaluation scheme admits numerically stable solutions.' The proof of Theorem 2.1 is only a product-counting argument; it does not prove existence. Table 1 shows that for Taylor exponential coefficients, which have p6s≠0, ε_double is 3.18e-14 at m=36 and 3.12e-12 at m=42, orders of magnitude above u≈1.11e-16, and ε_single is infinite at m=42. Thus the 'non-zero leading coefficient' condition is empirically false for degrees explicitly covered by the universal claim. The conclusion itself admits that 'the existence of numerically stable solutions is not guaranteed a priori.' The abstract must be revised to state the claim conditionally.
  2. [Section 2.2, Eq. (22); Algorithm 2 (Listing 4)] The automated 'stability verification' checks only whether the rounded coefficient set c reproduces the target scalar coefficients p to unit roundoff (internal reproduction). This is not a forward-stability certificate for evaluating the polynomial on a matrix. Moreover, the tool uses vpasolve's randomized search without any proof of convergence to a solution, let alone a stable one. The existence assumption in Theorem 2.1 therefore remains unverified for the general claim. To support the paper's central assertion, the authors need either a proof of existence of stable solutions for the claimed degrees, or a substantially broader empirical study over a wide class of coefficient vectors—not just the special Taylor/Euler cases shown.
  3. [Section 3] The numerical experiments are narrow: they cover Taylor approximations of the exponential and logarithm, and Euler polynomials for the exponential. The abstract and introduction promise a 'generalized approach for arbitrary polynomials,' but no tests are shown for arbitrary coefficient vectors, sparse or near-zero coefficients, or polynomials with large coefficient variation. Consequently, the experimental evidence does not support the claim that the 2M saving is systematically achievable for arbitrary polynomials. The paper should either restrict the claim to the tested classes or add experiments that sample a much broader spectrum of coefficient vectors.
minor comments (5)
  1. [Listing 5 (Algorithm 2)] The loop variable i in the nested loops is not used in the assignment; the intended operation appears to be cs{k1,k2}(i) = subs(...). This is likely a typographical or indexing error that should be corrected.
  2. [Eq. (18)] The notation \bar{y}_0 appears as a scalar in the expression for \bar{Y}_2, but from context it should refer to the polynomial \bar{Y}_0 (or its high-degree part). The inconsistency makes the formula hard to parse.
  3. [Section 2.2] The notation for polynomial coefficients is inconsistent: y0k, y1k, y2k are defined in (15), but later \bar{y}_0, \bar{y}_2, and related symbols are used without a clear mapping. Using a unified notation would improve readability.
  4. [Table 5] Only the total number of matrix products is reported for the testbed. To interpret the 22% reduction, the number of polynomials, their degrees, and the scaling/squaring parameters should be specified.
  5. [Section 2.1 and 2.2] The relationship between Algorithm 1, Algorithm 2, and Algorithm 3 is described only at a high level. A small worked example (e.g., s=3) would help readers verify the symbolic manipulation and the back-substitution steps.

Circularity Check

0 steps flagged

No circularity: the 2M saving is a conditional product-counting identity, and the stability filter is an internal round-trip check rather than a fitted prediction.

full rationale

Theorem 2.1's cost formula C_{Z_{2q}^s}(m)=s+2+ceil(q/s) is derived by counting matrix products in the nested evaluation (23); comparing it with the PS cost C_PS(m)=s+ceil(m/s)-2 is arithmetic, not a fitted prediction. The theorem is explicitly conditional ('Assume that the algebraic system defining the 6s evaluation scheme admits numerically stable solutions'), so the cost saving is not presented as an unconditional empirical claim. The stability check (22) reconstructs y2k from a rounded coefficient set and compares it with pk; this is an internal consistency check of the solve-and-round process, not a fitted parameter renamed as a prediction. The paper leans on earlier work by the same group ([8], [9], [10]) for the 6s structure and the reduction to s×s SNEs, but those are prior published mathematical results with stated assumptions (notably p6s≠0) and do not assume the present cost-saving conclusion; under the review rules they count as independent support rather than circular self-citation. The abstract's 'requiring primarily a non-zero leading coefficient' is broader than Theorem 2.1 and is contradicted for Taylor exponential coefficients at m=36 and 42 by the paper's own Table 1 (eps_double=3.18e-14 and 3.12e-12), but that is an overclaim/correctness risk, not a circular reduction: the paper's conclusion even admits 'the existence of numerically stable solutions is not guaranteed a priori for every target polynomial.' No equation or fitted parameter in the derivation reduces to its own input.

Axiom & Free-Parameter Ledger

0 free parameters · 4 axioms · 0 invented entities

No target-data fitting parameters are introduced. The coefficient sets are algebraic unknowns determined by the identity P_m=Y_2, not fitted to the numerical experiments. The main burden is the existence/stability assumption for the reduced system and reliance on same-group theoretical results [8,9,10].

axioms (4)
  • domain assumption Non-zero leading coefficient p_{6s} != 0 (condition (7))
    Needed to reduce the system from 6s+1 equations/7s+1 unknowns to an s×s system; polynomials with zero leading coefficient are excluded.
  • ad hoc to paper The reduced SNEs admit a numerically stable real solution for the target coefficients
    Theorem 2.1 assumes this. The 2M-saving guarantee and the practical stability of generated schemes depend on it; the paper does not prove it for arbitrary polynomials, and Table 1 shows failures at m=36 and m=42.
  • standard math PS cost formula C_PS(m)=s+ceil(m/s)-2 with optimal integer s
    Taken from [15] and used as the baseline for comparing costs.
  • domain assumption MATLAB vpasolve with random search finds a real solution when one exists
    The algorithm relies on numerical solver success; no formal convergence guarantee is given for the reduced nonlinear systems.

pith-pipeline@v1.3.0-alltime-deepseek · 16830 in / 14137 out tokens · 155211 ms · 2026-08-03T01:34:39.685008+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 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 2
Figure 2. Figure 2: Normwise relative errors for MATLAB functions [PITH_FULL_IMAGE:figures/full_fig_p021_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 ↗
Figure 3
Figure 3. Figure 3: Performance profile (right) for MATLAB functions [PITH_FULL_IMAGE:figures/full_fig_p021_3.png] view at source ↗

Review history (2 revisions) →

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

21 extracted references

  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) 2018

    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

    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

    arXiv 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. 25

    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

    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

    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

    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

    2018

  21. [21]

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

    2022