In April 2016 Manchester eScholar was replaced by the University of Manchester’s new Research Information Management System, Pure. In the autumn the University’s research outputs will be available to search and browse via a new Research Portal. Until then the University’s full publication record can be accessed via a temporary portal and the old eScholar content is available to search and browse via this archive.

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

Al Mohy, Awad

[Thesis]. Manchester, UK: The University of Manchester; 2010.

Access to files

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.

Bibliographic metadata

Type of resource:
Content type:
Form of thesis:
Type of submission:
Degree type:
Doctor of Philosophy
Degree programme:
PhD Mathematical Sciences
Publication date:
Location:
Manchester, UK
Total pages:
155
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.
Thesis main supervisor(s):
Thesis co-supervisor(s):
Thesis advisor(s):
Language:
en

Institutional metadata

University researcher(s):

Record metadata

Manchester eScholar ID:
uk-ac-man-scw:85423
Created by:
Almohy, Awad
Created:
3rd July, 2010, 18:14:44
Last modified by:
Almohy, Awad
Last modified:
18th May, 2012, 18:49:35

Can we help?

The library chat service will be available from 11am-3pm Monday to Friday (excluding Bank Holidays). You can also email your enquiry to us.