Categoría: Matrix functions

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.

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.

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.

New matrix series expansions for the matrix cosine approximations

New matrix series expansions for the matrix cosine approximations, E. Defez, J. Ibáñez, P. Alonso, J.M. Alonso J. Peinado, and P. Alonso-Jordá, in the International Conference Mathematical Modelling in Engineering and Human Behaviour 2019,  Mathematical Modelling Conference Series at the Institute for Multidisciplinary Mathematics, Universitat Politècnica de València,  July, 10-12, 2019, Valencia (Spain).

Fast Taylor polynomial evaluation for the matrix cosine

Fast Taylor polynomial evaluation for the matrix cosine, J. Sastre, J. Ibañez, P. Alonso, J. Peinado and E. Defez, 18th International Conference on Computational and Mathematical Methods in Science and Engineering, CMMSE 2018, July 2018, Rota (Cadiz)-SPAIN.

In this work we introduce a new method to compute the matrix cosine. It is based on recent new matrix polynomial evaluation methods for the Taylor approximation and forward and backward error analysis. The matrix polynomial evaluation methods allow to evaluate the Taylor polynomial approximation of the cosine function more efficiently than using Paterson-Stockmeyer method. A MATLAB implementation of the new algorithm is provided, giving better efficiency and accuracy than state-of-the-art algorithms.