A NEW SCALING AND SQUARING ALGORITHM FOR THE MATRIX EXPONENTIAL
A NEW SCALING AND SQUARING ALGORITHM FOR THE MATRIX EXPONENTIAL
复制标题
DOI：
10.1137/09074721x
复制
发表时间：
2009-01-01
影响因子：
1.5
通讯作者：
Higham, Nicholas J.
中科院分区：
文献类型：
--
作者：
Al-Mohy, Awad H.;Higham, Nicholas J.
矩阵指数的缩放和平方方法基于近似E（a）近似为（r（m）（2（-s）a））（2s），其中r（m）（m）（x）是[m /m]淡淡的E（X）和整数M和S的幻象。几位作者已经确定了所谓的过度尺度的现有缩放算法和平方算法的弱点，其中选择了S的值大得多，从而导致浮点算术的准确性损失。建立在Higham的缩放和平方算法上[Siam J. Matrix肛门。 Appl。，26（2005），第1179-1193页]，由MATLAB函数expm使用，我们得出了一种减轻过度尺度问题的新算法。采用了两个关键思想。第一个特定于三角形矩阵的是在平方阶段计算对角线元素为指数，而不是从r（m）的幂中计算出对角线。第二个想法是基于基础的向后误差分析，该分析是基于序列成员{平行于（k）平行于（1/k）}的算法的，而不是平行于平行的算法，因为对于非正态矩阵，平行于平行于（1/k）的A（k）比平行与平行的平行得多，实际上，当在现有算法中发生过度尺度时，这很可能。通过使用矩阵1-norm估计量与平行于（k）平行于（1/k）平行的形式的形式结合使用矩阵1-norm估计量来估算平行于（1/k）的AK的术语，而无需计算A的功率。
The scaling and squaring method for the matrix exponential is based on the approximation e(A) approximate to (r(m)(2(-s) A))(2s), where r(m)(x) is the [m/m] Pade approximant to e(x) 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 r(m). The second idea is to base the backward error analysis that underlies the algorithm on members of the sequence {parallel to A(k)parallel to(1/k)} instead of parallel to A parallel to, since for nonnormal matrices it is possible that parallel to A(k)parallel to(1/k) is much smaller than parallel to A parallel to, and indeed this is likely when overscaling occurs in existing algorithms. The terms parallel to Ak parallel to(1/k) are estimated without computing powers of A by using a matrix 1-norm estimator in conjunction with a bound of the form parallel to A(k)parallel to(1/k)