Algorithms for the Matrix Exponential and its Fr\'echet Derivative

Al-Mohy, Awad H. (2010) Algorithms for the Matrix Exponential and its Fr\'echet Derivative. Doctoral thesis, University of Manchester.

PDF almohy10.pdf Download (945kB)

Abstract

New algorithms for the matrix exponential and its Fr\'echet derivative are presented. First, we derive a new scaling and squaring algorithm (denoted $\mathrm{expm_{new}}$) for computing $e^A$, where $A$ is any square matrix, that mitigates the overscaling problem. The algorithm is built on the algorithm of Higham [{\em SIAM J. Matrix Anal. Appl.}, 26\penalty0(4):\penalty0 1179--1193, 2005] but improves on it by two key features. The first, specific to triangular matrices, is to compute the diagonal elements in the squaring phase as exponentials instead of powering them. The second is to base the backward error analysis that underlies the algorithm on members of the sequence $\{\|A^k\|^{1/k}\}$ instead of $\|A\|$. The terms $\|A^k\|^{1/k}$ are estimated without computing powers of $A$ by using a matrix 1-norm estimator. Second, a new algorithm is developed for computing the action of the matrix exponential on a matrix, $e^{tA}B$, where $A$ is an $n\times n$ matrix and $B$ is $n\times n_0$ with $n_0 \ll n$. The algorithm works for any $A$, its computational cost is dominated by the formation of products of $A$ with $n\times n_0$ matrices, and the only input parameter is a backward error tolerance. The algorithm can return a single matrix $e^{tA}B$ or a sequence $e^{t_kA}B$ on an equally spaced grid of points $t_k$. It uses the scaling part of the scaling and squaring method together with a truncated Taylor series approximation to the exponential. It determines the amount of scaling and the Taylor degree using the strategy of $\mathrm{expm_{new}}$. Preprocessing steps are used to reduce the cost of the algorithm. An important application of the algorithm is to exponential integrators for ordinary differential equations. It is shown that the sums of the form $\sum_{k=0}^p\varphi_k(A)u_k$ that arise in exponential integrators, where the $\varphi_k$ are related to the exponential function, can be expressed in terms of a single exponential of a matrix of dimension $n+p$ built by augmenting $A$ with additional rows and columns. Third, a general framework for simultaneously computing a matrix function, $f(A)$, and its Fr\'echet derivative in the direction $E$, $L_f(A,E)$, is established for a wide range of matrix functions. In particular, we extend the algorithm of Higham and $\mathrm{expm_{new}}$ to two algorithms that intertwine the evaluation of both $e^A$ and $L(A,E)$ at a cost about three times that for computing $e^A$ alone. These two extended algorithms are then adapted to algorithms that simultaneously calculate $e^A$ together with an estimate of its condition number. Finally, we show that $L_f(A,E)$, where $f$ is a real-valued matrix function and $A$ and $E$ are real matrices, can be approximated by $\Im f(A+ihE)/h$ for some suitably small $h$. This approximation generalizes the complex step approximation known in the scalar case, and is proved to be of second order in $h$ for analytic functions $f$ and also for the matrix sign function. It is shown that it does not suffer the inherent cancellation that limits the accuracy of finite difference approximations in floating point arithmetic. However, cancellation does nevertheless vitiate the approximation when the underlying method for evaluating $f$ employs complex arithmetic. The complex step approximation is attractive when specialized methods for evaluating the Fr\'echet derivative are not available.

Item Type:	Thesis (Doctoral)
Subjects:	MSC 2010, the AMS's Mathematics Subject Classification > 15 Linear and multilinear algebra; matrix theory MSC 2010, the AMS's Mathematics Subject Classification > 65 Numerical analysis
Depositing User:	Awad Al-Mohy
Date Deposited:	08 Jul 2010
Last Modified:	20 Oct 2017 14:12
URI:	https://eprints.maths.manchester.ac.uk/id/eprint/1494

Actions (login required)

View Item