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

Abstract

New algorithms for the matrix exponential and its Fr\'echetderivative are presented.First, we derive a new scaling and squaring algorithm (denoted$\mathrm{expm_{new}}$) for computing $e^A$, where $A$ is anysquare matrix, that mitigates the overscaling problem. Thealgorithm 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 diagonalelements in the squaring phase as exponentials instead ofpowering them. The second is to base the backward erroranalysis that underlies the algorithm on members of thesequence $\{\|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 usesthe scaling part of the scaling and squaring method togetherwith a truncated Taylor series approximation to theexponential. It determines the amount of scaling and the Taylordegree using the strategy of $\mathrm{expm_{new}}$.Preprocessing steps are used to reduce the cost of thealgorithm. An important application of the algorithm is toexponential integrators for ordinary differential equations. Itis shown that the sums of the form $\sum_{k=0}^p\varphi_k(A)u_k$ that arise in exponential integrators, wherethe $\varphi_k$ are related to the exponential function, can beexpressed in terms of a single exponential of a matrix ofdimension $n+p$ built by augmenting $A$ with additional rowsand columns.Third, a general framework for simultaneously computing amatrix function, $f(A)$, and its Fr\'echet derivative in thedirection $E$, $L_f(A,E)$, is established for a wide range ofmatrix functions. In particular, we extend the algorithm ofHigham and $\mathrm{expm_{new}}$ to two algorithms thatintertwine the evaluation of both $e^A$ and $L(A,E)$ at a costabout three times that for computing $e^A$ alone. These twoextended algorithms are then adapted to algorithms thatsimultaneously calculate $e^A$ together with an estimate of itscondition number.Finally, we show that $L_f(A,E)$, where $f$ is a real-valuedmatrix function and $A$ and $E$ are real matrices, can beapproximated by $\Im f(A+ihE)/h$ for some suitably small $h$.This approximation generalizes the complex step approximationknown in the scalar case, and is proved to be of second orderin $h$ for analytic functions $f$ and also for the matrix signfunction. It is shown that it does not suffer the inherentcancellation that limits the accuracy of finite differenceapproximations in floating point arithmetic. However,cancellation does nevertheless vitiate the approximation whenthe underlying method for evaluating $f$ employs complexarithmetic. The complex step approximation is attractive whenspecialized methods for evaluating the Fr\'echet derivative arenot available.

Date of Award	3 Jan 2011
Original language	English
Awarding Institution	The University of Manchester
Supervisor	Nicholas Higham (Supervisor) & Ronald Thatcher (Supervisor)