pith. sign in

arxiv: 2602.17516 · v2 · submitted 2026-02-19 · 🧮 math.NA · cs.NA

Computing the action of a matrix exponential on an interval via the star-product approach

Stefano Pozza , Shazma Zahid This is my paper

Pith reviewed 2026-05-15 20:49 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords matrix exponentialaction on vectororthogonal polynomialsStein matrix equationstar algebranumerical linear algebrainterval evaluation
0
0 comments X

The pith

A ⋆-algebra representation converts computing e^{At}v over an interval into solving one Stein-type matrix equation for orthogonal polynomial coefficients.

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

The paper introduces a method to compute the action of the matrix exponential e^{At}v for every t inside a prescribed bounded interval. It expands the solution as an orthogonal polynomial series whose coefficients come from a single linear system. The system arises from a new representation of the matrix exponential inside the ⋆-algebra of bivariate distributions and takes the form of a Stein matrix equation. This equation is solved once, either directly or with Krylov methods, after which the series supplies the value at any desired t. A reader would care because the approach avoids recomputing the exponential separately for each time point.

Core claim

By placing the matrix exponential inside the ⋆-algebra, the problem of evaluating e^{At}v on an interval reduces to a Stein-type matrix equation whose solution supplies the coefficients of an orthogonal polynomial expansion that recovers the action for every t in the interval.

What carries the argument

The ⋆-product representation of the matrix exponential inside the algebra of bivariate distributions, which produces a Stein matrix equation solved for the polynomial coefficients.

If this is right

  • One solve of the Stein equation supplies coefficients that evaluate the action at any number of times inside the interval.
  • Both direct solvers and Krylov subspace methods can be used to obtain the coefficients.
  • The approach is intended to be competitive in accuracy and speed with existing methods for the same task.

Where Pith is reading between the lines

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

  • The same algebraic construction might extend to computing actions of other matrix functions if analogous representations can be found.
  • The single-interval solve could be reused inside larger time-stepping schemes that require the exponential at many successive intervals.
  • Because the output is an explicit polynomial series, it supplies not only point values but also derivatives with respect to t at negligible extra cost.

Load-bearing premise

The new representation inside the ⋆-algebra must produce coefficients that, when inserted into the orthogonal polynomial series, recover the true action e^{At}v for every t inside the chosen interval.

What would settle it

Compute the series for a small matrix whose exponential is known exactly, such as a diagonal matrix, and check whether the recovered values deviate from the true e^{At}v by more than the expected truncation error at several interior points of the interval.

Figures

Figures reproduced from arXiv: 2602.17516 by Shazma Zahid, Stefano Pozza.

Figure 1
Figure 1. Figure 1: Comparison of solution norms ke Atvk for Examples 4 and 8 over their respective time intervals ([0, 4] and [0, 0.1]). The ⋆-method results are shown as curves and expmv_tspan values as discrete points overlaid on that curve. The solution obtained by the ⋆-approaches can be used to compute the exponential at every point of the given interval [PITH_FULL_IMAGE:figures/full_fig_p017_1.png] view at source ↗
read the original abstract

We present a new method for computing the action of the matrix exponential on a vector, \( e^{At}v \). The proposed approach efficiently evaluates the solution for all \( t \) within a prescribed bounded interval by expanding it into an orthogonal polynomial series. This method is derived from a new representation of the matrix exponential in the so-called \(\star\)-algebra, an algebra of bivariate distributions. The resulting formulation leads to a linear system equivalent to a matrix equation of Stein type, which can be solved by either direct or Krylov subspace methods. Numerical experiments demonstrate the accuracy and efficiency of the proposed approach in comparison to state-of-the-art techniques.

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

2 major / 2 minor

Summary. The manuscript presents a new method for computing the action of the matrix exponential e^{At}v for all t in a prescribed bounded interval. The approach expands the solution in an orthogonal polynomial series obtained from a novel representation of the matrix exponential inside the ⋆-algebra of bivariate distributions; this representation is shown to yield a Stein-type matrix equation whose solution supplies the series coefficients, which can then be evaluated directly or via Krylov subspace methods. Numerical experiments are included to illustrate accuracy and efficiency relative to existing techniques.

Significance. If the algebraic equivalence between the ⋆-representation and the Stein-equation coefficients is rigorously verified, the method would supply a practical route to dense-in-time evaluations of the exponential action without repeated matrix-function calls. The combination of distribution-algebra techniques with orthogonal expansions and iterative solvers could be useful in applications such as time-dependent PDE discretizations or control problems where interval-wide information is required.

major comments (2)
  1. [Section 3 (derivation of the Stein equation from the ⋆-representation)] The central claim that the solution of the Stein-type equation recovers coefficients c_k such that the orthogonal polynomial series exactly equals e^{At}v (or with controllable error) for every t in the interval is asserted but not demonstrated. The mapping from the ⋆-exponential definition to the Stein equation must be shown to be free of residual terms arising from the support of the bivariate distributions or from non-commutativity of A; without this verification the numerical accuracy reported in the experiments cannot be guaranteed to hold in general.
  2. [Section 4 (numerical analysis and error discussion)] No a-priori error bounds or convergence rates for the polynomial truncation are supplied, nor is it shown how the truncation error interacts with the solution error of the Stein equation when the latter is solved iteratively. This omission makes it impossible to assess the method’s reliability for arbitrary interval lengths or matrix spectra.
minor comments (2)
  1. [Section 2] The notation for the ⋆-product and the precise definition of the bivariate-distribution algebra should be recalled explicitly in the main text rather than assumed from prior work.
  2. [Section 5] The numerical experiments section would benefit from a table summarizing CPU times, residual norms, and maximum pointwise errors across the tested interval for each competing method.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the careful reading of our manuscript and the constructive comments. We agree that additional rigor in the derivation and error analysis will improve the paper. We address each major comment below and indicate the planned revisions.

read point-by-point responses

  1. Referee: [Section 3 (derivation of the Stein equation from the ⋆-representation)] The central claim that the solution of the Stein-type equation recovers coefficients c_k such that the orthogonal polynomial series exactly equals e^{At}v (or with controllable error) for every t in the interval is asserted but not demonstrated. The mapping from the ⋆-exponential definition to the Stein equation must be shown to be free of residual terms arising from the support of the bivariate distributions or from non-commutativity of A; without this verification the numerical accuracy reported in the experiments cannot be guaranteed to hold in general.

    Authors: We thank the referee for highlighting this gap. The derivation in Section 3 proceeds from the ⋆-algebra representation by direct substitution into the Stein equation, using the definition of the ⋆-product on bivariate distributions. In the revised version we will insert a complete verification: we expand both sides explicitly, confirm that support conditions eliminate boundary residuals, and show that non-commutativity of A is handled by the algebraic structure without introducing extra terms. This establishes exact recovery of the coefficients (up to truncation) and justifies the reported numerical accuracy. revision: yes

  2. Referee: [Section 4 (numerical analysis and error discussion)] No a-priori error bounds or convergence rates for the polynomial truncation are supplied, nor is it shown how the truncation error interacts with the solution error of the Stein equation when the latter is solved iteratively. This omission makes it impossible to assess the method’s reliability for arbitrary interval lengths or matrix spectra.

    Authors: We agree that a-priori error analysis is missing. In the revision we will add a dedicated subsection deriving truncation-error bounds from the analyticity of the matrix exponential on the given interval and the approximation properties of the chosen orthogonal polynomials. We will also provide an explicit interaction estimate between the polynomial truncation error and the residual of the (iteratively solved) Stein equation, expressed in terms of interval length and the spectral radius of A. These bounds will allow quantitative assessment of reliability for general cases. revision: yes

Circularity Check

0 steps flagged

No circularity: derivation starts from new ⋆-algebra representation and proceeds to independent Stein equation

full rationale

The paper defines a new representation of the matrix exponential inside the ⋆-algebra of bivariate distributions, then derives a Stein-type linear system whose solution supplies the coefficients of the orthogonal polynomial expansion for e^{At}v over the interval. No quoted equation reduces the claimed recovery of the action to a fitted parameter, a self-citation chain, or a renaming of an existing result; the Stein equation is presented as a direct algebraic consequence of the ⋆-product definition rather than an input that is merely relabeled. Numerical experiments are offered as external verification, keeping the central claim self-contained against the stated assumptions.

Axiom & Free-Parameter Ledger

0 free parameters · 0 axioms · 0 invented entities

Based solely on the abstract, no explicit free parameters, axioms, or invented entities are identifiable; the ⋆-algebra is introduced as a new representation but its precise definition and any background assumptions are not stated.

pith-pipeline@v0.9.0 · 5405 in / 1276 out tokens · 38049 ms · 2026-05-15T20:49:23.614839+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

24 extracted references · 24 canonical work pages

  1. [1]

    H., and Higham, N

    A l-Mohy, A. H., and Higham, N. J. A new scaling and squaring algorithm for the matrix exponential. SIAM Journal on Matrix Analysis and Applications 31 , 3 (2010), 970–989

    work page 2010

  2. [2]

    H., and Higham, N

    A l-Mohy, A. H., and Higham, N. J. Computing the action of the matrix exponential, with a n application to exponential integrators. SIAM journal on scientific computing 33 , 2 (2011), 488–511. 18

    work page 2011

  3. [3]

    A prahamian, M., and Higham, N. J. The matrix unwinding function, with an application to computing the matrix exponential. SIAM Journal on Matrix Analysis and Applications 35 , 1 (2014), 88–109

    work page 2014

  4. [4]

    An exact formulation of the time- ordered exponential using path-sums

    G iscard, P .-L., Lui, K., T hw aite, S., and Jaksch, D. An exact formulation of the time- ordered exponential using path-sums. Journal of Mathematical Physics 56 , 5 (2015)

    work page 2015

  5. [5]

    Lanczos-like algorithm for the time-ordered exponenti al: the ⋆-inverse problem

    G iscard, P .-L.,and Pozza, S. Lanczos-like algorithm for the time-ordered exponenti al: the ⋆-inverse problem. Applications of Mathematics 65 , 6 (2020), 807–827

    work page 2020

  6. [6]

    A Lanczos-like method for non-autonomous linear ordina ry differential equations

    G iscard, P .-L.,and Pozza, S. A Lanczos-like method for non-autonomous linear ordina ry differential equations. Bollettino dell’Unione Matematica Italiana 16 , 1 (2023), 81–102

    work page 2023

  7. [7]

    H., and Meurant, G

    G olub, G. H., and Meurant, G. Matrices, moments and quadrature with applications . Princeton University Press, 2009

    work page 2009

  8. [8]

    H igham, N. J. Accuracy and stability of numerical algorithms , second ed. SIAM, Philadel- phia, P A, 2002

    work page 2002

  9. [9]

    H igham, N. J. Functions of matrices: theory and computation . SIAM Review, 2008

    work page 2008

  10. [10]

    Exponential integrators

    H ochbruck, M., and Ostermann, A. Exponential integrators. Acta Numerica 19 (2010), 209–286

    work page 2010

  11. [11]

    Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later

    M oler, C., and Van Loan, C. Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. SIAM Review 45, 1 (2003), 3–49

    work page 2003

  12. [12]

    A new closed-form expression for the solution of ODEs in a ring of distributions and its connection with the matrix algebra

    P ozza, S. A new closed-form expression for the solution of ODEs in a ring of distributions and its connection with the matrix algebra. Linear Multilinear Algebra 73, 9 (2025), 2025– 2035

    work page 2025

  13. [13]

    On the stability of network indices defined by means of mat rix functions

    P ozza, S., and Tudisco, F. On the stability of network indices defined by means of mat rix functions. SIAM J. Matrix Anal. Appl. 39 , 4 (2018), 1521–1546

    work page 2018

  14. [14]

    A ⋆-product solver with spectral accuracy for non- autonomous ordinary di fferential equations

    P ozza, S., and Van Buggenhout, N. A ⋆-product solver with spectral accuracy for non- autonomous ordinary di fferential equations. PAMM 23, 1 (2023), e202200050

    work page 2023

  15. [15]

    A new Legendre polynomial-based approach for non- autonomous linear ODEs

    P ozza, S., and Van Buggenhout, N. A new Legendre polynomial-based approach for non- autonomous linear ODEs. Electronic Transactions on Numerical Analysis 60 (2024), 292– 326

    work page 2024

  16. [16]

    A new Legendre polynomial approach for computing the mat rix exponential action on a vector

    P ozza, S., and Zahid, S. A new Legendre polynomial approach for computing the mat rix exponential action on a vector. PAMM 24, 4 (2024), e202400049

    work page 2024

  17. [17]

    A Frechet Lie group on distributions

    R yckebusch, M., B ouhamidi, A., and Giscard, P .-L. A Frechet Lie group on distributions. Journal of Mathematical Analysis and Applications 546 , 1 (2025), 129195

    work page 2025

  18. [18]

    Analysis of some Krylov subspace approximations to the m atrix exponential operator

    S aad, Y . Analysis of some Krylov subspace approximations to the m atrix exponential operator. SIAM Journal on Numerical Analysis 29 , 1 (1992), 209–228

    work page 1992

  19. [19]

    S idje, R. B. Expokit: A software package for computing matrix expo nentials. ACM Trans- actions on Mathematical Software 24 , 1 (1998), 130–156. 19

    work page 1998

  20. [20]

    Computational methods for linear matrix equations

    S imoncini, V . Computational methods for linear matrix equations. SIAM Review 58 , 3 (2016), 377–441

    work page 2016

  21. [21]

    Functions of di fference matrices are Toeplitz plus Hankel

    S trang, G., and MacNamara, S. Functions of di fference matrices are Toeplitz plus Hankel. SIAM Review 56, 3 (2014), 525–546

    work page 2014

  22. [22]

    Orthogonal polynomials, fourth ed., vol

    S zeg˝o, G. Orthogonal polynomials, fourth ed., vol. V ol. XXIII of American Mathematical Society Colloquium Publications . American Mathematical Society, Providence, RI, 1975

    work page 1975

  23. [23]

    N., W eideman, J

    T refethen, L. N., W eideman, J. A. C., and Schmelzer, T. Talbot quadratures and Rational approximations. BIT Numerical Mathematics 46 , 3 (2006), 653–670

    work page 2006

  24. [24]

    On the convergence rates of Legendre approximation

    W ang, H., and Xiang, S. On the convergence rates of Legendre approximation. Mathemat- ics of Computation 81 , 278 (2012), 861–877. 20

    work page 2012