Efficient Evaluation of Matrix Polynomials beyond the Paterson-Stockmeyer Method

[EN] Recently, two general methods for evaluating matrix polynomials requiring one matrix product less than the Paterson-Stockmeyer method were proposed, 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.This research was partially funded by the European Regional Development Fund (ERDF) and the Spanish Ministerio de Economia y Competitividad grant TIN2017-89314-P, and by the Programa de Apoyo a la Investigacion y Desarrollo 2018 of the Universitat Politecnica de Valencia grant PAID-06-18-SP20180016.Sastre, J.; Ibáñez González, JJ. (2021). Efficient Evaluation of Matrix Polynomials beyond the Paterson-Stockmeyer Method. Mathematics. 9(14):1-23. https://doi.org/10.3390/math9141600S12391

Boosting the computation of the matrix exponential

[EN] 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 Pade algorithm for the computation of the matrix exponential, providing higher accuracy and cost performances.This work has been supported by Spanish Ministerio de Economia y Competitividad and European Regional Development Fund (ERDF) grant TIN2014-59294-P.Sastre, J.; Ibáñez González, JJ.; Defez Candel, E. (2019). Boosting the computation of the matrix exponential. Applied Mathematics and Computation. 340:206-220. https://doi.org/10.1016/j.amc.2018.08.017S20622034

Solving engineering models using hyperbolic matrix functions

In this paper a method for computing hyperbolic matrix functions based on Hermite matrix polynomial expansions is outlined. Hermite series truncation together with Paterson-Stockmeyer method allow to compute the hyperbolic matrix cosine efficiently. A theoretical estimate for the optimal value of its parameters is obtained. An efficient and highly-accurate Hermite algorithm and a MATLAB implementation have been developed. The MATLAB implementation has been compared with the MATLAB function funm on matrices of different dimensions, obtaining lower execution time and higher accuracy in most cases. To do this we used an NVIDIA Tesla K20 GPGPU card, the CUDA environment and MATLAB. With this implementation we get much better performance for large scale problems. (C) 2015 Elsevier Inc. All rights reserved.This work has been supported by Spanish Ministerio de Educacion TIN2014-59294-P.Defez Candel, E.; Sastre, J.; Ibáñez González, JJ.; Peinado Pinilla, J. (2016). Solving engineering models using hyperbolic matrix functions. Applied Mathematical Modelling. 40(4):2837-2844. https://doi.org/10.1016/j.apm.2015.09.050S2837284440

Solving initial value problems for ordinary differential equations by two approaches: BDF and Piecewise-linearized methods

Many scientific and engineering problems are described using Ordinary Differential Equations (ODEs), where the analytic solution is unknown. Much research has been done by the scientific community on developing numerical methods which can provide an approximate solution of the original ODE. In this work, two approaches have been considered based on BDF and Piecewise-linearized Methods. The approach based on BDF methods uses a Chord–Shamanskii iteration for computing the nonlinear system which is obtained when the BDF schema is used. Two approaches based on piecewise-linearized methods have also been considered. These approaches are based on a theorem proved in this paper which allows to compute the approximate solution at each time step by means of a block-oriented method based on diagonal Padé approximations. The difference between these implementations is in using or not using the scale and squaring technique. Five algorithms based on these approaches have been developed. MATLAB and Fortran versions of the above algorithms have been developed, comparing both precision and computational costs. BLAS and LAPACK libraries have been used in Fortran implementations. In order to compare in equality of conditions all implementations, algorithms with fixed step have been considered. Four of the five case studies analyzed come from biology and chemical kinetics stiff problems. Experimental results show the advantages of the proposed algorithms, especially when they are integrating stiff problems.Ibáñez González, JJ.; Hernández García, V.; Arias, E.; Ruíz Martínez, PA. (2009). Solving initial value problems for ordinary differential equations by two approaches: BDF and Piecewise-linearized methods. Computer Physics Communications. 180(5):712-723. doi:10.1016/j.cpc.2008.11.013S712723180

Approximating a Special Class of Linear Fourth-Order Ordinary Differential Problems

[EN] Differential matrix models are an important component of many interesting applications in science and engineering. This work elaborates a procedure to approximate the solutions of special non linear fourth-order matrix differential problems by suitable matrix splinesThis work has been supported by the Spanish Ministerio de Economía y Competitividad and the European Regional Development Fund (ERDF) under grant TIN2014-59294-PDefez Candel, E.; Tung, MM.; Ibáñez González, JJ.; Sastre, J. (2016). Approximating a Special Class of Linear Fourth-Order Ordinary Differential Problems. Springer. 577-584. https://doi.org/10.1007/978-3-319-63082-3_89S577584Defez, E., Tung, M.M., Ibáñez, J., Sastre, J.: Approximating and computing nonlinear matrix differential models. Math. Comput. Model. 55(7), 2012–2022 (2012)Famelis, I., Tsitouras, C.: On modifications of Runge–Kutta–Nyström methods for solving y (4) = f(x, y). Appl. Math. Comput. 273, 726–734 (2016)Golub, G.H., Loan, C.F.V.: Matrix Computations, 3rd edn. The Johns Hopkins University Press, Baltimore, MD (1996)Hussain, K., Ismail, F., Senu, N.: Two embedded pairs of Runge-Kutta type methods for direct solution of special fourth-order ordinary differential equations. Math. Probl. Eng. 2015 (2015). doi:10.1155/2015/196595Loscalzo, F.R., Talbot, T.D.: Spline function approximations for solutions of ordinary differential equations. SIAM J. Numer. Anal. 4(3), 433–445 (1967)Olabode, B., et al.: Implicit hybrid block Numerov-type method for the direct solution of fourth-order ordinary differential equations. Am. J. Comput. Appl. Math. 5(5), 129–139 (2015)Papakostas, S.N., Tsitmidelis, S., Tsitouras, C.: Evolutionary generation of 7th order Runge - Kutta - Nyström type methods for solving y (4) = f(x, y). In: American Institute of Physics Conference Series, vol. 1702 (2015). doi: 10.1063/1.493898

Approximating and computing nonlinear matrix differential models

NOTICE: this is the author’s version of a work that was accepted for publication in Mathematical and Computer Modelling. Changes resulting from the publishing process, such as peer review, editing, corrections, structural formatting, and other quality control mechanisms may not be reflected in this document. Changes may have been made to this work since it was submitted for publication. A definitive version was subsequently published in Mathematical and Computer Modelling Volume 55, Issues 7–8, April 2012, Pages 2012–2022 DOI: 10.1016/j.mcm.2011.11.060Differential matrix models are an essential ingredient of many important scientific and engineering applications. In this work, we propose a procedure to represent the solutions of first-order matrix differential equations Y(x) = f(x, Y(x)) with approximate matrix splines. For illustration of the method, we choose one scalar example, a simple vector model, and finally a Sylvester matrix differential equation as a test.This work has been supported by grant PAID-06-11-2020 from the Universitat Politecnica de Valencia, Spain.Defez Candel, E.; Tung ., MM.; Ibáñez González, JJ.; Sastre, J. (2012). Approximating and computing nonlinear matrix differential models. Mathematical and Computer Modelling. 55(7):2012-2022. https://doi.org/10.1016/j.mcm.2011.11.0602012202255

Accurate and efficient matrix exponential computation

[EN] This work gives a new formula for the forward relative error of matrix exponential Taylor approximation and proposes new bounds for it depending on the matrix size and the Taylor approximation order, providing a new efficient scaling and squaring Taylor algorithm for the matrix exponential. A Matlab version of the new algorithm is provided and compared with Pad´e state-of-the-art algorithms obtaining higher accuracy in the majority of tests at similar or even lower cost.This work has been supported by the Programa de Apoyo a la Investigacion y el Desarrollo of the Universitat Politecnica de Valencia grant PAID-06-11-2020Sastre, J.; Ibáñez González, JJ.; Ruiz Martínez, PA.; Defez Candel, E. (2014). Accurate and efficient matrix exponential computation. International Journal of Computer Mathematics. 91(1):97-112. https://doi.org/10.1080/00207160.2013.791392S97112911Al-Mohy, A. H., & Higham, N. J. (2010). A New Scaling and Squaring Algorithm for the Matrix Exponential. SIAM Journal on Matrix Analysis and Applications, 31(3), 970-989. doi:10.1137/09074721xArioli, M., Codenotti, B., & Fassino, C. (1996). The Padé method for computing the matrix exponential. Linear Algebra and its Applications, 240, 111-130. doi:10.1016/0024-3795(94)00190-1S. Blackford and J. Dongarra,Installation guide for LAPACK, LAPACK Working Note 411, Department of Computer Science, University of Tenessee, 1999.Dieci, L., & Papini, A. (2000). Padé approximation for the exponential of a block triangular matrix. Linear Algebra and its Applications, 308(1-3), 183-202. doi:10.1016/s0024-3795(00)00042-2Dieci, L., & Papini, A. (2001). Numerical Algorithms, 28(1/4), 137-150. doi:10.1023/a:1014071202885Dolan, E. D., & Moré, J. J. (2002). Benchmarking optimization software with performance profiles. Mathematical Programming, 91(2), 201-213. doi:10.1007/s101070100263C. Fassino,Computation of matrix functions, Ph.D. thesis TD-7/93, Università di Pisa, Genova, 1993.Higham, N. J. (2002). Accuracy and Stability of Numerical Algorithms. doi:10.1137/1.9780898718027Higham, N. J. (2005). The Scaling and Squaring Method for the Matrix Exponential Revisited. SIAM Journal on Matrix Analysis and Applications, 26(4), 1179-1193. doi:10.1137/04061101xHigham, N. J. (2008). Functions of Matrices. doi:10.1137/1.9780898717778Higham, N. J., & Tisseur, F. (2000). A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra. SIAM Journal on Matrix Analysis and Applications, 21(4), 1185-1201. doi:10.1137/s0895479899356080Moler, C., & Van Loan, C. (2003). Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later. SIAM Review, 45(1), 3-49. doi:10.1137/s00361445024180Paterson, M. S., & Stockmeyer, L. J. (1973). On the Number of Nonscalar Multiplications Necessary to Evaluate Polynomials. SIAM Journal on Computing, 2(1), 60-66. doi:10.1137/0202007Sastre, J., Ibáñez, J., Defez, E., & Ruiz, P. (2011). Accurate matrix exponential computation to solve coupled differential models in engineering. Mathematical and Computer Modelling, 54(7-8), 1835-1840. doi:10.1016/j.mcm.2010.12.04

Accurate Approximation of the Matrix Hyperbolic Cosine Using Bernoulli Polynomials

[EN] This paper presents three different alternatives to evaluate the matrix hyperbolic cosine using Bernoulli matrix polynomials, comparing them from the point of view of accuracy and computational complexity. The first two alternatives are derived from two different Bernoulli series expansions of the matrix hyperbolic cosine, while the third one is based on the approximation of the matrix exponential by means of Bernoulli matrix polynomials. We carry out an analysis of the absolute and relative forward errors incurred in the approximations, deriving corresponding suitable values for the matrix polynomial degree and the scaling factor to be used. Finally, we use a comprehensive matrix testbed to perform a thorough comparison of the alternative approximations, also taking into account other current state-of-the-art approaches. The most accurate and efficient options are identified as results.This research was supported by the Vicerrectorado de Investigacion de la Universitat Politecnica de Valencia (PAID-11-21).Alonso Abalos, JM.; Ibáñez González, JJ.; Defez Candel, E.; Alvarruiz Bermejo, F. (2023). Accurate Approximation of the Matrix Hyperbolic Cosine Using Bernoulli Polynomials. Mathematics. 11(3):1-22. https://doi.org/10.3390/math1103052012211

Numerical solutions of matrix differential models using higher-order matrix splines

The final publication is available at http://link.springer.com/article/10.1007%2Fs00009-011-0159-zThis paper deals with the construction of approximate solution of first-order matrix linear differential equations using higher-order matrix splines. An estimation of the approximation error, an algorithm for its implementation and some illustrative examples are included. © 2011 Springer Basel AG.Defez Candel, E.; Hervás Jorge, A.; Ibáñez González, JJ.; Tung, MM. (2012). Numerical solutions of matrix differential models using higher-order matrix splines. Mediterranean Journal of Mathematics. 9(4):865-882. doi:10.1007/s00009-011-0159-zS86588294Al-Said E.A., Noor M.A.: Cubic splines method for a system of third-order boundary value problems. Appl. Math. Comput. 142, 195–204 (2003)Ascher U., Mattheij R., Russell R.: Numerical solutions of boundary value problems for ordinary differential equations. Prentice Hall, New Jersey, USA (1988)Barnett S.: Matrices in Control Theory. Van Nostrand, Reinhold (1971)Blanes S., Casas F., Oteo J.A., Ros J.: Magnus and Fer expansion for matrix differential equations: the convergence problem. J. Phys. Appl. 31, 259–268 (1998)Boggs P.T.: The solution of nonlinear systems of equations by a-stable integration techniques. SIAM J. Numer. Anal. 8(4), 767–785 (1971)Defez E., Hervás A., Law A., Villanueva-Oller J., Villanueva R.: Matrixcubic splines for progressive transmission of images. J. Math. Imaging Vision 17(1), 41–53 (2002)Defez E., Soler L., Hervás A., Santamaría C.: Numerical solutions of matrix differential models using cubic matrix splines. Comput. Math. Appl. 50, 693–699 (2005)Defez E., Soler L., Hervás A., Tung M.M.: Numerical solutions of matrix differential models using cubic matrix splines II. Mathematical and Computer Modelling 46, 657–669 (2007)Mazzia F., Trigiante A.S., Trigiante A.S.: B-spline linear multistep methods and their conitinuous extensions. SIAM J. Numer. Anal. 44(5), 1954–1973 (2006)Faddeyev L.D.: The inverse problem in the quantum theory of scattering. J. Math. Physics 4(1), 72–104 (1963)Flett, T.M.: Differential Analysis. Cambridge University Press (1980)Golub G.H., Loan C.F.V.: Matrix Computations, second edn. The Johns Hopkins University Press, Baltimore, MD, USA (1989)Graham A.: Kronecker products and matrix calculus with applications. John Wiley & Sons, New York, USA (1981)Jódar L., Cortés J.C.: Rational matrix approximation with a priori error bounds for non-symmetric matrix riccati equations with analytic coefficients. IMA J. Numer. Anal. 18(4), 545–561 (1998)Jódar L., Cortés J.C., Morera J.L.: Construction and computation of variable coefficient sylvester differential problems. Computers Maths. Appl. 32(8), 41–50 (1996)Jódar, L., Ponsoda, E.: Continuous numerical solutions and error bounds for matrix differential equations. In: Int. Proc. First Int. Colloq. Num. Anal., pp. 73–88. VSP, Utrecht, The Netherlands (1993)Jódar L., Ponsoda E.: Non-autonomous riccati-type matrix differential equations: Existence interval, construction of continuous numerical solutions and error bounds. IMA J. Numer. Anal. 15(1), 61–74 (1995)Loscalzo F.R., Talbot T.D.: Spline function approximations for solutions of ordinary differential equations. SIAM J. Numer. Anal. 4(3), 433–445 (1967)Marzulli P.: Global error estimates for the standard parallel shooting method. J. Comput. Appl. Math. 34, 233–241 (1991)Micula G., Revnic A.: An implicit numerical spline method for systems for ode’s. Appl. Math. Comput. 111, 121–132 (2000)Reid, W.T.: Riccati Differential Equations. Academic Press (1972)Rektorys, K.: The method of discretization in time and partial differential equations. D. Reidel Pub. Co., Dordrecht (1982)Scott, M.: Invariant imbedding and its Applications to Ordinary Differential Equations. Addison-Wesley (1973

Computing Matrix Trigonometric Functions with GPUs through Matlab

[EN] This paper presents an implementation of one of the most up-to-day algorithms proposed to compute the matrix trigonometric functions sine and cosine. The method used is based on Taylor series approximations which intensively uses matrix multiplications. To accelerate matrix products, our application can use from one to four NVIDIA GPUs by using the NVIDIA cublas and cublasXt libraries. The application, implemented in C++, can be used from the Matlab command line thanks to the mex files provided. We experimentally assess our implementation in modern and very high-performance NVIDIA GPUs.This work has been supported by Spanish Ministerio de Economia y Competitividad and the European Regional Development Fund (ERDF) Grants TIN2014-59294-P and TEC2015-67387-C4-1-RAlonso-Jordá, P.; Peinado Pinilla, J.; Ibáñez González, JJ.; Sastre, J.; Defez Candel, E. (2019). Computing Matrix Trigonometric Functions with GPUs through Matlab. The Journal of Supercomputing. 75(3):1227-1240. https://doi.org/10.1007/s11227-018-2354-1S12271240753Serbin SM (1979) Rational approximations of trigonometric matrices with application to second-order systems of differential equations. Appl Math Comput 5(1):75–92Serbin Steven M, Blalock Sybil A (1980) An algorithm for computing the matrix cosine. SIAM J Sci Stat Comput 1(2):198–204Hargreaves GI, Higham NJ (2005) Efficient algorithms for the matrix cosine and sine. Numer Algorithms 40:383–400Al-Mohy Awad H, Higham Nicholas J (2009) A new scaling and squaring algorithm for the matrix exponential. SIAM J Matrix Anal Appl 31(3):970–989Defez E, Sastre J, Ibáñez Javier J, Ruiz Pedro A (2011) Computing matrix functions arising in engineering models with orthogonal matrix polynomials. Math Comput Model 57:1738–1743Sastre J, Ibáñez J, Ruiz P, Defez E (2013) Efficient computation of the matrix cosine. Appl Math Comput 219:7575–7585Al-Mohy Awad H, Higham Nicholas J, Relton Samuel D (2015) New algorithms for computing the matrix sine and cosine separately or simultaneously. SIAM J Sci Comput 37(1):A456–A487Alonso P, Ibáñez J, Sastre J, Peinado J, Defez E (2017) Efficient and accurate algorithms for computing matrix trigonometric functions. J Comput Appl Math 309(1):325–332CUBLAS library (2017) http://docs.nvidia.com/cuda/cublas/index.html . Accessed May 2017Alonso Jordá P, Boratto M, Peinado Pinilla J, Ibáñez González JJ, Sastre Martínez J (2014) On the evaluation of matrix polynomials using several GPGPUs. Universitat Politècnica de València, 2014. http://hdl.handle.net/10251/39615 . Accessed Sept 2017Boratto Murilo, Alonso Pedro, Giménez Domingo, Lastovetsky Alexey L (2017) Automatic tuning to performance modelling of matrix polynomials on multicore and multi-gpu systems. J Supercomput 73(1):227–239Alonso P, Peinado J, Ibáñez J, Sastre J, Defez E (2017) A fast implementation of matrix trigonometric functions sine and cosine. In: Proceedings of the 17th International Conference on Computational and Mathematical Methods in Science and Engineering (CMMSE 2017), pp 51–55, Costa Ballena, Rota, Cadiz (Spain), July 4th–8thSastre Jorge, Ibáñez Javier, Alonso Pedro, Peinado Jesús, Defez Emilio (2017) Two algorithms for computing the matrix cosine function. Appl Math Comput 312:66–77Paterson Michael S, Stockmeyer Larry J (1973) On the number of nonscalar multiplications necessary to evaluate polynomials. SIAM J Comput 2(1):60–66Higham Nicholas J (2008) Functions of matrices: theory and computation. SIAM, PhiladelphiaSastre J, Ibáñez Javier J, Defez E, Ruiz Pedro A (2011) Efficient orthogonal matrix polynomial based method for computing matrix exponential. Appl Math Comput 217:6451–6463Sastre J, Ibáñez Javier J, Defez E, Ruiz Pedro A (2015) Efficient scaling-squaring Taylor method for computing matrix exponential. SIAM J Sci Comput 37(1):A439–455Higham NJ, Tisseur F (2000) A block algorithm for matrix 1-norm estimation, with an application to 1-norm pseudospectra. SIAM J Matrix Anal Appl 21:1185–1201Demmel JW (1987) A counterexample for two conjectures about stability. IEEE Trans Autom Control 32:340–343Wright Thomas G (2002) EigTool library. http://www.comlab.ox.ac.uk/pseudospectra/eigtool/ . Accessed May 201