# Jordan–Schur Algorithms for Computing the Matrix Exponential

Authors: Petko H. Petkov
Journal: International Journal of Mathematics and Mathematical Sciences
 (2023)

---

## Abstract

In this paper, two new versions of the Schur algorithm for computing the matrix exponential of ann×n complex matrix A are presented. Instead of the Schur form, these algorithms use the Jordan–Schur form of A. The Jordan–Schur form is found by less computation and it is determined more reliable than the reduction to Jordan form since it is obtained using only unitary similarity transformations. In contrast to the known methods, the diagonal blocks of the matrix exponential are obtained by using finite Taylor series. This improves the accuracy and avoids the decisions made about the termination of the series expansion. The off-diagonal blocks of the exponential are determined by modifications of the Schur–Parlett or Schur–Fréchet method, which takes advantage of the Jordan–Schur form of the matrix. The numerical features of the new algorithms are discussed, revealing their advantages and disadvantages in comparison with the other methods for computing the matrix exponential. Computational experiments show that using the new algorithms, the matrix exponential is determined in certain cases with higher accuracy than some widely used methods, however, at the price of an increase in the computational cost which is of order n4. It is shown that the Jordan–Schur algorithms for computing the matrix exponential are appropriate for matrices with multiple eigenvalues and are especially efficient in cases of large Weyr characteristics.


---

## Body

## 1. Introduction

The computation of the matrix exponentialeA of a real or complex n×n matrix A is an important task in the solution of linear differential equations, obtaining the time response and discretization of control systems, computing the integral of the exponential, and several other problems arising in applied mathematics. As it is well known [Moler and Van Loan [1, 2]], there exist at least nineteen “dubious” methods for computing the matrix exponential and there is no single “best” method for solving this important numerical problem. Among the good methods for computing the exponential that works in many cases, is the scaling and squaring method proposed in [Ward [3]], [Van Loan [4]] and improved in [Higham [5]]. An efficient method in this group combines the computation of the matrix exponential by using Padé approximations with scaling and squaring to improve the accuracy; see, for instance, [Higham [5]], [Al-Mohy and Higham [6]], and [Dieci and Papini [7]]. Alternative methods of the same group are based on the Taylor series, see [Caliari and Zivcovich [8]] and [Sastre et al. [9]].The methods to compute the matrix exponential of non-normal matrices that are based on the Schur form of the matrixA [Golub and Van Loan [10]; Ch. 9] and [Higham [11]; Ch. 10] are frequently preferred due to the numerical stability of the QR algorithm for computing the Schur form. Among the efficient methods for this aim is the Schur–Parlett algorithm [Davis and Higham [12]], which has comparable performance with the scaling and squaring method. The numerical properties of this algorithm are investigated in [Higham [13]; Ch. 9], see also Algorithm 10.23 in [Higham [11]; Ch. 10]. Another opportunity to compute the off-diagonal blocks of the matrix exponential by using the Schur form is provided by the Schur–Fréchet algorithm developed and analyzed in [Najfeld and Havel [14]], [Kenney and Laub [15]], and [Higham [13]; Ch. 9].The difficulty with the Schur–Parlett and Schur–Fréchet methods comes from the necessity to evaluate the exponential of a triangular matrix which may lead to large errors in some cases. Due to the occasional failures of these methods, there is a need for better versions that are more robust from the numerical point of view. In most cases, the failures occur in the presence of multiple eigenvalues ofA associated with complicated Jordan structures, so it is natural to make an attempt to use the Jordan form of the original matrix. However, determining the Jordan form of a matrix requires the use of nonunitary transformations which may cause numerical instability. More importantly, the computing of the Jordan form is an ill-posed computational problem since infinitely small perturbations of A may lead to abrupt changes in the Jordan structure, see [Golub and Wilkinson [16]] and [Demmel [17]]. That is why it seems advantageous to implement some other form of A which carries information about the Jordan structure of A but can be obtained by well-conditioned similarity transformations and is relatively insensitive to perturbations of the original matrix. Such form is the Jordan–Schur form of a square matrix, described in [Bai et al. [18]; Sect. 2.5], [Kågström and Wiberg [19]], and [Elmroth et al. [20]]. The Jordan–Schur can be considered as an intermediate form between Schur and Jordan form, and it is possible to compute it by a minor modification of a part of the algorithm of Kågström and Ruhe [Kågström and Ruhe [21, 22]] intended for finding the numerical Jordan form of a matrix. It should be noted that the reduction to the Jordan–Schur form is done by using only unitary similarity transformations and thus it is more reliable numerically than the reduction to the Jordan form. The use of the Jordan–Schur form allows for deriving more accurate algorithms for computing the matrix exponential, however at the price of an increase in the volume of computational work.In this paper, two new versions of the Schur method for computing the matrix exponential are presented. Instead of the Schur form, these algorithms use the Jordan–Schur form of a given complex matrix, which leads to some advantages in computing the exponential. In particular, the exponentials of the diagonal blocks are determined by using the finite Taylor series which improves the accuracy, and the off-diagonal blocks are found by modification of the Schur–Parlett or Schur–Fréchet method that takes advantage of the Jordan–Schur form of the matrix. It is confirmed by numerical examples that the new method may produce more accurate results than some widely used implementations of the scaling and squaring and Schur method.The paper is organized as follows. The Weyr and Segre characteristics, along with the Jordan–Schur form are introduced briefly in Section2. The new algorithms for computing the matrix exponential, based on the Jordan–Schur form are described in Section 3 and their numerical properties are discussed in Section 4. Four numerical experiments illustrating the behavior of the proposed algorithms are presented in Section 5, and some conclusions are drawn in Section 6.A notation will be the following.C is the set of complex numbers; AT is the transposed of A; AH=A¯T is the Hermitian transpose (the complex conjugate transpose) of A; In is the unit n×n matrix; Amij is the ij th block of the m th power of the block matrix A; A=Aiji=1:m,j=1:r is an m-by-r block matrix; λiA is the i th eigenvalue of A; A2, AF and A∞ are the spectral norm, the Frobenius norm, and the infinity norm of A, respectively; dimX is the dimension of the subspace X; nullA is the null space of A; rankA is the rank of A; ΔA is a perturbation of A; condFA is the condition number of A with respect to the inversion in the Frobenius norm; eA=expA is the exponential of A; condE is the condition number of the matrix exponential; A⊗B is the Kronecker product of A and B; sepA,B is the separation between square matrices A and B; condsylv is the condition number of the Sylvester equation.


## 2. The Jordan–Schur Form

The Jordan–Schur form is closely related to Jordan and Weyr canonical forms [Horn and Johnson [23], sect. 3.4]. To present this form, one needs to introduce the Segre and Weyr characteristics associated with multiple eigenvalues.


### 2.1. Segre and Weyr Characteristics

Consider the general case when a matrixA∈Cn×n has p≤n distinct eigenvalues λ1,λ2,…,λp with algebraic multiplicities.(1)ν1,ν2,…,νp,∑i=1pνi=n,respectively.Denote by,(2)Si=s1i,s2i,…,sqi,i,s1i≥s2i≥…≥sqi,i,∑j=1qisji=νi.and(3)Wi=wi1,wi2,…,wi,hi,wi1≥wi2≥⋯≥wi,hi,∑j=1hiwij=νi.The Segre and Weyr characteristics of A [Horn and Johnson [23]; p. 170], [Shapiro [24]], respectively, are associated with λi, where qi=wi1 and hi=s1i,i=1,2,…,p. The elements of the Segre characteristic Si are equal to the sizes of the Jordan blocks associated with λi in the Jordan form of A. The quantity qi,1≤qi≤νi is equal to the number of the linearly independent eigenvectors of A corresponding to λi, i.e., to the number of the Jordan blocks associated with λi. In turn, the elements of the Weyr characteristic Wi satisfy the relationships.(4)wij=dimnullA−λiInj−dimnullA−λiInj−1=rankA−λiInj−1−rankA−λiInj,j=1,2,…,hi.Note that,(5)h=maxihi=maxjs1j,is the index of nilpotency of A, i.e., the smallest positive integer for which.(6)dimnullA−λInhn.Letk1≥k2≥⋯≥0 be a partition of the positive integer k (i.e., k=k1+k2+⋯) and let mj denotes the number of the integers ki, which are greater or equal to j. Then, the integers mj form a partition of k, known as a conjugate partition of k. The relationship between the Segre and Weyr characteristics of A associated with an eigenvalue λi with multiplicity νi, is revealed by the fact that these characteristics are conjugated partitions of νi [Horn and Johnson [23]; p. 172]. Thus, the Weyr characteristic determines uniquely the Segre characteristic of A associated with λi and vice versa. Note that, the largest element s1i (the size of the largest Jordan block) of the Segre characteristic associated with λi is equal to the number hi of the elements of the corresponding Weyr characteristic. In turn, the first element wi1 of the Weyr characteristic (the geometric multiplicity of the eigenvalue λi) is equal to the number qi of the elements of the Segre characteristic (the number of Jordan blocks) associated with λi.


### 2.2. Reduction to Jordan–Schur Form

Using unitary similarity transformations with a matrixV, the matrix A∈Cn×n is reduced to the upper block triangular form (the Jordan–Schur form).(7)JS=VHAV=SW1,λ1T˜12…T˜1pSW2,λ2…T˜2p⋱⋮SWp,λp,where the diagonal blocks(8)SWi,λi=λiIwi1N12i…N1,hi−1iN1,hiiλiIwi2…N2,hi−1iN2,hii⋱⋮⋮λiIwi,hi−1Nhi−1,hiiλiIwi,hi∈ℂνi×νi,i=1,2,…,p,are in the so-called staircase form and the off-diagonal blocks T˜ij∈ℂwi×wj have no special form. Thus, the structure of the Jordan–Schur form is determined by the Weyr characteristics associated with λi,i=1,2,…,p. These characteristics can be used to find easily the corresponding Segre characteristics of A.A possible way to determine the Jordan–Schur form of a matrix is to implement a modification of the algorithm of Kågström and Ruhe [Kågström and Ruhe [21, 22]] intended to find the Jordan form of a matrix. This algorithm is based on the idea that using a sequence of unitary transformations it is possible to find the dimensions of the null spaces nullA−λIj (and hence the elements of the Weyr characteristic) without computing the matrix powers A−λIj, and consequently, the size of the Jordan blocks associated with multiple eigenvalues. As a result, the matrix is reduced to a staircase form and the corresponding algorithm is called a staircase algorithm because of the shape of the resulting matrix (for a detailed history of this algorithm, see [Edelman and Ma [25]]).The reduction of ann-th order complex matrix A to Jordan–Schur form can be carried out by a minor modification of the algorithm of Kågström and Ruhe [Petkov [26]].Concerning the computational cost, the reduction to the Jordan–Schur form requires more operations than the reduction to the Schur form. This is mainly due to the reduction of the diagonal blocks to the staircase form, which needs a large volume of operations if these blocks have a high index of nilpotence, i.e., if their Jordan forms involve high-order blocks. According to the estimates given in [Kågström and Ruhe [22]], in addition to the reduction into Schur form, the sorting of the eigenvalues, the clustering of multiple eigenvalues, and the reduction of the diagonal blocks of the matrix A to staircase form require 6n3+21/3n3hflops where one flop consists of one addition and one multiplication of floating point numbers. (note that, one flop in complex arithmetic is equivalent to four ordinary flops and the corresponding number of operations should be multiplied by four.) For h=n (nonderogatory but defective matrix with one Jordan block), this volume of work is of order 7/3n4+6n3flops and for h=1 (derogatory but nondefective matrix with scalar Jordan blocks) it is approximately 81/3n3flops. This way, for matrices with large Jordan blocks, the reduction to Jordan–Schur form requires a volume of work of order n4. Some increase of the efficiency of the reduction to staircase form is proposed, for instance, in [Guglielmi et al. [27]].Example 1.
Consider the matrix(9)A=1−57−1090−1−12−2−36−1012−110−23−652−57−108.
The following matrix has the Jordan form:(10)J=diag−11−1,−21−2−2.
The eigenvalueλ1=−1 has a Weyr characteristic equal to 1,1 and Segre characteristic equal to (2). The eigenvalue λ2=−2 has a Weyr characteristic equal to 2,1 and the same Segre characteristic.
Using a unitary similarity transformation with a matrix (all results are rounded to eight decimal digits),(11)V=0.500000000.56694671−0.18273332−0.628633410.00000000−0.166666670.566946710.754888190.15931703−0.23570226−0.666666670.000000000.11922195−0.564906590.471404520.16666667−0.566946710.42117723−0.50117978−0.471404520.50000000−0.188982240.452932910.095590220.70710678.
The matrixA is reduced by the modification of the algorithm of Kågström and Ruhe to the Jordan–Schur form.(12)JS=VHAV=−14.53557368−6.29342726−0.4431080930.05203820−1−0.84090739−0.95805781−0.13363062−203.06799138−2−1.85595575−2,
The reduction is carried out with a relative backward error.(13)A−VJSVHFAF=6.9⋅10−16.


## 2.1. Segre and Weyr Characteristics

Consider the general case when a matrixA∈Cn×n has p≤n distinct eigenvalues λ1,λ2,…,λp with algebraic multiplicities.(1)ν1,ν2,…,νp,∑i=1pνi=n,respectively.Denote by,(2)Si=s1i,s2i,…,sqi,i,s1i≥s2i≥…≥sqi,i,∑j=1qisji=νi.and(3)Wi=wi1,wi2,…,wi,hi,wi1≥wi2≥⋯≥wi,hi,∑j=1hiwij=νi.The Segre and Weyr characteristics of A [Horn and Johnson [23]; p. 170], [Shapiro [24]], respectively, are associated with λi, where qi=wi1 and hi=s1i,i=1,2,…,p. The elements of the Segre characteristic Si are equal to the sizes of the Jordan blocks associated with λi in the Jordan form of A. The quantity qi,1≤qi≤νi is equal to the number of the linearly independent eigenvectors of A corresponding to λi, i.e., to the number of the Jordan blocks associated with λi. In turn, the elements of the Weyr characteristic Wi satisfy the relationships.(4)wij=dimnullA−λiInj−dimnullA−λiInj−1=rankA−λiInj−1−rankA−λiInj,j=1,2,…,hi.Note that,(5)h=maxihi=maxjs1j,is the index of nilpotency of A, i.e., the smallest positive integer for which.(6)dimnullA−λInhn.Letk1≥k2≥⋯≥0 be a partition of the positive integer k (i.e., k=k1+k2+⋯) and let mj denotes the number of the integers ki, which are greater or equal to j. Then, the integers mj form a partition of k, known as a conjugate partition of k. The relationship between the Segre and Weyr characteristics of A associated with an eigenvalue λi with multiplicity νi, is revealed by the fact that these characteristics are conjugated partitions of νi [Horn and Johnson [23]; p. 172]. Thus, the Weyr characteristic determines uniquely the Segre characteristic of A associated with λi and vice versa. Note that, the largest element s1i (the size of the largest Jordan block) of the Segre characteristic associated with λi is equal to the number hi of the elements of the corresponding Weyr characteristic. In turn, the first element wi1 of the Weyr characteristic (the geometric multiplicity of the eigenvalue λi) is equal to the number qi of the elements of the Segre characteristic (the number of Jordan blocks) associated with λi.


## 2.2. Reduction to Jordan–Schur Form

Using unitary similarity transformations with a matrixV, the matrix A∈Cn×n is reduced to the upper block triangular form (the Jordan–Schur form).(7)JS=VHAV=SW1,λ1T˜12…T˜1pSW2,λ2…T˜2p⋱⋮SWp,λp,where the diagonal blocks(8)SWi,λi=λiIwi1N12i…N1,hi−1iN1,hiiλiIwi2…N2,hi−1iN2,hii⋱⋮⋮λiIwi,hi−1Nhi−1,hiiλiIwi,hi∈ℂνi×νi,i=1,2,…,p,are in the so-called staircase form and the off-diagonal blocks T˜ij∈ℂwi×wj have no special form. Thus, the structure of the Jordan–Schur form is determined by the Weyr characteristics associated with λi,i=1,2,…,p. These characteristics can be used to find easily the corresponding Segre characteristics of A.A possible way to determine the Jordan–Schur form of a matrix is to implement a modification of the algorithm of Kågström and Ruhe [Kågström and Ruhe [21, 22]] intended to find the Jordan form of a matrix. This algorithm is based on the idea that using a sequence of unitary transformations it is possible to find the dimensions of the null spaces nullA−λIj (and hence the elements of the Weyr characteristic) without computing the matrix powers A−λIj, and consequently, the size of the Jordan blocks associated with multiple eigenvalues. As a result, the matrix is reduced to a staircase form and the corresponding algorithm is called a staircase algorithm because of the shape of the resulting matrix (for a detailed history of this algorithm, see [Edelman and Ma [25]]).The reduction of ann-th order complex matrix A to Jordan–Schur form can be carried out by a minor modification of the algorithm of Kågström and Ruhe [Petkov [26]].Concerning the computational cost, the reduction to the Jordan–Schur form requires more operations than the reduction to the Schur form. This is mainly due to the reduction of the diagonal blocks to the staircase form, which needs a large volume of operations if these blocks have a high index of nilpotence, i.e., if their Jordan forms involve high-order blocks. According to the estimates given in [Kågström and Ruhe [22]], in addition to the reduction into Schur form, the sorting of the eigenvalues, the clustering of multiple eigenvalues, and the reduction of the diagonal blocks of the matrix A to staircase form require 6n3+21/3n3hflops where one flop consists of one addition and one multiplication of floating point numbers. (note that, one flop in complex arithmetic is equivalent to four ordinary flops and the corresponding number of operations should be multiplied by four.) For h=n (nonderogatory but defective matrix with one Jordan block), this volume of work is of order 7/3n4+6n3flops and for h=1 (derogatory but nondefective matrix with scalar Jordan blocks) it is approximately 81/3n3flops. This way, for matrices with large Jordan blocks, the reduction to Jordan–Schur form requires a volume of work of order n4. Some increase of the efficiency of the reduction to staircase form is proposed, for instance, in [Guglielmi et al. [27]].Example 1.
Consider the matrix(9)A=1−57−1090−1−12−2−36−1012−110−23−652−57−108.
The following matrix has the Jordan form:(10)J=diag−11−1,−21−2−2.
The eigenvalueλ1=−1 has a Weyr characteristic equal to 1,1 and Segre characteristic equal to (2). The eigenvalue λ2=−2 has a Weyr characteristic equal to 2,1 and the same Segre characteristic.
Using a unitary similarity transformation with a matrix (all results are rounded to eight decimal digits),(11)V=0.500000000.56694671−0.18273332−0.628633410.00000000−0.166666670.566946710.754888190.15931703−0.23570226−0.666666670.000000000.11922195−0.564906590.471404520.16666667−0.566946710.42117723−0.50117978−0.471404520.50000000−0.188982240.452932910.095590220.70710678.
The matrixA is reduced by the modification of the algorithm of Kågström and Ruhe to the Jordan–Schur form.(12)JS=VHAV=−14.53557368−6.29342726−0.4431080930.05203820−1−0.84090739−0.95805781−0.13363062−203.06799138−2−1.85595575−2,
The reduction is carried out with a relative backward error.(13)A−VJSVHFAF=6.9⋅10−16.


## 3. Description of the Algorithms Proposed

The two algorithms proposed in this paper are modifications of the Schur–Parlett and Schur–Fréchet methods for computing the matrix exponential which instead of the Schur form uses the Jordan–Schur form. The implementation of the Jordan–Schur form(14)JS=VHAV=SW1,λ1T12…T1pSW2,λ2…T2p⋱⋮SWp,λp∈ℂn×n,in which the diagonal blocks SWi,λi have distinct eigenvalues λi,i=1,2,…,p, simplifies the corresponding computations and increases the accuracy. More specifically, the benefit of using the Jordan–Schur form is twofold.Firstly, the diagonal blocks have the following structure:(15)SW,λ=λIw1N12…N1,h−1N1,hλIw2…N2,h−1N2,h⋱⋮⋮λIwh−1Nh−1,hλIwh∈ℂν×ν,Nij∈ℂwi×wj,where to simplify the notation, the Weyr characteristic associated with λ is denoted by(16)W=WS,λ=w1,w2,…,wh.The corresponding Segre characteristic by(17)S=SS,λ=s1,s2,…,sq,ν=∑k=1hwk=∑ℓ=1qsℓ,and h=s1,q=w1. Due to the structure of these blocks, their exponentials can be found by finite Taylor series accessing only the nonzero blocks Nkℓ.Secondly, the Jordan–Schur form is in upper block triangular form with the numerical multiple eigenvalues already clustered into blocks. As shown in the following theory, this facilitates significantly the usage of the Schur–Parlett and Schur–Fréchet method. Moreover, due to the special structure of the matricesTii and Tjj, the solution of the Sylvester matrix equations arising in determining the off-diagonal blocks of eJS by both methods, can be carried out more efficiently in comparison with the original methods.


### 3.1. Determine the Exponentials of the Diagonal Blocks

Denote byN=SW,λ−λIν∈ℂν×ν the nilpotent part of a diagonal block, which has an index of nilpotence equal to 1≤h=s1≤ν. To determine the exponentials of the diagonal blocks, one exploits the following relationship:(18)expSW,λ=expλexpN.Thus, the problem of computingexpSW,λ reduces to the evaluation of the exponential of a nilpotent matrix N with a staircase structure. The Taylor series expansion,(19)expN=Iν+N1!+N22!+⋯+Nh−1h−1!,used to determine the exponential contains only the first h−1 powers of N (if one doesn’t count N0=Iν). Note that, h=ν only in the extreme case of one Jordan block associated with λ of order ν when WS,λ=1,1,…,1 and SS,λ=ν (nonderogatory but defective matrix). In the other extreme case of SS,λ=1,1,…,1,WS,λ=ν (derogatory but nondefective matrix), the index of nilpotence h is equal to 1 and the exponential is found directly from(20)expSW,λ=expλIν.Forh=2, one has that expN=Iν+N so that the raising of N to power is not necessary. In the generic case, it is fulfilled that 1<h<ν which allows obtaining the exponential with a few terms of the Taylor series expansion. Apart from increasing efficiency, this helps avoid some numerical difficulties related to the decision about the termination of the expansion. Since the nilpotent matrix is in staircase form, it is possible to achieve further reduction of the number of computational operations in evaluating the powers of this matrix. Denoting by Fm the m th power of N, one has that(21)Fm=NFm−1=0N12N13…N1h0N23…N2h0⋱⋮⋱Nh−1,h00Fm−11m…Fm−11h0…Fm−12h⋱⋮0,m=2,3,…,h−1,F1=N.Hence, the powerFm,m=2,3,…,h−1 of the nilpotent matrix N with staircase structure can be evaluated from Fm−1 for h>2 using only the nonzero blocks Nij of N and Yij=Fm−1ij of Fm−1. In computing the powers of the nilpotent matrix, one can use the scaling and squaring technique which makes the norm of N less than one. The scaling is necessary in case of the large norm of N and small values of the index of nilpotence h when the norm of powers Nm decreases quickly with m which is related to cancelations in the floating point arithmetic and the introduction of large rounding errors. The scaling and squaring exploit the relationship eN=eN/qq where q=2s and s is chosen so that ‖N/q‖∞<1. Hence,(22)eN=eN/2s2s.Note that, the squaring of the matrixeN can be carried out taking into account that this matrix is upper triangular.Under the assumptions for a matrix with equal eigenvalue multiplicitiesν=n/p and equal elements of the Weyr characteristic wj=v/h, one obtains the following overestimate of the number of operations necessary to compute the exponentials of pν×ν diagonal blocks,(23)ph4+31/3h3−29h2+482/3h−2424n3p3h3flops.The maximum volume of operations is obtained forp=1. Hence, for h=ν=n,wj=1 (this is the case of a single Jordan block of size n), the estimated volume of operations is about 1/24n4+5/36n3flops and for h=3wj=n/3 it is equal to 4/3n3flops. Note that, the computation of the first n−1 powers of a full n×n triangular matrix requires n4/6+n3/6flops.


### 3.2. Determining the Off-Diagonal Blocks by the Schur–Parlett Method

In [Parlett [28]], Parlett developed a recurrent method for finding the off-diagonal elements of a matrix function F=fT using the upper triangular (Schur) form T=QHAQ of a given matrix A, where Q is unitary. The method is based on the commutativity relation FT=TF and requires knowing the diagonal elements fii=ftii of F. Since T is upper triangular, so is F which allows to calculate a superdiagonal of F at a time starting with fii. A disadvantage of the Parlett recurrence is that it breaks down in the presence of equal diagonal elements (eigenvalues) tii and tjj of T and can produce large errors in case of close eigenvalues. That is why the idea is to use the block-diagonal form of T=Tij where the diagonal blocks Tii have distinct eigenvalues. In such case, the off-diagonal blocks of F can be computed by using the following block-recurrence.The Schur–Parlett method [Davis and Higham [12]] is as follows:(24)TiiFij−FijTjj=FiiTij−TijFjj+∑d=i+1j−1FidTdj−TidFdj,j=2,3,…,p,i=j−1,j−2,…,1.The recurrence (24) represents a Sylvester matrix equation with respect to Fij. To have a well-conditioned Sylvester equation (24), the eigenvalues of the blocks Tii and Tjj should be well separated for each i,j. Similarly to the scalar case, it is necessary to know in advance the diagonal blocks Fii=fTii. The computation of the exponential of the diagonal blocks can be carried out efficiently by using Padé approximations combined with scaling and squaring [Higham [5]], [Al-Mohy and Higham [6]].The use of the Jordan–Schur form ofA facilitates significantly the usage of Parlett block recurrence (24). The exponential of the Jordan–Schur form JS is determined as(25)expJS=E11E12…E1pE22…E2p⋱⋮Epp,where the diagonal blocks Eii=expSWi,λi,i=1,…,p are already evaluated.The computation of the off-diagonal blocksEij of the matrix exponential implementing the Parlett block recurrence involves the solution of the Sylvester equations.(26)SiiEij−EijSjj=EiiTij−TijEjj+Z,j=2,3,…,p,i=j−1,j−2,…,1,where the blocks Sii and Sjj are of the form (15), Tij are the superdiagonal blocks of the Jordan–Schur form (14) and(27)Z=∑d=i+1j−1EidTdj−TidEdj.Due to the specific structure of the blocksSii and Sjj, equation (26) will be called the staircase Sylvester equation. It can be solved efficiently taking into account the structure of the diagonal blocks Sii and Sjj.Consider the staircase Sylvester equation(28)SAX+XSB=C,where for simplification we use the following notation.(29)SA=λIwa1N12…N1,haλIwa2…N2,ha⋱⋮λIwaha,SB=μIwb1R12…R1,hbμIwb2…R2,hb⋱⋮μIwbhb,and SA∈ℂκ×κ,SB∈ℂν×ν. It is assumed that SA and SB have eigenvalues λ and μ, respectively, with the associated Weyr characteristics.(30)WSA,λ=wa1,wa2,…,waha,WSB,μ=wb1,wb2,…,wbhb.The staircase Sylvester equation (28) can be solved by the Bartels–Stewart algorithm [Bartels and Stewart [29]]. However, in the given case, it is advantageous to exploit the staircase structure of the matrices SA and SB using Algorithm 1 described in [Petkov [26]]. Briefly, this algorithm works as follows. Partitioning the right-hand side C∈ℂκ×ν and the unknown matrix X∈ℂκ×ν as(31)C=Ckℓ,X=Xkℓ,where Ckℓ,Xkℓ∈ℂwak×wbℓ, the blocks of the matrix X are obtained as follows:(32)Xkℓ=1λ−μCkℓ−∑f=k+1haNkfXfℓ+∑g=1ℓ−1XkgRgℓ,ℓ=1,2,…,hb,k=ha,ha−1,…,1.If the elementswak and wbℓ of the Weyr characteristics are large, then this algorithm finds at once a large block Xkℓ instead of a block of size 1×1. Hence, in the cases of multiple eigenvalues and large elements of the Weyr characteristics, the staircase algorithm is very efficient. In fact, the solution of the staircase Sylvester equation (28) in which the matrix SA is κ×κ, and the matrix SB is ν×ν by the staircase algorithm requires(33)12ha−1κ2νha+hb−1κν2hbflops.For large indices of nilpotency,ha=κ,hb=ν all elements of the Weyr characteristics are equal to 1 and the staircase algorithm reduces to the Bartels–Stewart algorithm which requires(34)κκ−1ν+κνν−12flops.By decreasing the index of nilpotency and increasing the elements of the Weyr characteristics, the computational work associated with the staircase algorithm decreases.Finally, having computedexpJS, the exponential of A is determined as(35)expA=VexpJSVH,where V is the matrix of the unitary transformation of A to JS.The volume of the computational work associated with the determination of the off-diagonal blocks ofexpJS is estimated as follows. Assuming again equal eigenvalue multiplicities κi=n/p,νi=n/p and equal Weyr elements of the characteristics w=wai=wbj=v/h,h=ha=hb, and using the estimate (33), the necessary computational work for solving all pp+1/2 Sylvester equations of the type (26) is estimated as(36)p3−pν33h3flops.The maximum of this volume of work is obtained forp=n,ν=1,h=1,w=1 and for large n it is approximately equal to 1/3n3.Example 2.
Consider again the matrix(37)A=1−57−1090−1−12−2−36−1012−110−23−652−57−108.whose Jordan–Schur form was given in Example 1. This form has a 2×2 diagonal block(38)S11=−14.53557368−1,
associated withλ1=−1, a 3×3 diagonal block,(39)S22=−203.06799138−2−1.85595575−2,
associated withλ2=−2, and a 2×3 superdiagonal block,(40)S12=−6.29342726−0.4431080930.05203820−0.84090739−0.95805781−0.13363062.
The exponential of the first diagonal block is computed by the algorithm presented in Section3.1 with scaling equal to 23. The result is(41)expS11=0.367879441.6685443100.36787944.
The exponential of the second diagonal block, evaluated with scaling22, is(42)expS22=0.1353352800.415207480.13533528−0.251176300.13533528.
The solution of the staircase Sylvester (24), carried out by the algorithm described above, is(43)E12=−1.97966816−0.691120264.97076879−0.19554810−0.22279075−0.10901551.
The computed matrix exponentialE^ is(44)E^=1.23897361−0.968303041.23897361−1.936606081.10363832−0.135335280.36787944−0.232544160.46508832−0.32975303−1.238973611.20084720−1.336182482.40169440−1.433391360.13533528−0.367879440.50321472−0.600423600.735758880.87109417−0.968303041.23897361−1.936606081.47151776.
The exponential is computed with a relative error.(45)E^−expAFexpAF=9.12×10−15.


### 3.3. Determining the Off-Diagonal Blocks Using the Schur–Fréchet Method

Consider a block-triangular matrix of the form(46)T=T11T12T22,and assume that the blocks T11∈ℂn1×n1 and T22∈Csn2×n2 are scaled so that their norms are less than 1. The exponential of T has the following form:(47)E=E11E12E22.Suppose that one has already calculatedE11=eT11 and E22=eT22. Then, the block E12 can be computed by the following algorithm, implementing an 8th-order Padé approximation to the function xcothx. (See [Higham [5]] for a detailed exposition of the corresponding method.)The Schur–Fréchet method [Kenney and Laub [15]].(1)
Letα1,α2,…,α8 be the zeros of the polynomial px=x8+990x6+135135x4+4729725x2+34459425(2)
Letβ1,β2,…,β8 be the zeros of the polynomial qx=45x8+13860x6+945945x4+16216200x2+34459425(3)
SetD12=T12,M11=T11,M22=T22(4)
Evaluate the “Sylvester cascade”fori=1,2,…,8Q=In1/2+M11/αiD12+D12In2/2−M22/αiSolve forD12inIn1/2+M11/βiD12+D12In2/2−M22/αi=Qend(5)
D12=E11D12+D12E22/2(6)
E12=D12In exact arithmetic, this algorithm produces an approximationE˜12 satisfying E˜12−E12F<10−16 [Kenney and Laub [15]]. Note that, the zeros αi and βi are purely imaginary. The algorithm requires the solution of eight Sylvester equations.Thus, the Schur–Fréchet method determines recursively an approximation of the off-diagonal blockE12 of the exponential. The recursion involves the solution of the Sylvester equations.(48)In1/2+M11/βiD12+D12In2/2−M22/αi=Q,i=1,2,…,8,with matrices M11=T11, M22=T22 and the initial approximation D12=T12. Note that, T11 is a single diagonal block of order νi (the multiplicity of the eigenvalue λi) but the matrix T22 may have on its diagonal several blocks corresponding to the distinct eigenvalues μj of multiplicities vj.The key observation here is that the matricesM11, M22, In1/2+M11/βi, and In2/2+M22/βi, respectively, are in Jordan–Schur form. This allows for solving equation (48) efficiently taking into account the structure of these matrices. More specifically, in simplified notation, one has to solve the Sylvester equation.(49)FX+XSWg1,μ1G12…G1pSWg2,μ2…G2p⋱⋮SWgp,μp=H,where F=λInf+M∈ℂnf×nf and M is a staircase matrix with Weyr characteristic Wf associated with λ, the diagonal block SWgj,μj=μjIνj+Nj∈ℂνj×νj corresponds to a distinct numerically multiple eigenvalue μj and Nj is a staircase matrix with Weyr characteristic Wgj associated with μj,j=1,…,p.Representing the matricesH and X as(50)H=H1,H2,…,Hp,Hj∈ℂnf×νj,X=X1,X2,…,Xp,Xj∈ℂnf×νj.One obtains the system of staircase Sylvester equations.(51)FXj+XjSWgj,μj=Hj−∑d=1j−1XdGdj,j=1,2,…,p.It is important that (51) is a staircase Sylvester equation of the same type as (28).Having an algorithm to solve (48), it is possible to find an off-diagonal block E12 of the matrix exponential of T utilizing the Schur–Fréchet method.Consider now, ann×n matrix(52)JS=SW1,λ1T12…T1pSW2,λ2…T2p⋱⋮SWp,λp∈ℂn×n,which is in Jordan–Schur form and assumes that the exponentials eSWi,λi,i=1,2,…,p of the diagonal blocks are already determined. Using the Schur–Fréchet method, the off-diagonal blocks of(53)eJS=E11…E1,p−2E1,p−1E1,p⋱⋮⋮⋮Ep−2,p−2Ep−2,p−1Ep−2,pEp−1,p−1Ep−1,pEpp,can be evaluated recursively in the following way. First one can find the block Ep−1,p using the blocks Tp−1,p−1 and Tpp. Then, using the diagonal blocks Tp−2,p−2 and(54)Tp−1,p−1Tp−1,pTpp.It is possible to evaluate the off-diagonal blocksEp−2,p−1,Ep−2,p and so on, finding ultimately the off-diagonal blocks in the first row of eJS.The volume of computational work necessary to determine the off-diagonal blocks by the Schur–Fréchet method is estimated as follows. Assuming as in the case of the Schur–Parlett method, equal eigenvalue multiplicitiesκi=n/p,νi=n/p, equal elements of the Weyr characteristics w=wfi=wgj=v/h,h=hf=hg, and using the estimate (33), the necessary computational work for solving eight times the Sylvester equation of the type (48) by Schur–Fréchet algorithm is estimated as(55)4p3−3p2ν33h3flops.The maximum of this volume of work is obtained forp=n,ν=1,w=1 and for large n it is approximately equal to 4/3n3. Thus, the Schur–Fréchet method requires four times more work than the Schur–Parlett method provided the Jordan–Schur form is known. Note that, the reduction to this form may require a volume of work proportional to n4, so that the difference between the number of operations for the two methods may be negligible in the total cost.Example 3.
Consider again the matrixA used in Examples 1 and 2. Using the Schur–Fréchet algorithm implementing the Jordan–Schur form, the matrix exponential eA is computed with a relative error.(56)E^−expAFexpAF=6.80×10−15.


## 3.1. Determine the Exponentials of the Diagonal Blocks

Denote byN=SW,λ−λIν∈ℂν×ν the nilpotent part of a diagonal block, which has an index of nilpotence equal to 1≤h=s1≤ν. To determine the exponentials of the diagonal blocks, one exploits the following relationship:(18)expSW,λ=expλexpN.Thus, the problem of computingexpSW,λ reduces to the evaluation of the exponential of a nilpotent matrix N with a staircase structure. The Taylor series expansion,(19)expN=Iν+N1!+N22!+⋯+Nh−1h−1!,used to determine the exponential contains only the first h−1 powers of N (if one doesn’t count N0=Iν). Note that, h=ν only in the extreme case of one Jordan block associated with λ of order ν when WS,λ=1,1,…,1 and SS,λ=ν (nonderogatory but defective matrix). In the other extreme case of SS,λ=1,1,…,1,WS,λ=ν (derogatory but nondefective matrix), the index of nilpotence h is equal to 1 and the exponential is found directly from(20)expSW,λ=expλIν.Forh=2, one has that expN=Iν+N so that the raising of N to power is not necessary. In the generic case, it is fulfilled that 1<h<ν which allows obtaining the exponential with a few terms of the Taylor series expansion. Apart from increasing efficiency, this helps avoid some numerical difficulties related to the decision about the termination of the expansion. Since the nilpotent matrix is in staircase form, it is possible to achieve further reduction of the number of computational operations in evaluating the powers of this matrix. Denoting by Fm the m th power of N, one has that(21)Fm=NFm−1=0N12N13…N1h0N23…N2h0⋱⋮⋱Nh−1,h00Fm−11m…Fm−11h0…Fm−12h⋱⋮0,m=2,3,…,h−1,F1=N.Hence, the powerFm,m=2,3,…,h−1 of the nilpotent matrix N with staircase structure can be evaluated from Fm−1 for h>2 using only the nonzero blocks Nij of N and Yij=Fm−1ij of Fm−1. In computing the powers of the nilpotent matrix, one can use the scaling and squaring technique which makes the norm of N less than one. The scaling is necessary in case of the large norm of N and small values of the index of nilpotence h when the norm of powers Nm decreases quickly with m which is related to cancelations in the floating point arithmetic and the introduction of large rounding errors. The scaling and squaring exploit the relationship eN=eN/qq where q=2s and s is chosen so that ‖N/q‖∞<1. Hence,(22)eN=eN/2s2s.Note that, the squaring of the matrixeN can be carried out taking into account that this matrix is upper triangular.Under the assumptions for a matrix with equal eigenvalue multiplicitiesν=n/p and equal elements of the Weyr characteristic wj=v/h, one obtains the following overestimate of the number of operations necessary to compute the exponentials of pν×ν diagonal blocks,(23)ph4+31/3h3−29h2+482/3h−2424n3p3h3flops.The maximum volume of operations is obtained forp=1. Hence, for h=ν=n,wj=1 (this is the case of a single Jordan block of size n), the estimated volume of operations is about 1/24n4+5/36n3flops and for h=3wj=n/3 it is equal to 4/3n3flops. Note that, the computation of the first n−1 powers of a full n×n triangular matrix requires n4/6+n3/6flops.


## 3.2. Determining the Off-Diagonal Blocks by the Schur–Parlett Method

In [Parlett [28]], Parlett developed a recurrent method for finding the off-diagonal elements of a matrix function F=fT using the upper triangular (Schur) form T=QHAQ of a given matrix A, where Q is unitary. The method is based on the commutativity relation FT=TF and requires knowing the diagonal elements fii=ftii of F. Since T is upper triangular, so is F which allows to calculate a superdiagonal of F at a time starting with fii. A disadvantage of the Parlett recurrence is that it breaks down in the presence of equal diagonal elements (eigenvalues) tii and tjj of T and can produce large errors in case of close eigenvalues. That is why the idea is to use the block-diagonal form of T=Tij where the diagonal blocks Tii have distinct eigenvalues. In such case, the off-diagonal blocks of F can be computed by using the following block-recurrence.The Schur–Parlett method [Davis and Higham [12]] is as follows:(24)TiiFij−FijTjj=FiiTij−TijFjj+∑d=i+1j−1FidTdj−TidFdj,j=2,3,…,p,i=j−1,j−2,…,1.The recurrence (24) represents a Sylvester matrix equation with respect to Fij. To have a well-conditioned Sylvester equation (24), the eigenvalues of the blocks Tii and Tjj should be well separated for each i,j. Similarly to the scalar case, it is necessary to know in advance the diagonal blocks Fii=fTii. The computation of the exponential of the diagonal blocks can be carried out efficiently by using Padé approximations combined with scaling and squaring [Higham [5]], [Al-Mohy and Higham [6]].The use of the Jordan–Schur form ofA facilitates significantly the usage of Parlett block recurrence (24). The exponential of the Jordan–Schur form JS is determined as(25)expJS=E11E12…E1pE22…E2p⋱⋮Epp,where the diagonal blocks Eii=expSWi,λi,i=1,…,p are already evaluated.The computation of the off-diagonal blocksEij of the matrix exponential implementing the Parlett block recurrence involves the solution of the Sylvester equations.(26)SiiEij−EijSjj=EiiTij−TijEjj+Z,j=2,3,…,p,i=j−1,j−2,…,1,where the blocks Sii and Sjj are of the form (15), Tij are the superdiagonal blocks of the Jordan–Schur form (14) and(27)Z=∑d=i+1j−1EidTdj−TidEdj.Due to the specific structure of the blocksSii and Sjj, equation (26) will be called the staircase Sylvester equation. It can be solved efficiently taking into account the structure of the diagonal blocks Sii and Sjj.Consider the staircase Sylvester equation(28)SAX+XSB=C,where for simplification we use the following notation.(29)SA=λIwa1N12…N1,haλIwa2…N2,ha⋱⋮λIwaha,SB=μIwb1R12…R1,hbμIwb2…R2,hb⋱⋮μIwbhb,and SA∈ℂκ×κ,SB∈ℂν×ν. It is assumed that SA and SB have eigenvalues λ and μ, respectively, with the associated Weyr characteristics.(30)WSA,λ=wa1,wa2,…,waha,WSB,μ=wb1,wb2,…,wbhb.The staircase Sylvester equation (28) can be solved by the Bartels–Stewart algorithm [Bartels and Stewart [29]]. However, in the given case, it is advantageous to exploit the staircase structure of the matrices SA and SB using Algorithm 1 described in [Petkov [26]]. Briefly, this algorithm works as follows. Partitioning the right-hand side C∈ℂκ×ν and the unknown matrix X∈ℂκ×ν as(31)C=Ckℓ,X=Xkℓ,where Ckℓ,Xkℓ∈ℂwak×wbℓ, the blocks of the matrix X are obtained as follows:(32)Xkℓ=1λ−μCkℓ−∑f=k+1haNkfXfℓ+∑g=1ℓ−1XkgRgℓ,ℓ=1,2,…,hb,k=ha,ha−1,…,1.If the elementswak and wbℓ of the Weyr characteristics are large, then this algorithm finds at once a large block Xkℓ instead of a block of size 1×1. Hence, in the cases of multiple eigenvalues and large elements of the Weyr characteristics, the staircase algorithm is very efficient. In fact, the solution of the staircase Sylvester equation (28) in which the matrix SA is κ×κ, and the matrix SB is ν×ν by the staircase algorithm requires(33)12ha−1κ2νha+hb−1κν2hbflops.For large indices of nilpotency,ha=κ,hb=ν all elements of the Weyr characteristics are equal to 1 and the staircase algorithm reduces to the Bartels–Stewart algorithm which requires(34)κκ−1ν+κνν−12flops.By decreasing the index of nilpotency and increasing the elements of the Weyr characteristics, the computational work associated with the staircase algorithm decreases.Finally, having computedexpJS, the exponential of A is determined as(35)expA=VexpJSVH,where V is the matrix of the unitary transformation of A to JS.The volume of the computational work associated with the determination of the off-diagonal blocks ofexpJS is estimated as follows. Assuming again equal eigenvalue multiplicities κi=n/p,νi=n/p and equal Weyr elements of the characteristics w=wai=wbj=v/h,h=ha=hb, and using the estimate (33), the necessary computational work for solving all pp+1/2 Sylvester equations of the type (26) is estimated as(36)p3−pν33h3flops.The maximum of this volume of work is obtained forp=n,ν=1,h=1,w=1 and for large n it is approximately equal to 1/3n3.Example 2.
Consider again the matrix(37)A=1−57−1090−1−12−2−36−1012−110−23−652−57−108.whose Jordan–Schur form was given in Example 1. This form has a 2×2 diagonal block(38)S11=−14.53557368−1,
associated withλ1=−1, a 3×3 diagonal block,(39)S22=−203.06799138−2−1.85595575−2,
associated withλ2=−2, and a 2×3 superdiagonal block,(40)S12=−6.29342726−0.4431080930.05203820−0.84090739−0.95805781−0.13363062.
The exponential of the first diagonal block is computed by the algorithm presented in Section3.1 with scaling equal to 23. The result is(41)expS11=0.367879441.6685443100.36787944.
The exponential of the second diagonal block, evaluated with scaling22, is(42)expS22=0.1353352800.415207480.13533528−0.251176300.13533528.
The solution of the staircase Sylvester (24), carried out by the algorithm described above, is(43)E12=−1.97966816−0.691120264.97076879−0.19554810−0.22279075−0.10901551.
The computed matrix exponentialE^ is(44)E^=1.23897361−0.968303041.23897361−1.936606081.10363832−0.135335280.36787944−0.232544160.46508832−0.32975303−1.238973611.20084720−1.336182482.40169440−1.433391360.13533528−0.367879440.50321472−0.600423600.735758880.87109417−0.968303041.23897361−1.936606081.47151776.
The exponential is computed with a relative error.(45)E^−expAFexpAF=9.12×10−15.


## 3.3. Determining the Off-Diagonal Blocks Using the Schur–Fréchet Method

Consider a block-triangular matrix of the form(46)T=T11T12T22,and assume that the blocks T11∈ℂn1×n1 and T22∈Csn2×n2 are scaled so that their norms are less than 1. The exponential of T has the following form:(47)E=E11E12E22.Suppose that one has already calculatedE11=eT11 and E22=eT22. Then, the block E12 can be computed by the following algorithm, implementing an 8th-order Padé approximation to the function xcothx. (See [Higham [5]] for a detailed exposition of the corresponding method.)The Schur–Fréchet method [Kenney and Laub [15]].(1)
Letα1,α2,…,α8 be the zeros of the polynomial px=x8+990x6+135135x4+4729725x2+34459425(2)
Letβ1,β2,…,β8 be the zeros of the polynomial qx=45x8+13860x6+945945x4+16216200x2+34459425(3)
SetD12=T12,M11=T11,M22=T22(4)
Evaluate the “Sylvester cascade”fori=1,2,…,8Q=In1/2+M11/αiD12+D12In2/2−M22/αiSolve forD12inIn1/2+M11/βiD12+D12In2/2−M22/αi=Qend(5)
D12=E11D12+D12E22/2(6)
E12=D12In exact arithmetic, this algorithm produces an approximationE˜12 satisfying E˜12−E12F<10−16 [Kenney and Laub [15]]. Note that, the zeros αi and βi are purely imaginary. The algorithm requires the solution of eight Sylvester equations.Thus, the Schur–Fréchet method determines recursively an approximation of the off-diagonal blockE12 of the exponential. The recursion involves the solution of the Sylvester equations.(48)In1/2+M11/βiD12+D12In2/2−M22/αi=Q,i=1,2,…,8,with matrices M11=T11, M22=T22 and the initial approximation D12=T12. Note that, T11 is a single diagonal block of order νi (the multiplicity of the eigenvalue λi) but the matrix T22 may have on its diagonal several blocks corresponding to the distinct eigenvalues μj of multiplicities vj.The key observation here is that the matricesM11, M22, In1/2+M11/βi, and In2/2+M22/βi, respectively, are in Jordan–Schur form. This allows for solving equation (48) efficiently taking into account the structure of these matrices. More specifically, in simplified notation, one has to solve the Sylvester equation.(49)FX+XSWg1,μ1G12…G1pSWg2,μ2…G2p⋱⋮SWgp,μp=H,where F=λInf+M∈ℂnf×nf and M is a staircase matrix with Weyr characteristic Wf associated with λ, the diagonal block SWgj,μj=μjIνj+Nj∈ℂνj×νj corresponds to a distinct numerically multiple eigenvalue μj and Nj is a staircase matrix with Weyr characteristic Wgj associated with μj,j=1,…,p.Representing the matricesH and X as(50)H=H1,H2,…,Hp,Hj∈ℂnf×νj,X=X1,X2,…,Xp,Xj∈ℂnf×νj.One obtains the system of staircase Sylvester equations.(51)FXj+XjSWgj,μj=Hj−∑d=1j−1XdGdj,j=1,2,…,p.It is important that (51) is a staircase Sylvester equation of the same type as (28).Having an algorithm to solve (48), it is possible to find an off-diagonal block E12 of the matrix exponential of T utilizing the Schur–Fréchet method.Consider now, ann×n matrix(52)JS=SW1,λ1T12…T1pSW2,λ2…T2p⋱⋮SWp,λp∈ℂn×n,which is in Jordan–Schur form and assumes that the exponentials eSWi,λi,i=1,2,…,p of the diagonal blocks are already determined. Using the Schur–Fréchet method, the off-diagonal blocks of(53)eJS=E11…E1,p−2E1,p−1E1,p⋱⋮⋮⋮Ep−2,p−2Ep−2,p−1Ep−2,pEp−1,p−1Ep−1,pEpp,can be evaluated recursively in the following way. First one can find the block Ep−1,p using the blocks Tp−1,p−1 and Tpp. Then, using the diagonal blocks Tp−2,p−2 and(54)Tp−1,p−1Tp−1,pTpp.It is possible to evaluate the off-diagonal blocksEp−2,p−1,Ep−2,p and so on, finding ultimately the off-diagonal blocks in the first row of eJS.The volume of computational work necessary to determine the off-diagonal blocks by the Schur–Fréchet method is estimated as follows. Assuming as in the case of the Schur–Parlett method, equal eigenvalue multiplicitiesκi=n/p,νi=n/p, equal elements of the Weyr characteristics w=wfi=wgj=v/h,h=hf=hg, and using the estimate (33), the necessary computational work for solving eight times the Sylvester equation of the type (48) by Schur–Fréchet algorithm is estimated as(55)4p3−3p2ν33h3flops.The maximum of this volume of work is obtained forp=n,ν=1,w=1 and for large n it is approximately equal to 4/3n3. Thus, the Schur–Fréchet method requires four times more work than the Schur–Parlett method provided the Jordan–Schur form is known. Note that, the reduction to this form may require a volume of work proportional to n4, so that the difference between the number of operations for the two methods may be negligible in the total cost.Example 3.
Consider again the matrixA used in Examples 1 and 2. Using the Schur–Fréchet algorithm implementing the Jordan–Schur form, the matrix exponential eA is computed with a relative error.(56)E^−expAFexpAF=6.80×10−15.


## 4. Numerical Considerations

The error in the evaluation of the matrix exponential depends on the conditioning of the exponential and the round-off properties of the numerical algorithm used.The sensitivity of the matrix exponential to perturbations inA is studied by several authors, see, for instance, [Van Loan [30]], [Kenney and Laub [31]], and [Mathias [32]]. The relative condition number of the matrix exponential that characterizes its sensitivity to small perturbations ΔA of A, is defined as [Higham [11]; Ch. 3].(57)condE≔limε⟶0supΔA≤εAeA+ΔA−eAeA.This condition number is determined by using the Fréchet derivative of the exponential whose computation can be carried out by Algorithm 10.27 presented in [Higham [13]; Ch. 10], or by the algorithms given in [Mathias [33]] and [Al-Mohy and Higham [34]].


### 4.1. The Reduction to Jordan–Schur Form

As mentioned earlier, the reduction to the Jordan–Schur form is numerically stable due to the implementation of unitary transformations. However, it should be emphasized that the numerical stability does not guarantee that the clustering of numerically multiple eigenvalues and the determination of the size of the Jordan blocks are always done in an optimal way. The extensive experiments with this algorithm show that it works reliably except in the case of a very ill-conditioned system of eigenvectors and principal vectors when any clustering method will have serious difficulties in grouping the eigenvalues.


### 4.2. Computing the Off-Diagonal Blocks by the Schur–Parlett Method

In the present version of the Schur–Parlett algorithm, the exponentials of the diagonal blocks are found by using the finite term Taylor series combined with scaling and squaring. As concluded by Higham [Higham [11]; Ch. 10], no clear general conclusions can be drawn about the stability of the scaling and squaring algorithm. Nevertheless, this algorithm is generally considered stable and is implemented, for instance, in the MATLAB™ function expm.The next step of the algorithm is the Parlett block recursion which requires the solution of a Sylvester equation in the form (26) with respect to the unknown off-diagonal block Eij. In the case of close eigenvalues of the blocks Tii and Tjj, this equation can be ill-conditioned. More specifically, let the Sylvester equation be written as(58)SAX−XSB=C,where SA∈ℂκ×κ,SB∈ℂν×ν,C∈ℂκ×ν and let Ω=Iν⊗SA−SBT⊗Iκ. Then, the relative change ΔXF/XF in the solution of this equation due to small relative perturbations ΔSA/SAF,ΔSB/SBF and ΔC/CF, is given by [Higham [35]; Ch. 16].(59)ΔXFXF≤3condsylvη,where(60)condsylv=Ω−1FSAFXT⊗Iκ,−SBFIν⊗X,−CFIκ⋅ν2XF,is the relative condition number of the Sylvester equation and(61)η=maxΔSAFSAF,ΔSBFSBF,ΔCFCF,is the maximum relative perturbation in the data. Now, if the separation,(62)sepSA,SB=minX≠0SAX−XSBFXF,of the matrices SA and SB [Stewart [36]] is small, then the condition number of the Sylvester equation is large and the sensitivity of this equation can be large even in the case of a well-conditioned matrix exponential. (This follows from the fact that Ω−12=1/sepSA,SB). Hence, the Jordan–Schur algorithm implementing the Schur–Parlett method may behave unstably in the presence of blocks in the Jordan–Schur form that has a small separation between them. Note that, this instability is common for all methods implementing the Parlett block recurrence and is not specific to the use of the Jordan–Schur form itself.


### 4.3. Computing the Off-Diagonal Blocks by the Schur–Fréchet Method

The main computational task in the Schur–Fréchet method is the solution of the system of staircase Sylvester equations (the “Sylvester cascade”) (48). Unlike the Sylvester equation (26), solved by the Schur–Parlett method, the sensitivity of (48) does not depend on the closeness of the eigenvalues of Jordan–Schur blocks. Moreover, the separation of matrices F and −SWgj,μj in (51) is always of order 1 which leads to a low condition number of the Sylvester equation. For this reason, the Schur–Fréchet algorithm almost always produces more accurate results than the Schur–Parlett algorithm. The Schur–Fréchet algorithm is especially appropriate in cases when the blocks of the Jordan–Schur form have close eigenvalues.


### 4.4. Computational Complexity

The number of computational operations in complex flops required by the Schur–Parlett and Schur–Fréchet algorithm is summarized in Tables1 and 2. The estimates of the computation cost are obtained under the assumptions ν=n/p,wj=ν/hs.Table 1
Number of computational operations for reduction into Jordan–Schur form and computing the diagonal blocks of the matrix exponential.
OperationNumber of flopsReduction ofA into Schur form T=UHAU≈25n3Reduction ofT into Jordan–Schur form JS6n3+21/3n3hCase ofh=181/3n3Case ofh=n21/3n4+6n3Computing the exponentials of the diagonal blocksh4+31/3h3−29h2+482/3h−24n3/24h3Case ofh=1—Case ofh=n1/24n4+5/36n3−29/24n2+21/36n−1Total: reduction into Jordan–Schur form and computing the exponentials of the diagonal blocks57h4+31/3h3−29h2+482/3h−24n3/24h3+31n3Case ofh=1331/3n3Case ofh=n29/24n4+315/36n3−29/24n2+21/36n−1Table 2
Number of computational operations required by the Schur–Parlett and Schur–Fréchet algorithm.
OperationNumber of flopsTotal number of flopsComputing the exponential of the Jordan–Schur form by the Schur–Parlett methodp3−pν3/3h357h4+31/3h3−29h2+482/3h−24n3/24h3+31n3+p3−pν3/3h3Caseh=1p=1,ν=nn3−n/3332/3n3Caseh=np=n,ν=1—29/24n4+315/36n3Computing the exponential of the Jordan–Schur form by the Schur–Fréchet method4p3−3p2ν3/3h357h4+31/3h3−29h2+482/3h−24n3/24h3+31n3+4p3−3p2ν3/3h3Caseh=1p=1,ν=n4n3−3n2/3342/3n3Caseh=np=n,ν=1—29/24n4+315/36n3In Table1, we give the approximate number of operations necessary for reduction into Jordan–Schur form and computation of the diagonal blocks of the exponential. These computations are common for both methods and the corresponding volume of work is given for the two extreme cases:Derogatory but nondefective matrix with scalar Jordan blocksh=1Nonderogatory but defective matrix with a single Jordan blockh=nAs seen from Table1, the number of operations in these cases varies between 331/3n3 and 29/24n4+315/36n3+On2.In Table2, we give the number of operations for the Schur–Parlett and Schur–Fréchet algorithm including the extreme cases h=1 and h=n. The third column of the table contains the total sum of operations necessary for reduction into the Jordan–Schur form, computing the exponentials of the diagonal blocks and computing the off-diagonal blocks by the corresponding method. Clearly, the work necessary for the reduction into Jordan–Schur form dominates the total operation count which is approximately the same for both methods. In the most difficult case h=n, both methods require 29/24n4+315/36n3 flops.


## 4.1. The Reduction to Jordan–Schur Form

As mentioned earlier, the reduction to the Jordan–Schur form is numerically stable due to the implementation of unitary transformations. However, it should be emphasized that the numerical stability does not guarantee that the clustering of numerically multiple eigenvalues and the determination of the size of the Jordan blocks are always done in an optimal way. The extensive experiments with this algorithm show that it works reliably except in the case of a very ill-conditioned system of eigenvectors and principal vectors when any clustering method will have serious difficulties in grouping the eigenvalues.


## 4.2. Computing the Off-Diagonal Blocks by the Schur–Parlett Method

In the present version of the Schur–Parlett algorithm, the exponentials of the diagonal blocks are found by using the finite term Taylor series combined with scaling and squaring. As concluded by Higham [Higham [11]; Ch. 10], no clear general conclusions can be drawn about the stability of the scaling and squaring algorithm. Nevertheless, this algorithm is generally considered stable and is implemented, for instance, in the MATLAB™ function expm.The next step of the algorithm is the Parlett block recursion which requires the solution of a Sylvester equation in the form (26) with respect to the unknown off-diagonal block Eij. In the case of close eigenvalues of the blocks Tii and Tjj, this equation can be ill-conditioned. More specifically, let the Sylvester equation be written as(58)SAX−XSB=C,where SA∈ℂκ×κ,SB∈ℂν×ν,C∈ℂκ×ν and let Ω=Iν⊗SA−SBT⊗Iκ. Then, the relative change ΔXF/XF in the solution of this equation due to small relative perturbations ΔSA/SAF,ΔSB/SBF and ΔC/CF, is given by [Higham [35]; Ch. 16].(59)ΔXFXF≤3condsylvη,where(60)condsylv=Ω−1FSAFXT⊗Iκ,−SBFIν⊗X,−CFIκ⋅ν2XF,is the relative condition number of the Sylvester equation and(61)η=maxΔSAFSAF,ΔSBFSBF,ΔCFCF,is the maximum relative perturbation in the data. Now, if the separation,(62)sepSA,SB=minX≠0SAX−XSBFXF,of the matrices SA and SB [Stewart [36]] is small, then the condition number of the Sylvester equation is large and the sensitivity of this equation can be large even in the case of a well-conditioned matrix exponential. (This follows from the fact that Ω−12=1/sepSA,SB). Hence, the Jordan–Schur algorithm implementing the Schur–Parlett method may behave unstably in the presence of blocks in the Jordan–Schur form that has a small separation between them. Note that, this instability is common for all methods implementing the Parlett block recurrence and is not specific to the use of the Jordan–Schur form itself.


## 4.3. Computing the Off-Diagonal Blocks by the Schur–Fréchet Method

The main computational task in the Schur–Fréchet method is the solution of the system of staircase Sylvester equations (the “Sylvester cascade”) (48). Unlike the Sylvester equation (26), solved by the Schur–Parlett method, the sensitivity of (48) does not depend on the closeness of the eigenvalues of Jordan–Schur blocks. Moreover, the separation of matrices F and −SWgj,μj in (51) is always of order 1 which leads to a low condition number of the Sylvester equation. For this reason, the Schur–Fréchet algorithm almost always produces more accurate results than the Schur–Parlett algorithm. The Schur–Fréchet algorithm is especially appropriate in cases when the blocks of the Jordan–Schur form have close eigenvalues.


## 4.4. Computational Complexity

The number of computational operations in complex flops required by the Schur–Parlett and Schur–Fréchet algorithm is summarized in Tables1 and 2. The estimates of the computation cost are obtained under the assumptions ν=n/p,wj=ν/hs.Table 1
Number of computational operations for reduction into Jordan–Schur form and computing the diagonal blocks of the matrix exponential.
OperationNumber of flopsReduction ofA into Schur form T=UHAU≈25n3Reduction ofT into Jordan–Schur form JS6n3+21/3n3hCase ofh=181/3n3Case ofh=n21/3n4+6n3Computing the exponentials of the diagonal blocksh4+31/3h3−29h2+482/3h−24n3/24h3Case ofh=1—Case ofh=n1/24n4+5/36n3−29/24n2+21/36n−1Total: reduction into Jordan–Schur form and computing the exponentials of the diagonal blocks57h4+31/3h3−29h2+482/3h−24n3/24h3+31n3Case ofh=1331/3n3Case ofh=n29/24n4+315/36n3−29/24n2+21/36n−1Table 2
Number of computational operations required by the Schur–Parlett and Schur–Fréchet algorithm.
OperationNumber of flopsTotal number of flopsComputing the exponential of the Jordan–Schur form by the Schur–Parlett methodp3−pν3/3h357h4+31/3h3−29h2+482/3h−24n3/24h3+31n3+p3−pν3/3h3Caseh=1p=1,ν=nn3−n/3332/3n3Caseh=np=n,ν=1—29/24n4+315/36n3Computing the exponential of the Jordan–Schur form by the Schur–Fréchet method4p3−3p2ν3/3h357h4+31/3h3−29h2+482/3h−24n3/24h3+31n3+4p3−3p2ν3/3h3Caseh=1p=1,ν=n4n3−3n2/3342/3n3Caseh=np=n,ν=1—29/24n4+315/36n3In Table1, we give the approximate number of operations necessary for reduction into Jordan–Schur form and computation of the diagonal blocks of the exponential. These computations are common for both methods and the corresponding volume of work is given for the two extreme cases:Derogatory but nondefective matrix with scalar Jordan blocksh=1Nonderogatory but defective matrix with a single Jordan blockh=nAs seen from Table1, the number of operations in these cases varies between 331/3n3 and 29/24n4+315/36n3+On2.In Table2, we give the number of operations for the Schur–Parlett and Schur–Fréchet algorithm including the extreme cases h=1 and h=n. The third column of the table contains the total sum of operations necessary for reduction into the Jordan–Schur form, computing the exponentials of the diagonal blocks and computing the off-diagonal blocks by the corresponding method. Clearly, the work necessary for the reduction into Jordan–Schur form dominates the total operation count which is approximately the same for both methods. In the most difficult case h=n, both methods require 29/24n4+315/36n3 flops.


## 5. Numerical Experiments

In assessing the numerical behavior of the algorithms for computing the matrix exponential, the authors commonly use the test matrices included in the sets described in [Higham [37]], [Fasi and Higham [38]], and [Wright [39]]. Many of the matrices from these sets have distinct eigenvalues so that their Jordan–Schur form coincides with the Schur form. In such cases, the Jordan–Schur algorithms for computing the exponential do not have some advantages over the Schur algorithms. That is why in this section we present four experiments involving matrices with different Jordan structures so that the Jordan–Schur form differs significantly from the Schur form.The computations in the paper are carried out with MATLAB™ Version 9.9 (R2020b) [40] using IEEE double precision arithmetic. The tolerances ein and tol, used in the reduction to Jordan–Schur form, are both taken equal to 10−14. The relative condition number of the matrix exponential is determined by the function expm_cond from the matrix function toolbox [Higham [13]]. The computations are done using the Schur–Parlett algorithm, the Schur–Fréchet algorithm, the function expm in MATLAB™ implementing Padé approximations combined with scaling and squaring [Al-Mohy and Higham [6]], and the function funm in MATLAB™ which implements the Schur–Parlett algorithm [Davis and Higham [12]].Experiment 1.
The first experiment demonstrates the behavior of the two versions of the Jordan–Schur algorithm for a family of problems with nonderogatory matrices and with increasing conditioning of the exponential. In the given case, the 8th-order matrixA is constructed as(63)A=QJQ−1,where(64)J=diag−210000−210000−210000−210000−2,410041004,is the desired Jordan form and Q is a nonsingular matrix. The matrix Q is taken as(65)Q=W⋅W−T,Q−1=WT⋅W−1,where the matrix W is chosen as(66)W=1−10010−1τ10−1001010010−11100110101011−11.
The parameterτ varies between 0 and 30 and the condition number of the matrix exponential varies between 1.08⋅103 and 4.74⋅109, respectively. An accurate approximation of the exact exponential of A is determined as E=QexpJQ−1, where expJ is computed using the corresponding analytical expression.
The results of this experiment are shown in Figure1. The Schur–Parlett and the Schur–Fréchet algorithm produce almost identical results, close to the results obtained by the function funm.Figure 1
Relative errors for Experiment1.Experiment 2.
This experiment illustrates the difficulties with the Schur–Parlett version of the algorithm in case of a small separation between the diagonal blocks of the Jordan–Schur form.
The matrixA in the experiment is constructed as in Experiment 1,(67)A=QJQ−1,Q=W⋅W−T,with the matrix W taken as(68)W=1−10011−1111−100111001−1−1110011−1101011−11.
Note that, the matricesQ and Q−1=WT⋅W−1 have integer entries. The condition number with respect to the inversion of the similarity transformation matrix Q is condFQ=333.40. The Jordan form J is chosen with two blocks. The first block is of 5th order and its eigenvalue λ1 varies between 3 and 3.9, while the second block is of 3rd order with an eigenvalue which is constant and equal to 4. As in Example 1, the exact exponential is obtained using the analytical expression of the exponential expJ. With the increasing of λ1, the separation between the two blocks T11 and T22 decreases quickly, which leads to large errors in the solution of Sylvester (26) appearing in the block-recurrence. The relative error in computing the exponential along with the error bound condE⋅eps where eps is the MATLAB™ function eps, eps=2−52, as a function of τ=λ1−λ2 is shown in Figure 2.
The experiment presented shows the potential instability of the Schur–Parlett method and its Jordan–Schur version in the case of blocks with close eigenvalues. The function funm also produces large errors, as a result of implementing the Schur–Parlett method. For such matrices, it is appropriate to use the Schur–Fréchet algorithm, the Padé approximation combined with scaling and squaring, or the alternative methods that are based on the Taylor series [Caliari and Zivcovich [8]].Figure 2
Relative errors for Experiment2.Experiment 3.
This experiment illustrates the properties of the algorithms proposed in the case of matrices with defective eigenvalues and large elements of the Weyr characteristic. The matrixA is of even order and is taken as(69)A=QJQ−1,where the Jordan form.(70)J=diagλ1λ,λ1λ,…,λ1λ.
Consists of2×2 blocks with equal eigenvalues and the nonsingular transformation matrix Q is constructed as proposed in [Bavely and Stewart [41]],(71)Q=H2SH1,Q−1=H1S−1H2,H1=In−2uuTn,H2=In−2vvTn,u=1,1,1,…,1T,v=1,−1,1,…,−1n−1T,S=diag1,σ,σ2,…,σn−1,where H1,H2 are elementary reflections that are orthogonal and symmetric matrices [Stewart [42]]. The condition number of Q with respect to the inversion is controlled by the variable σ and is equal to σn−1. Note that, the matrices Q and Q−1 are obtained with small rounding errors. An accurate approximation of the exact exponential of A is determined from E=QexpJQ−1, where expJ is evaluated using its analytical expression. In the given case, we take λ=0.1 and the parameter σ controlling the conditioning of Q is chosen equal to 1.25. The Segre characteristic of A associated with λ is equal to SA,λ=2,2,2,…,2 and the corresponding Weyr characteristic is WA,λ=n/2,n/2.
The relative errors of the computed exponential forn=4,8,12,…,80 are shown in Figure 3 for the four methods used earlier. The function expm exhibits the worst performance and for n=68 the relative error of the computed result is equal to 4.50⋅10−1. At the same time, the relative errors of the exponential computed by the two versions of the Jordan–Schur algorithm are identical and for n=68 are approximately equal to 6.67⋅10−4. Even for n=80, the relative errors of these algorithms are smaller than 1. The good results for the Jordan–Schur algorithms are due to the simple structure of the Jordan–Schur form.(72)JS=λIn/2N12λIn/2,which facilities very much the computation of the exponential. The function funm also produces accurate results for this example.Figure 3
Relative errors for Experiment3.Experiment 4.
In this experiment, the matrixA is constructed in the same way as in Experiment 3 but in the Jordan form.(73)J=diagλ1λ,…,λ1λ,λ1λ⋱⋱1λ.
It consists ofn/4 blocks of order 2 whose eigenvalues are equal to λ and one n/2×n/2 single Jordan block that has same eigenvalue. In the given case, we take λ=5.0 and the parameter σ controlling the conditioning of the nonsingular transformation matrix Q is chosen equal to 1.2.
The relative error of the computed exponential for differentn is shown in Figure 4. In the given case, the Jordan–Schur algorithms and the function expm produce results that are bounded by the quantity condE⋅eps while the function funm exhibits instability for n>28. Note that, the Segre characteristic of A associated with λ is n/2,2,2,…,2 and the corresponding Weyr characteristic is n/4+1,n/4,1,1,…,1.Figure 4
Relative errors for Experiment4.


## 6. Conclusions

The two algorithms proposed in this paper take advantage of using the Jordan–Schur form for finding the matrix exponential in the case of multiple eigenvalues. The algorithms are more accurate than some known methods for computing the matrix exponential in certain cases at the price of increasing the computational cost of the overall algorithm due to the reduction to the Jordan–Schur form. The numerical experiments confirm that the Jordan–Schur algorithms for computing the matrix exponential are appropriate for matrices with multiple eigenvalues and are especially efficient in cases of large elements of the Weyr characteristics.



---

*Source: 4908185-2023-02-20.xml*