Autor: Jorge Sastre

A Charlier-Polynomial Approach to the Matrix Exponential

The matrix exponential is one of the most frequently evaluated matrix functions, and small improvements in its numerical computation can have a significant impact on large-scale scientific workloads. In a paper published in the International Journal of Computer Mathematics, J. M. Alonso, E. Defez, J. Ibáñez, J. Sastre, and A. Ibáñez introduce a new method based on Charlier polynomials in combination with the scaling-and-squaring technique.

The approximating matrix polynomial is evaluated through an efficient polynomial evaluation scheme that keeps the computational cost competitive with existing state-of-the-art codes. The authors perform a comprehensive numerical comparison across different matrix types. In their testbed, the variant expm_charlier_1 consistently appears among the most accurate and stable methods, while other variants and a reference code (expm_new) show higher error levels in a substantial fraction of the cases.

A distinctive aspect of this work is the use of a discrete orthogonal polynomial family in the matrix-exponential context. To our knowledge, discrete orthogonal polynomials, like the case of Charlier polynomials, have not been previously employed to address the matrix exponential computation.The article thus contributes a new polynomial family to the toolbox of matrix-exponential algorithms and provides empirical evidence of its numerical behaviour, without relying on asymptotic arguments alone.

Publication Details

  • Title: Computing Matrix Exponentials via Charlier Method
  • Authors: J. M. Alonso, E. Defez, J. Ibáñez, J. Sastre, A. Ibáñez
  • Journal: International Journal of Computer Mathematics, 1-20, 2026
  • DOI: 10.1080/00207160.2026.2697770

Main Contributions

  • A new algorithm for the matrix exponential based on Charlier matrix polynomialsand scaling-and-squaring.
  • Efficient evaluation of the approximating matrix polynomial via the Paterson–Stockmeyer method.
  • A numerical study showing that the expm_charlier_1 variant achieves competitive accuracy and stability across different matrix classes.
  • The first use, to the authors’ knowledge, of discrete orthogonal polynomials (Charlier) for matrix-exponential computation.

Why This Matters

Matrix exponentials appear in the solution of differential equations, in control theory, and in many scientific computing pipelines. Extending the set of reliable polynomial-based methods gives practitioners more options when selecting an algorithm for a specific matrix structure or accuracy requirement.

New Generating Functions for Chebyshev, Gegenbauer, and Jacobi Matrix Polynomials

Generating functions are a fundamental tool in the analysis of orthogonal polynomials, but their extension to the matrix setting has remained fragmented across different polynomial families. In a new article published in Mediterranean Journal of Mathematics, Emilio Defez and Michael M. Tung present a unified approach to deriving generating functions for the matrix analogues of the Chebyshev, Gegenbauer, and Jacobi polynomials.

The work does not stop at the matrix level. Alongside the matrix generating functions, the authors establish several explicit identities for the corresponding scalar polynomial families. These results provide new analytical tools for a class of polynomials that play a growing role in the numerical solution of matrix differential equations, the approximation of inverse Laplace transforms, and the high-accuracy computation of matrix functions such as the exponential.

By bringing three classical families under a single derivation framework, the article offers a compact reference for researchers working with matrix orthogonal polynomials and opens the door to analogous results for other polynomial sequences.

Publication Details

  • Title: New Generating Functions for Chebyshev, Gegenbauer, and Jacobi Matrix Polynomials
  • Authors: E. Defez, M. M. Tung
  • Journal: Mediterranean Journal of Mathematics, Vol. 23, 2026, Article 118
  • DOI: 10.1007/s00009-026-03119-x

Main Contributions

  • A unified derivation of generating functions for Chebyshev, Gegenbauer, and Jacobi matrix polynomials.
  • New explicit identities for the corresponding scalar polynomial families.
  • A consistent analytical framework that connects three classical polynomial sequences in the matrix setting.

Why This Matters

Matrix orthogonal polynomials are increasingly used in numerical linear algebra and in the approximation of matrix functions. Having a unified source of generating functions simplifies the derivation of new series expansions and can reduce the analytical effort required when adapting classical polynomial techniques to the matrix case.

Beyond Paterson–Stockmeyer: Advancing Matrix Polynomial Computation

For over fifty years, the Paterson–Stockmeyer method has been considered the benchmark for efficient matrix polynomial evaluation. In our recent open access article, we provide a summary of recent advances in this area and present a constructive scheme that evaluates a degree‑20 matrix polynomial using only 5 matrix multiplications—two fewer than Paterson–Stockmeyer.

We also show how the coefficients of this scheme can be derived from the solutions of a single equation involving one coefficient, and we include the full process in our supplementary materials.

Publication Details

  • Title: Beyond Paterson–Stockmeyer: Advancing Matrix Polynomial Computation
  • Authors: J. Sastre, J. Ibáñez, J. M. Alonso, E. Defez
  • Journal: WSEAS Transactions on Mathematics, Vol. 24, pp. 684–693, 2025
  • Conference: 5th Int. Conf. on Applied Mathematics, Computational Science and Systems Engineering (AMCSE), Paris, France, April 14–16, 2025
  • Open Access: https://doi.org/10.37394/23206.2025.24.68
  • Supplementary Material:

Main Contributions

  • Survey of recent advances in matrix polynomial evaluation.
  • Constructive result: A method to compute a degree‑20 matrix polynomial with just 5 matrix multiplications, improving efficiency over Paterson–Stockmeyer (needing 7 matrix products).
  • Coefficient derivation: All coefficients can be obtained by solving an equation in one unknown, documented step by step in the .txt file.
  • Generalization: We propose a framework for evaluation formulas of the type yk2(A)y_{k2}(A), see with Ck2C_k^2​ available variables, and set two conjectures for future research.

Why This Matters

Reducing matrix multiplications significantly lowers computational cost, which is crucial for:

  • Large-scale scientific computing
  • Numerical linear algebra
  • AI and machine learning models involving matrix functions

Access and Resources

Next Steps

If you work with matrix functions or large-scale computations:

  • Try the 5-multiplication scheme for degree‑20 polynomials.
  • Benchmark against Paterson–Stockmeyer.
  • Explore adapting the rational-coefficient approach to other degrees.

We welcome collaboration on proving the conjectures and extending these ideas to broader polynomial families.

Polynomial approximations for the matrix logarithm with computation graphs

Polynomial approximations for the matrix logarithm with computation graphs, E. Jarlebring, J. Sastre, J. Ibáñez, Linear Algebra Applications, in Press (open access), 2024. https://doi.org/10.1016/j.laa.2024.10.024, https://arxiv.org/abs/2401.10089, code.

In this article the matrix logarithm is computed by using matrix polynomial approximations evaluated by using matrix polynomial multiplications and additions. The most popular method for computing the matrix logarithm is a combination of the inverse scaling and squaring method in conjunction with a Padé approximation, sometimes accompanied by the Schur decomposition. The main computational effort lies in matrix-matrix multiplications and left matrix division. In this work we illustrate that the number of such operations can be substantially reduced, by using a graph based representation of an efficient polynomial evaluation scheme. A technique to analyze the rounding error is proposed, and backward error analysis is adapted. We provide substantial simulations illustrating competitiveness both in terms of computation time and rounding errors.

An Improved Taylor Algorithm for Computing the Matrix Logarithm

An Improved Taylor Algorithm for Computing the Matrix Logarithm, J. Ibáñez, J. Sastre, P. Ruiz, J.M. Alonso, E. Defez ), Mathematics. Vol. 9 (17), 2021, pp. 2018.

The matrix logarithm is the used in many applications of science and engineering [2], such as machine learning [7,8,9,10], computer-aided design (CAD) [19], computer graphics [17], the analysis of the topological distances between networks [23], graph theory [11,12], quantum chemistry and mechanics [3,4], buckling simulation [5], biomolecular dynamics [6], the study of Markov chains [13], sociology [14], optics [15], mechanics [16], control theory [18], optimization [20], the study of viscoelastic fluids [21,22], the study of brain–machine interfaces [24], and also in statistics and data analysis [25], among other areas.

The most popular method for computing the matrix logarithm is a combination of the inverse scaling and squaring method in conjunction with a Padé approximation, sometimes accompanied by the Schur decomposition. In this work, we present a Taylor series algorithm, based on the free-transformation approach of the inverse scaling and squaring technique, that uses recent matrix polynomial formulas for evaluating the Taylor approximation of the matrix logarithm more efficiently than the Paterson–Stockmeyer method. Two MATLAB implementations of this algorithm, related to relative forward or backward error analysis, were developed and compared with different state-of-the art MATLAB functions. Numerical tests showed that the new implementations are generally more accurate than the previously available codes, with an intermediate execution time among all the codes in comparison.

Efficient Evaluation of Matrix Polynomials beyond the Paterson–Stockmeyer Method

Efficient Evaluation of Matrix Polynomials beyond the Paterson–Stockmeyer Method (License CC BY 4.0), J. Sastre, J. Ibáñez, Mathematics 2021, 9(14), 1600; https://doi.org/10.3390/math9141600, Researchgate link

Two general methods for evaluating matrix polynomials requiring one matrix product less than the Paterson–Stockmeyer method were proposed in J. Sastre, Efficient evaluation of Matrix Polynomials, Linear Algebra Applications, where the cost of evaluating a matrix polynomial is given asymptotically by the total number of matrix product evaluations.

An analysis of the stability of those methods was given and the methods have been applied to Taylor-based implementations for computing the exponential, the cosine and the hyperbolic tangent matrix functions. Moreover, a particular example for the evaluation of the matrix exponential Taylor approximation of degree 15 requiring four matrix products was given, whereas the maximum polynomial degree available using Paterson–Stockmeyer method with four matrix products is 9.

Based on this example, a new family of methods for evaluating matrix polynomials more efficiently than the Paterson–Stockmeyer method was proposed, having the potential to achieve a much higher efficiency, i.e., requiring less matrix products for evaluating a matrix polynomial of certain degree, or increasing the available degree for the same cost. However, the difficulty of these family of methods lies in the calculation of the coefficients involved for the evaluation of general matrix polynomials and approximations.

In this paper, we provide a general matrix polynomial evaluation method for evaluating matrix polynomials requiring two matrix products less than the Paterson-Stockmeyer method for degrees higher than 30. Moreover, we provide general methods for evaluating matrix polynomial approximations of degrees 15 and 21 with four and five matrix product evaluations, respectively, whereas the maximum available degrees for the same cost with the Paterson–Stockmeyer method are 9 and 12, respectively. Finally, practical examples for evaluating Taylor approximations of the matrix cosine and the matrix logarithm accurately and efficiently with these new methods are given.

Advances in the Approximation of the Matrix Hyperbolic Tangent

Advances in the Approximation of the Matrix Hyperbolic Tangent (License CC BY 4.0), J. Ibáñez, J.M. Alonso, J. Sastre, E. Defez, P.A. Alonso, Mathematics 2021, 9(11), 1219; https://doi.org/10.3390/math9111219, Researchgate Link.

In this paper, we introduce two approaches to compute the matrix hyperbolic tangent. While one of them is based on its own definition and uses the matrix exponential, the other one is focused on the expansion of its Taylor series. For this second approximation, we analyse two different alternatives to evaluate the corresponding matrix polynomials. This resulted in three stable and accurate codes, which we implemented in MATLAB and numerically and computationally compared by means of a battery of tests composed of distinct state-of-the-art matrices.

Our results show that the Taylor series-based methods were more accurate, although somewhat more computationally expensive, compared with the approach based on the exponential matrix. To avoid this drawback, we propose the use of a set of formulas that allows us to evaluate polynomials in a more efficient way compared with that of the traditional Paterson–Stockmeyer method, thus, substantially reducing the number of matrix products (practically equal in number to the approach based on the matrix exponential), without penalising the accuracy of the result.

Simulation of harmonic oscillators on the lattice

Simulation of harmonic oscillators on the lattice, M.Tung, J. Ibáñez, E. Defez, J. Sastre, Mathematical Methods in the Applied Sciences 43(14), May 2020. Researchgate Link

This work deals with the simulation of a two‐dimensional ideal lattice having simple tetragonal geometry. The harmonic character of the oscillators give rise to a system of second‐order linear differential equations, which can be recast into matrix form. The explicit solutions which govern the dynamics of this system can be expressed in terms of matrix trigonometric functions.

For the derivation we employ the Lagrangian formalism to determine the correct solutions, which extremize the underlying action of the system. In the numerical evaluation we develop diverse state‐of‐the‐art algorithms which efficiently tackle equations with matrix sine and cosine functions. For this purpose, we introduce two special series related to trigonometric functions. They provide approximate solutions of the system through a suitable combination.

For the final computation an algorithm based on Taylor expansion with forward and backward error analysis for computing those series had to be devised. We also implement several MATLAB programs which simulate and visualize the two‐dimensional lattice and check its energy conservation.

Boosting the computation of the matrix exponential

Boosting the computation of the matrix exponential, J. Sastre, J. Ibáñez, E. Defez, Appl. Math. Comput. in press, 2018, doi:10.1016/j.amc.2018.08.017, Preprint, Matlab code expmpol.m.

This paper presents new Taylor algorithms for the computation of the matrix exponential based on recent new matrix polynomial evaluation methods. Those methods are more efficient than the well known Paterson–Stockmeyer method. The cost of the proposed algorithms is reduced with respect to previous algorithms based on Taylor approximations. Tests have been performed to compare the MATLAB implementations of the new algorithms to a state-of-the-art Padé algorithm for the computation of the matrix exponential, providing higher accuracy and cost performances.

First article with an application of the new matrix polynomial evaluation methods from J. Sastre, Efficient evaluation of matrix polynomials, Linear Algebra Appl. 539, (2018) 229-250. With the new matrix polynomial evaluation methods, Taylor approximation methods are more efficient than Padé approximant based methods.