A New Scaling and Squaring Algorithm for the Matrix Exponential 论文

摘要

Abstract. The scaling and squaring method for the matrix exponential is based on the approximation eA ≈ (rm(2−sA)) 2s,whererm(x) isthe[m/m] Padéapproximant to ex and the integers m and s are to be chosen. Several authors have identified a weakness of existing scaling and squaring algorithms termed overscaling, in which a value of s much larger than necessary is chosen, causing a loss of accuracy in floating point arithmetic. Building on the scaling and squaring algorithm of Higham [SIAM J. Matrix Anal. Appl., 26 (2005), pp. 1179–1193], which is used by the MATLAB function expm, we derive a new algorithm that alleviates the overscaling problem. Two key ideas are employed. The first, specific to triangular matrices, is to compute the diagonal elements in the squaring phase as exponentials instead of from powers of rm. The second idea is to base the backward error analysis that underlies the algorithm on members of the sequence {‖Ak‖1/k} instead of ‖A‖, since for nonnormal matrices it is possible that ‖Ak‖1/k is much smaller than ‖A‖, andindeed this is likely when overscaling occurs in existing algorithms. The terms ‖Ak‖1/k are estimated without computing powers of A by using a matrix 1-norm estimator in conjunction with a bound of the form ‖Ak‖1/k ≤ max ( ‖Ap‖1/p, ‖Aq‖1/q) that holds for certain fixed p and q less than k. The improvements to the truncation error bounds have to be balanced by the potential for a large ‖A‖