Computing the matrix exponential with the double exponential formula

Fuminori Tatsuoka; Tomohiro Sogabe; Tomoya Kemmochi; Shao-Liang Zhang

doi:10.1515/spma-2024-0013

Artikel Open Access

Computing the matrix exponential with the double exponential formula

Fuminori Tatsuoka , Tomohiro Sogabe , Tomoya Kemmochi und Shao-Liang Zhang

Veröffentlicht/Copyright: 26. Oktober 2024

Veröffentlicht von

Veröffentlichen auch Sie bei De Gruyter Brill

Manuskript einreichen Informationen für Autor*innen Erkunden Sie dieses Fachgebiet

Aus der Zeitschrift Special Matrices Band 12 Heft 1

Abstract

This article considers the computation of the matrix exponential e A with numerical quadrature. Although several quadrature-based algorithms have been proposed, they focus on (near) Hermitian matrices. In order to deal with non-Hermitian matrices, we use another integral representation including an oscillatory term and consider applying the double exponential (DE) formula specialized to Fourier integrals. The DE formula transforms the given integral into another integral whose interval is infinite, and therefore, it is necessary to truncate the infinite interval. In this article, to utilize the DE formula, we analyze the truncation error and propose two algorithms. The first one approximates e A with the fixed mesh size, which is a parameter in the DE formula affecting the accuracy. The second one computes e A based on the first one with automatic selection of the mesh size depending on the given error tolerance.

Keywords: matrix function; matrix exponential; numerical quadrature; double exponential formula

MSC 2010: 65F60; 65D30

1 Introduction

The exponential of a matrix A ∈ C n × n is defined as follows:

e A ≔ I + A + 1 2 ! A 2 + 1 3 ! A 3 + ⋯ ,

and it arises in several situations in scientific computing. One of the applications is the exponential integrator, a class of numerical solvers for stiff ordinary differential equations [13]. In recent years, e A also arises in the analysis of directed graphs [4]. Thus, several classes of computational methods of e A and e A b ( b ∈ C n ) have been proposed. For example, methods based on Padé approximation of e z at z = 0 [1,2,7], methods based on projections onto Krylov-like subspaces [8,18], methods based on the best approximation of e z on ( − ∞ , 0 ] [19], and methods based on numerical quadrature [5,24]. For a more detailed review of the computational methods of the matrix exponential, see, e.g. [14,15].

In general, quadrature-based methods have two advantages: these algorithms can compute e A b without computing e A itself when A is large and either sparse or structured, and it is possible to easily make these algorithms parallel in the sense that the integrand can be computed independently on each abscissa [23, Sect. 18]. For similar reasons, quadrature-based algorithms for other matrix functions have been developed in recent years [21,22]. Thus, we consider quadrature-based algorithms in this article.

Quadrature-based algorithms in [5,24] focus on (nearly) Hermitian matrices. In detail, they are proposed for the efficient computation of Bromwich integrals whose integrand has singularities on or near the negative real axis. For the computation of e A , (all) the eigenvalues of A are the singularities of the integrand. Hence, when A has eigenvalues with large imaginary parts, the singularities are also located far from the negative real axis, and these algorithms would be inaccurate or may not converge.

The motivation of this study is to construct quadrature-based algorithms that can be used for non-Hermitian matrices. Here, we consider another integral representation

(1) e A = 2 π ∫ 0 ∞ x sin ( x ) ( x 2 I + A 2 ) − 1 d x .

The derivation of (1) is presented in A. The integral representation (1) holds when all the eigenvalues of A lie in the open left half plane { z ∈ C : Re ( z ) < 0 } , but this condition can be assumed without loss of generality because e A + s I = e s e A ( s ∈ C ) . Unlike the Cauchy integral used in [5,24], the integral representation (1) allows us to avoid considering its integral path even when A is a non-Hermitian matrix.

Since the integrand in (1) includes an oscillatory term sin ( x ) and the interval is half infinite, it is difficult to compute (1) by using typical quadrature formula such as Gaussian quadrature. To compute (1), we consider applying the double exponential (DE) formula specialized to Fourier integrals [17].

In this article, we propose algorithms using the DE formula to compute e A b (or e A ) for non-Hermitian matrices with numerical quadrature. The DE formula transforms a given integral into another integral suited for the trapezoidal rule. Because the transformed integral interval is infinite, we need to truncate the infinite sum of the discretized integral appropriately into a finite one. Thus, we propose a truncation method of the infinite sum based on error analysis. In addition, we show an automatic quadrature algorithm that selects the mesh size of the trapezoidal rule depending on the given tolerance.

The organization of this article is as follows: The DE formula used in this article is introduced in Section 2. In Section 3, we analyze the truncation error and propose algorithms. Numerical results are presented in Section 4, and we conclude this article in Section 5.

2 The DE formula for Fourier integrals

The DE formula exploits the fact that the trapezoidal rule for the integrals of analytic functions on the real line converges exponentially [20]. Several types of change of variable are proposed to deal with different types of integrals. In [16], integral forms of

(2) ℐ = ∫ 0 ∞ g ( x ) sin ( x ) d x

are considered, where g is a scalar function decaying polynomially as x → ∞ .

The DE formula for (2) is as follows. We first select the mesh size h > 0 to be used for the trapezoidal rule. Next, we apply a change of variable x = x h ( t ) such that x h ′ ( t ) decays double exponentially as t → − ∞ and x h ( t ) − π t ⁄ h decays double exponentially as t → ∞ . Then, we compute the transformed integral with the trapezoidal rule:

ℐ = ∫ − ∞ ∞ x h ′ ( t ) g ( x h ( t ) ) sin ( x h ( t ) ) d t ≈ h ∑ k = − ∞ ∞ x h ′ ( k h ) g ( x h ( k h ) ) sin ( x h ( k h ) ) .

The summand decays double exponentially as ∣ k ∣ → ∞ , and therefore, the infinite sum can be approximated by a sum of not too many terms. Examples of such change of variable are as follows:

x h ( t ) = π h t 1 − exp ( − α sinh t ) , ( α = 6 ) ,

which is proposed in [16], or

(3) x h ( t ) = π h t 1 − exp ( − 2 t − α ( 1 − e − t ) − β ( e t − 1 ) ) β = 1 4 , α = β 1 + log ( 1 + π ⁄ h ) ⁄ 4 h ,

which is proposed in [17]. The implementation notes for (3) is presented in Appendix B.

The DE formula can be applied to (1) because the integrand is a rational function of A in the sense of [12, Definition 1.2]. See [21, Section 1.1] for more details of the discussion.

3 Computing e A with the DE formula

In this section, we propose algorithms for e A based on the DE formula:

(4) e A = ∫ − ∞ ∞ F h ( t , A ) d t ≈ ∑ k = l r h F h ( k h , A ) ,

where

F h ( t , A ) ≔ x h ′ ( t ) sin ( x h ( t ) ) G ( x h ( t ) , A ) , G ( x , A ) ≔ 2 π x ( x 2 I + A 2 ) − 1 .

We may not write the second argument of F h and G when it is obvious from the context.

The error of the DE formula can be divided into the discretization error and the truncation error:

∫ − ∞ ∞ F h ( t , A ) d t − ∑ k = l r h F h ( k h , A ) = ∫ − ∞ ∞ F h ( t ) d t − ∑ k = − ∞ ∞ h F h ( k h ) + ∑ k = − ∞ ∞ h F h ( k h ) − ∑ k = l r h F h ( k h ) .

We employ a posteriori error estimation technique to presume the discretization error and provide an upper bound on the truncation error. In Section 3.1, we analyze the truncation error for the given interval, and in Section 3.2, we present the two algorithms; notably, one algorithm uses the technique to achieve the required accuracy.

3.1 Truncation error

We show an upper bound on the truncation error as follows:

Proposition 1

Suppose that all the eigenvalues of A ∈ C n × n lie in the left half plane. Let u ( t ) = 1 − π t ⁄ h x h ( t ) , where x h ( t ) is a DE transformation for Fourier integrals. Given l , r ∈ Z such that l < r , x h ( l h ) ≤ min { 1 ⁄ 2 ‖ A − 2 ‖ , π } , and x h ′ ( k h ) ≤ 2 π ⁄ h for k ≥ r + 1 , we have

(5) ∑ k = − ∞ ∞ h F h ( k h ) − ∑ k = l r h F h ( k h ) ≤ 2 h π ∑ k = − ∞ l − 1 x h ′ ( k h ) + 2 π ∑ k = r + 1 ∞ k u ( k h ) 1 − u ( k h ) ( ∥ ( A + i x h ( k h ) I ) − 1 ∥ + ∥ ( A − i x h ( k h ) I ) − 1 ∥ ) ,

where ‖ ⋅ ‖ is any subordinate norm.

Proof

From the triangle inequality, it holds that

∑ k = − ∞ ∞ h F h ( k h ) − ∑ k = l r h F h ( k h ) ≤ ∑ k = − ∞ l − 1 h F h ( k h ) + ∑ k = r + 1 ∞ h F h ( k h ) .

In this proof, we show (5) that

(6) ∑ k = − ∞ l − 1 h F h ( k h ) < 2 h π ∑ k = − ∞ l − 1 x h ′ ( k h ) ,

(7) ∑ k = r + 1 ∞ h F h ( k h ) ≤ 2 π ∑ k = r + 1 ∞ k u ( k h ) 1 − u ( k h ) ( ∥ ( A + i x h ( k h ) I ) − 1 ∥ + ∥ ( A − i x h ( k h ) I ) − 1 ∥ ) .

We first show (6). The assumption ( 0 < ) x h ( l h ) ≤ 1 ⁄ 2 ‖ A − 2 ‖ gives that for k ≤ l , ‖ x h ( k h ) 2 A − 2 ‖ ≤ 1 ⁄ 2 ( < 1 ) . Hence, by using the Neumann expansion, we have

∥ x h ( k h ) 2 ( x h ( k h ) 2 I + A 2 ) − 1 ∥ = ∥ x h ( k h ) 2 A − 2 [ − ( − x h ( k h ) 2 A − 2 ) + I ] − 1 ∥ = ∑ m = 0 ∞ [ − x h ( k h ) 2 A − 2 ] m + 1 ≤ ∑ m = 1 ∞ ∥ x h ( k h ) 2 A − 2 ∥ m ≤ 1 .

In addition, because ( 0 < ) x h ( k h ) ≤ π for k < l , we have

∑ k = − ∞ l − 1 h F h ( k h ) = 2 h π ∑ k = − ∞ l − 1 x h ′ ( k h ) sin ( x h ( k h ) ) x h ( k h ) x h ( k h ) 2 ( x h ( k h ) 2 I + A 2 ) − 1 ≤ 2 h π ∑ k = − ∞ l − 1 x h ′ ( k h ) sin ( x h ( k h ) ) x h ( k h ) ≤ 2 h π ∑ k = − ∞ l − 1 x h ′ ( k h ) .

Next, we show (7). It is easily seen that

∑ k = r + 1 ∞ h F h ( k h ) ≤ h π ∑ k = r + 1 ∞ x h ′ ( k h ) ∣ sin ( x h ( k h ) ) ∣ ∥ 2 x h ( k h ) ( x h ( k h ) 2 I + A 2 ) − 1 ∥ = h π ∑ k = r + 1 ∞ x h ′ ( k h ) ∣ sin ( x h ( k h ) ) ∣ [ ∥ ( A + i x h ( k h ) I ) − 1 ∥ + ∥ ( A − i x h ( k h ) I ) − 1 ∥ ] .

Here, because x h ′ ( k h ) ≤ 2 π ⁄ h for k ≥ r + 1 and

∣ sin ( x h ( k h ) ) ∣ = sin k π 1 − u ( k h ) = sin k π + k π u ( k h ) 1 − u ( k h ) ≤ k π u ( k h ) 1 − u ( k h ) ,

we have (7), and (5) is proved.□

We remark that the transformation (3) satisfies x h ′ ( t ) < 2 π ⁄ h for t ≥ t 0 , where t 0 is the solution of − 2 t − α ( 1 − e − t ) − β ( e t − 1 ) = log ( 1 − 1 ⁄ 2 ) . For more details, see Appendix B.

3.2 Algorithms

In this subsection, we propose two algorithms. The first one computes e A by using the DE formula with the given mesh size h and the second one computes e A with automatic mesh size selection.

In the first algorithm, the input matrix A is shifted so that all the eigenvalues of A lie in the left half plain. Subsequently, the algorithm selects the truncation point of the discretized integral so that the truncation error in 2-norm would be less than the given error tolerance. Then the algorithm calculates the exponential of the shifted matrix. Finally, we obtain e A by scaling the exponential of the shifted matrix exponential. The details of the algorithm are given in Algorithm 1.

Algorithm 1 Computing e A with the DE formula with given mesh size
	Input A ∈ C n × n , mesh size h > 0 , tolerance for the truncation error ε > 0 , shift parameter σ < 0
1:	Compute the right-most eigenvalue of A , λ right
2:	A ˜ = A + ( σ − λ right ) I , ε ˜ = ε ⁄ ∣ e λ right − σ ∣ ⊳ Because of the scaling of A , the tolerance also needs to be scaled.
3:	l , r = GetInterval ( ε ˜ , h , A ˜ )
4:	X ˜ = h ∑ k = l r F h ( k h , A ˜ )
	Output X = e λ right − σ X ˜ ≈ e A

5:	function GetInterval ( ε , h , A ) ⊳ Compute an interval whose truncation error is approximately smaller than ε
6:	Find the maximum l ∈ Z satisfying 2 h ( ∑ k = − ∞ l − 1 x ′ ( k h ) ) ⁄ π ≤ ε ⁄ 2
7:	Find the minimum r ∈ Z satisfying 4 π ‖ A − 1 ‖ 2 ∑ k = r + 1 ∞ k u ( k h ) ⁄ ( 1 − u ( k h ) ) ≤ ε ⁄ 2
8:	return l , r
9:	end function

Algorithm 1 requires the shift parameter σ , and this parameter must be chosen appropriately. Indeed, when σ is positive, the integral representation (1) does not make sense. In addition, when σ is large in the negative direction, it makes inaccurate results because the condition number for e A becomes large as ‖ A ‖ is large, see [9, Section 1]. From the numerical examples in Section 4.2, it may be better to choose σ from [ − 5 , 0 ) .

The computation of the resolvent ( x h ( k h ) 2 I + A ˜ 2 ) − 1 is necessary in Step 4 of Algorithm 1, where A ˜ is the shifted matrix computed at Step 2. The condition number for the resolvent can be as large as the square of the condition number of A ˜ , and it may result in low accuracy. In this case, although the computational cost increases, we can evaluate it as follows:

( x h ( k h ) 2 I + A ˜ 2 ) − 1 = i 2 x h ( k h ) [ ( i x h ( k h ) I + A ˜ ) − 1 − ( − i x h ( k h ) I + A ˜ ) − 1 ]

for accuracy.

To obtain an upper bound on the 2-norm of the truncation error using Proposition 1, it is necessary to compute the norm ‖ ( A ˜ + i x h ( k h ) I ) − 1 ‖ 2 for k ≥ r + 1 . As the computation of the norm at each k is challenging, the algorithm assumes that the upper bound on the norm can be approximated by ‖ A ˜ − 1 ‖ 2 . In other words, the upper bound on the truncation error ‖ ∑ k = r + 1 ∞ h F h ( k h ) ‖ 2 is approximated by

4 π ‖ A ˜ − 1 ‖ 2 ∑ k = r + 1 ∞ k u ( k h ) 1 − u ( k h ) .

Hence, for highly non-normal matrices such that ‖ A ˜ − 1 ‖ 2 is much larger than ‖ ( A ˜ + i x h ( k h ) I ) − 1 ‖ 2 , the computational cost of the algorithm may be large because of the overestimation of the truncation error.

The infinite sums in steps 6 and 7 of Algorithm 1 can be approximated without too much computations because the summands x h ′ ( k h ) and k u ( k h ) ⁄ ( 1 − u ( k h ) ) decay double exponentially. For example, Figure 1 illustrates the value of ∣ x h ′ ( k h ) ∣ and ∣ k u ( k h ) ⁄ ( 1 − u ( k h ) ) ∣ for h = 0.05 , and it shows that these sums can be presumed by computing the first 50 summands.

$Figure 1 The absolute value of the summands in the upper bound on the truncation error. The figure on the left illustrates ∣ x ′ h ( k h ) ∣ | x{^{\prime} }_{h}\left(kh)| in (5), and the one on the right illustrates ∣ k u ( k h ) ⁄ ( 1 − u ( k h ) ) ∣ | ku\left(kh)/\left(1-u\left(kh))| . The transformation (3) is used for this calculation.$

Figure 1

The absolute value of the summands in the upper bound on the truncation error. The figure on the left illustrates ∣ x ′ h ( k h ) ∣ in (5), and the one on the right illustrates ∣ k u ( k h ) ⁄ ( 1 − u ( k h ) ) ∣ . The transformation (3) is used for this calculation.

We next propose an automatic quadrature algorithm. Ooura and Mori proposed an automatic quadrature algorithm in [17, Section 5] that is applicable when the convergence rate of the DE formula, the constant ρ > 0 such that ∣ ℐ − h ∑ k = − ∞ ∞ x h ′ ( k h ) g ( x h ( k h ) ) sin ( x h ( k h ) ) ∣ = O ( exp ( − ρ ⁄ h ) ) , is known.

In the computation of e A ˜ , the error ‖ e A ˜ − h ∑ k = − ∞ ∞ F h ( k h , A ˜ ) ‖ can be bounded by γ exp ( − ρ ⁄ h ) with some positive constants γ and ρ because any element of F h ( t ) is an analytic function on a strip region { z ∈ C : ∣ Im ( z ) ∣ < d } for some d > 0 . However, the convergence rate ρ varies with the eigenvalue distribution of A ˜ because the singularities of the integrand in (4) are solutions z of x h ( z ) 2 = − λ ˜ 2 , where λ ˜ is any eigenvalue of A ˜ . Thus, we add a presumption of the convergence rate to the algorithm. Once we presume γ and ρ , we can select an appropriate mesh size h by solving ∣ e λ right − σ ∣ ε = γ exp ( − ρ ⁄ h ) for h with a given tolerance ε .

In the algorithm, we first select three mesh sizes h 1 > h 2 > h 3 > 0 . Let X ˜ i ( i = 1 , 2 , 3 ) be the computational result of the DE formula for e A ˜ with the mesh size h i . If X ˜ 3 is much more accurate than X ˜ 1 and X ˜ 2 , the errors of X ˜ 1 and X ˜ 2 can be approximated as follows:

‖ X ˜ i − e A ˜ ‖ ≈ ε ˜ i ≔ ‖ X ˜ i − X ˜ 3 ‖ ( i = 1 , 2 ) .

By solving ε ˜ i = γ exp ( − ρ ⁄ h i ) ( i = 1 , 2 ) for ρ and γ , we obtain

ρ = h 1 h 2 h 1 − h 2 log ε ˜ 1 ε ˜ 2 , γ = ε ˜ 1 exp ρ h 1 ,

and we can approximate ‖ X ˜ 3 − e A ˜ ‖ ≈ ε ˜ 3 ≔ γ exp ( − ρ ⁄ h 3 ) . Because ε ˜ 3 is just an approximation, we select a safety parameter η > 0 and stop the algorithm if ε ˜ 3 ≤ η ε ˜ to obtain X ˜ 3 ( ≈ e A ˜ ) . If ε ˜ 3 > η ε ˜ , we set h 4 = ρ ⁄ log ( γ ⁄ η ε ˜ ) and evaluate the DE formula (4) once more. When ε ˜ 1 ≲ ε ˜ 2 , it means that the initial mesh size is too large to assume the exponentially convergence of the DE formula (4). For this case, we set h i ← h i + 1 ( i = 1 , 2 ) and repeat the procedure. The detail of the algorithm is given in Algorithm 2.

Algorithm 2 Automatic quadrature algorithm for e A based on the DE formula
	Input A ∈ C n × n , tolerance for the error ε > 0 , shift parameter σ < 0 , initial mesh size h 1 > 0 , safety parameter η > 0 , the minimum mesh size h min > 0
	Output X ≈ e A
1 :	Compute the rightmost eigenvalue of A , λ right
2 :	A ˜ = A + ( σ − λ right ) I , ε ˜ = ε ⁄ ∣ e λ right − σ ∣
3 :	h 2 = h 1 ⁄ 2 , h 3 = h 1 ⁄ 4
4:	l i , r i = GetInterval ( ε ˜ ⁄ 2 , h i , A ˜ ) ( i = 1 , 2 , 3 ) ⊳ Set ε ˜ ⁄ 2 for the input of `GetInterval` to bound the truncation error by ε ˜ ⁄ 2
5:	X ˜ i = h i ∑ k = l i r i F h i ( k h i , A ˜ ) ( i = 1 , 2 , 3 )
6:	ε ˜ i = ‖ X ˜ i − X ˜ 3 ‖ 2 ( i = 1 , 2 )
7:	ρ = h 1 h 2 log ( ε ˜ 1 ⁄ ε ˜ 2 ) ⁄ ( h 1 − h 2 ) , γ = ε ˜ 1 exp ( ρ ⁄ h 1 )
8:	ε ˜ 3 = γ exp ( − ρ ⁄ h 3 )
9:	if ε ˜ 3 < ε ˜ ⁄ η then
10:	return X = e λ right − σ X ˜ 3
11:	else
12:	h 4 = ρ ⁄ log ( γ η ⁄ ε ˜ )
13:	if ε ˜ 1 ≲ ε ˜ 2 or h 4 < h min then
14:	Set h i ← h i + 1 ( i = 1 , 2 , 3 ) and X ˜ i ← X ˜ i + 1 ( i = 1 , 2 )
15:	l 3 , r 3 = GetInterval ( ε ˜ ⁄ 2 , h 3 , A ˜ )
16:	X ˜ 3 = h 3 ∑ k = l 3 r 3 F h 3 ( k h 3 , A ˜ )
17:	go to Step 6.
18:	else
19:	l 4 , r 4 = GetInterval ( ε ˜ ⁄ 2 , h 4 , A ˜ )
20:	X ˜ 4 = h 4 ∑ k = l 4 r 4 F h 4 ( k h 4 , A ˜ )
21:	return X = e λ right − σ X ˜ 4
22:	end if
23:	end if

4 Numerical examples

The computation is carried out with Julia 1.10.0 on an Intel Core i5-9600K CPU with 32 GB RAM. The IEEE double-precision arithmetic is used unless otherwise stated. The reference solutions in Section 4.5 are computed by [ 13 ⁄ 13 ] -type Padé approximation with the BigFloat data type of Julia, which gives roughly 77 significant decimal digits. The Padé approximation is performed after scaling so that the 1 norm of the matrix is less than 5.37. The programs for these experiments are available on GitHub^[1]. For the DE formula, we use the change of variable in (3), and the infinite sums in GetInterval are truncated to 50 terms.

4.1 Computing the scalar exponential with the DE formula

By using the result in [3], the upper bound on the error of the DE formula can be obtained by ‖ e A − h ∑ k = l r F h ( k h , A ) ‖ 2 ≤ ( 1 + 2 ) max z ∈ W ( A ) ∣ e z − h ∑ k = l r F h ( k h , z ) ∣ , where W ( A ) ≔ { x H A x : x ∈ C n , ‖ x ‖ 2 = 1 } is the field-of-value of A . Hence, it is worth observing the error of the DE formula for the scalar exponential ∣ e z − h ∑ k = l r F h ( k h , z ) ∣ . Figure 2 visualizes the error, and we can notice the followings:

When z is on or near the negative real axis, the error of the DE formula with h = 0.1 is about 1 0 − 16 .
When − Re ( z ) is large, the error of the DE formula with h = 0.2 is about 1 0 − 16 .
When z is in a sector region { z ∈ C : ∣ arg ( − z ) ∣ < π ⁄ 4 } , the error of the DE formula with h = 0.05 is about 1 0 − 16 .

From these findings, when A is a near Hermitian matrix, the error in 2-norm of the DE formula for e A with h = 0.1 will be about 1 0 − 16 .

$Figure 2 The error of the double exponential formula for the scalar exponential e z {{\rm{e}}}^{z} . The mesh size is set to h = 0.2 , 0.1 , 0.05 h=0.2,0.1,0.05 . The left column shows the error for the area { z ∈ C : ‒ 30 ≤ Re ( z ) ≤ 10 , ‒ 20 ≤ Im ( z ) ≤ 20 } \left\{z\in {\mathbb{C}}:‒30\le {\rm{Re}}\left(z)\le 10,‒20\le {\rm{Im}}\left(z)\le 20\right\} and the right one shows the error for { z ∈ C : ‒ 5,000 ≤ Re ( z ) ≤ 0 , ‒ 2,500 ≤ Im ( z ) ≤ 2,500 } \left\{z\in {\mathbb{C}}:‒\hspace{0.1em}\text{5,000}\hspace{0.1em}\le {\rm{Re}}\left(z)\le 0,‒\hspace{0.1em}\text{2,500}\hspace{0.1em}\le {\rm{Im}}\left(z)\hspace{0.33em}\le \hspace{0.33em}\hspace{0.1em}\text{2,500}\hspace{0.1em}\right\} .$

Figure 2

The error of the double exponential formula for the scalar exponential e z . The mesh size is set to h = 0.2 , 0.1 , 0.05 . The left column shows the error for the area { z ∈ C : ‒ 30 ≤ Re ( z ) ≤ 10 , ‒ 20 ≤ Im ( z ) ≤ 20 } and the right one shows the error for { z ∈ C : ‒ 5,000 ≤ Re ( z ) ≤ 0 , ‒ 2,500 ≤ Im ( z ) ≤ 2,500 } .

4.2 Accuracy dependence on the shift parameter σ

Algorithm 1 requires the parameter σ for the input. For reference, we see the accuracy dependence of the DE formula on σ . We consider two test matrices A k = Z D k Z − 1 ( k = 1 , 2 ) , where Z ∈ C 50 × 50 is generated by using a function randsvd in MatrixDepot.jl [25] so that κ ( Z ) = 1 0 2 , and D k ∈ C 50 × 50 is a diagonal matrix whose diagonals are d i = 1 − 1 0 2 k ( i − 1 ) ⁄ 49 + i ν i ⁄ 20 ( i = 1 , … , 50 , ν i ∼ N ( 1 , 0 ) ) .

Figure 3 shows the error of the DE formula for σ ∈ [ − 10 , 5 ] . When σ ≥ 0 , the error of the DE formula is large for both matrices because A does not satisfy the condition for (1). Thus, σ must be smaller than 0. On the other hand, as σ becomes small, the computational result becomes inaccurate. Hence, it would be better to set σ ∈ [ − 5 , 0 ) .

$Figure 3 Accuracy dependence of the double exponential formula on the shift parameter σ \sigma .$

Figure 3

Accuracy dependence of the double exponential formula on the shift parameter σ .

4.3 Accuracy control of automatic quadrature

Here, we see the accuracy control of Algorithm 2. Figure 4 shows the error of Algorithm 2, where the test matrices are the same as in Section 4.2, and the safe parameter is set to η = 2 . We note that these test matrices have field-of-value that crosses the imaginary axis: max z ∈ W ( A 1 ) Re ( z ) ≈ 612 , and max z ∈ W ( A 2 ) Re ( z ) ≈ 61,403 .

$Figure 4 The accuracy of Algorithm 2. The left figure shows the accuracy dependence on the input tolerance ε \varepsilon and the error of the computational results. The right figure shows the error and its presumed value of the case when A = A 1 A={A}_{1} and ε = 1 0 ‒ 10 \varepsilon =1{0}^{‒10} . The horizontal axis is the inverse of the selected mesh size, and the vertical axis is the error.$

Figure 4

The accuracy of Algorithm 2. The left figure shows the accuracy dependence on the input tolerance ε and the error of the computational results. The right figure shows the error and its presumed value of the case when A = A 1 and ε = 1 0 ‒ 10 . The horizontal axis is the inverse of the selected mesh size, and the vertical axis is the error.

The left figure shows the accuracy dependence on the input tolerance ε , and it shows that the accuracy of Algorithm 2 could be controlled by ε . The right figure shows the error and its presumed value of the computational results for the selected mesh size for the case, where A = A 1 and ε = 1 0 − 10 . From the figure, it is found that the algorithm could presume the error.

4.4 Accuracy of Algorithm 2

In this example, we compute the exponential of 51 test matrices from MatrixDepot.jl [25]. The size of these matrices is 10 × 10 , and they are shifted as A ′ = A − λ right I before the computation. We computed e A ′ by using Algorithm 2 with ε = 1 0 − 8 , 1 0 − 16 . For reference, we also computed e A ′ with exp function in Julia, which uses the scaling and squaring algorithm [11].

Figure 5 shows the error of this example. It is observed that the error of Algorithm 2 is less than 1 0 − 8 for most matrices when ε = 1 0 − 8 . This indicates that Algorithm 2 could control the error. In addition, when ε = 1 0 − 16 , the accuracy of Algorithm 2 is comparable to that of exp function.

$Figure 5 The error ‖ X ‒ e A ′ ‖ 2 \Vert X‒{{\rm{e}}}^{A^{\prime} }{\Vert }_{2} of Algorithm 2 and exp function in Julia, which uses the scaling and squaring algorithm in [11].$

Figure 5

The error ‖ X ‒ e A ′ ‖ 2 of Algorithm 2 and exp function in Julia, which uses the scaling and squaring algorithm in [11].

4.5 A comparison of the DE formula with the Cauchy integral for non-Hermitian matrices

In this example, we compute e A for non-Hermitian matrices by using Algorithm 1 (denoted by DE) and an algorithm based on the Cauchy integral.

The matrices are generated from the convection-diffusion problems

∂ u ∂ t = d Δ u − c ⊤ ∇ u in Ω = ( 0 , 1 ) 2 , u = 0 on ∂ Ω , u ( 0 , x ) = u 0 ( x ) ,

where d = 0.001 , 0.01 , and c = [ 0.2 , 0.2 ] , [ 0.4 , 0.4 ] . By using the finite element method, we have a linear evolution equation M u ′ ( t ) = K u ( t ) whose solution is exp ( t M − 1 K ) u ( 0 ) . We use FreeFEM++ [10] for the discretization, and the eigenvalues of M − 1 K is illustrated in Figure 6.

$Figure 6 Eigenvalues of the test matrix M ‒ 1 K {M}^{‒1}K in Section 4.5.$

Figure 6

Eigenvalues of the test matrix M ‒ 1 K in Section 4.5.

For DE, we set σ = − 2.5 and ε ≈ 2.2 × 1 0 − 16 . For the Cauchy integral:

e A = 1 2 π i ∫ Γ e z ( z I + A ) − 1 d z ,

where Γ is a contour enclosing the eigenvalue of A , we used the algorithm (denoted by Talbot) proposed in [5]. Talbot employs the midpoint rule on the Talbot contour

Γ = m ( − σ ˜ + μ ˜ θ cot ( α ˜ θ ) + ν ˜ i θ ) : − π 2 ≤ θ ≤ − π 2 ,

where m is the number of abscissas, and σ ˜ , μ ˜ , α ˜ , ν ˜ is selected to minimize the error of the midpoint rule. Although Talbot is not designed for non-Hermitian matrices, it is not entirely inapplicable to them. Because of the optimality of the integral path of Talbot, its convergence is faster than that of DE, provided that the eigenvalues do not cross the contour. Here, we intend to show an example illustrating the limitations of Talbot and its performance compared to DE.

Figure 7 illustrates the convergence histories of these algorithms for t = 1 . For d = 0.01 , both DE and Talbot give accurate results, and Talbot converges faster than DE. For d = 0.001 , DE still gives accurate results while Talbot does not converge. Hence, our algorithm can be a choice of quadrature-based algorithms for non-Hermitian matrices having non-trivial eigenvalue distribution.

$Figure 7 Convergence histories of DE and Talbot for the problem exp ( M ‒ 1 K ) \exp \left({M}^{‒1}K) in Section 4.5.$

Figure 7

Convergence histories of DE and Talbot for the problem exp ( M ‒ 1 K ) in Section 4.5.

5 Conclusion

The DE formula was considered for the computation of e A . To utilize the DE formula, we analyzed the truncation error and proposed algorithms. Numerical results showed the validity of the algorithm.

Future work includes improvement of the change of variable of the DE formula for e A and comprehensive performance evaluations of the proposed algorithms on parallel computers for practical problems.

Acknowledgements

The authors are grateful to the two anonymous referees for the careful reading and the comments that substantially enhanced the quality of the manuscript.

Funding information: This work was supported by JSPS KAKENHI Grant Numbers 20H00581 and 20K20397. The authors would like to thank Enago (www.enago.jp) for English language editing.
Author contributions: Fuminori Tatsuoka: conceptualization, experimentation, writing – original draft; Tomohiro Sogabe: writing – review and editing, supervision; Tomoya Kemmochi: writing – review and editing, supervision; Shao-Liang Zhang: writing – review and editing, supervision.
Conflict of interest: The authors declare no conflicts of interest.
Data availability statement: All data that support this study are available at GitHub https://github.com/f-ttok/article-expmde.

Appendix A Derivation of the integral representation (1)

From [6, p. 758], it holds that

(A1) e − θ b 1 ⁄ 2 − 1 b = ∫ − ∞ 0 sin ( θ ( − u ) 1 ⁄ 2 ) π u 1 b − u d u ,

where b 1 ⁄ 2 is the principal square root. Because (A1) is based on the Cauchy integral for ( e − θ b 1 ⁄ 2 − 1 ) ⁄ b , see the proof of [6, Prop. 1], we have

B − 1 ( e − θ B 1 ⁄ 2 − I ) = ∫ − ∞ 0 sin ( θ ( − u ) 1 ⁄ 2 ) π u ( B − u I ) − 1 d u ,

where all the eigenvalues of B 1 ⁄ 2 are in the open right half plane. By differentiating both sides twice with respect to θ and then substituting θ = 1 , we have

(A2) e − B 1 ⁄ 2 = 1 π ∫ − ∞ 0 sin ( ( − u ) 1 ⁄ 2 ) ( B − u I ) − 1 d u = 2 π ∫ 0 ∞ x sin ( x ) ( x 2 I + B ) − 1 d x ,

where u = − x 2 . We have (1) by substituting A = − B 1 ⁄ 2 into (A2).

B Derivative of the change of variable

Let v ( t ) = − 2 t − α ( 1 − e − t ) − β ( e t − 1 ) . Then, the transformation (3) for the DE formula is x h ( t ) = π t ⁄ h ( 1 − exp ( v ( t ) ) ) and its derivative is

x h ′ ( t ) = π h 1 − e v ( t ) + t v ′ ( t ) e v ( t ) ( e v ( t ) − 1 ) 2 .

Because exp ( v ( t ) ) → 1 as t → 0 , we would use x h ( 0 ) = π ⁄ h ( α + β + 2 ) and

x h ′ ( 0 ) = π 2 h α 2 + 2 α β + 5 α + β 2 + 3 β + 4 α 2 + 2 α β + 4 α + β 2 + 4 β + 4 .

When v ( t ) ≤ log ( 1 − 1 ⁄ 2 ) ( < 0 ) , the transformation (3) holds the inequality x h ′ ( t ) < 2 π ⁄ h that is used in Proposition 1. First, we note that v ( t ) monotonically decreases for t > 0 because v ′ ( t ) = − 2 − α e − t − β e t < 0 , where α > 0 and β > 0 . Therefore, for t ≥ t 0 , where t 0 ( > 0 ) be such that v ( t 0 ) = log ( 1 − 1 ⁄ 2 ) , it leads 1 ⁄ ( 1 − e v ( t ) ) 2 ≤ 2 . Finally, we have

h π x h ′ ( t ) = 1 − e v ( t ) + t v ′ ( t ) e v ( t ) ( 1 − e v ( t ) ) 2 ≤ 2 ( 1 − e v ( t ) + t v ′ ( t ) e v ( t ) ) < 2 .

Experimentally, the value t 0 will be about 0.5 as in Table A1, and the right truncation point t = r h is larger than 3. Hence, the transformation (3) will satisfy the assumption of Proposition 1 in practical cases.

Table A1

The value of t 0 , which is the solution of v ( t ) = log ( 1 ‒ 1 ⁄ 2 ) , for several h

h	t 0
1 0 0	0.493
1 0 ‒ 1	0.514
1 0 ‒ 2	0.524
1 0 ‒ 3	0.526

References

[1] A. H. Al-Mohy and N. J. Higham, A new scaling and squaring algorithm for the matrix exponential, SIAM J. Matrix Anal. Appl. 31 (2009), no. 3, 970–989. 10.1137/09074721XSuche in Google Scholar

[2] A. H. Al-Mohy and N. J. Higham, Computing the action of the matrix exponential, with an application to exponential integrators, SIAM J. Sci. Comput. 33 (2011), no. 2, 488–511. 10.1137/100788860Suche in Google Scholar

[3] M. Crouzeix and C. Palencia, The numerical range is a (1+2)-spectral set, SIAM J. Matrix Anal. Appl. 38 (2017), no. 2, 649–655. 10.1137/17M1116672Suche in Google Scholar

[4] O. De la Cruz Cabrera, M. Matar, and L. Reichel, Analysis of directed networks via the matrix exponential, J. Comput. Appl. Math. 355 (2019), 182–192. 10.1016/j.cam.2019.01.015Suche in Google Scholar

[5] B. Dingfelder and J. A. C. Weideman, An improved Talbot method for numerical Laplace transform inversion, Numer Algor, 68 (2015), no. 1, 167–183. 10.1007/s11075-014-9895-zSuche in Google Scholar

[6] V. Druskin and L. Knizhnerman, Extended Krylov subspaces: Approximation of the matrix square root and related functions, SIAM J. Matrix Anal. Appl. 19 (1998), no. 3, 755–771. 10.1137/S0895479895292400Suche in Google Scholar

[7] M. Fasi and N. J. Higham, An arbitrary precision scaling and squaring algorithm for the matrix exponential, SIAM J. Matrix Anal. Appl. 40 (2019), no. 4, 1233–1256. 10.1137/18M1228876Suche in Google Scholar

[8] T. Göckler and V. Grimm, Uniform approximation of ϕ-functions in exponential integrators by a rational krylov subspace method with simple poles, SIAM J. Matrix Anal. Appl. 35 (2014), no. 4, 1467–1489. 10.1137/140964655Suche in Google Scholar

[9] S. Güttel and Y. Nakatsukasa, Scaled and squared subdiagonal Padé approximation for the matrix exponential, SIAM J. Matrix Anal. Appl. 37 (2016), no. 1, 145–170. 10.1137/15M1027553Suche in Google Scholar

[10] F. Hecht, New development in FreeFem++, J. Numer. Math. 20 (2012), no. 3–4, 251–265. 10.1515/jnum-2012-0013Suche in Google Scholar

[11] N. J. Higham, The scaling and squaring method for the matrix exponential revisited, SIAM J. Matrix Anal. Appl. 26 (2005), no. 4, 1179–1193. 10.1137/04061101XSuche in Google Scholar

[12] N. J. Higham, Functions of Matrices: Theory and Computation, SIAM, Philadelphia, PA, 2008. 10.1137/1.9780898717778Suche in Google Scholar

[13] M. Hochbruck and A. Ostermann, Exponential integrators, Acta Numer. 19 (2010), 209–286. 10.1017/S0962492910000048Suche in Google Scholar

[14] C. Moler and C. Van Loan, Nineteen dubious ways to compute the exponential of a matrix, SIAM Review 20 (1978), no. 4, 801–836. 10.1137/1020098Suche in Google Scholar

[15] C. Moler and C. Van Loan, Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later, SIAM Review 45 (2003), no. 1, 3–49. 10.1137/S00361445024180Suche in Google Scholar

[16] T. Ooura and M. Mori, The double exponential formula for oscillatory functions over the half infinite interval, J. Comput. Appl. Math. 38 (1991), no. 1, 353–360. 10.1016/0377-0427(91)90181-ISuche in Google Scholar

[17] T. Ooura and M. Mori, A robust double exponential formula for Fourier-type integrals, J. Comput. Appl. Math. 112 (1999), no. 1–2, 229–241. 10.1016/S0377-0427(99)00223-XSuche in Google Scholar

[18] Y. Saad, Analysis of some Krylov subspace approximations to the matrix exponential operator, SIAM J. Numer. Anal. 29 (1992), no. 1, 209–228. 10.1137/0729014Suche in Google Scholar

[19] T. Schmelzer and L. N. Trefethen, Evaluating matrix functions for exponential integrators via Carathéodory-Fejér approximation and contour integrals, Electron. Trans. Numer. Anal. 29 (2007), 1–18. Suche in Google Scholar

[20] H. Takahasi and M. Mori, Double exponential formulas for numerical integration, Publ. Res. Inst. Math. Sci. 9 (1974), 721–741. 10.2977/prims/1195192451Suche in Google Scholar

[21] F. Tatsuoka, T. Sogabe, Y. Miyatake, T. Kemmochi, and S.-L. Zhang, Computing the matrix fractional power with the double exponential formula, Electron. Trans. Numer. Anal. 54 (2021), 558–580. 10.1553/etna_vol54s558Suche in Google Scholar

[22] F. Tatsuoka, T. Sogabe, Y. Miyatake, and S.-L. Zhang, Algorithms for the computation of the matrix logarithm based on the double exponential formula, J. Comput. Appl. Math. 373 (2020), 112396. 10.1016/j.cam.2019.112396Suche in Google Scholar

[23] L. N. Trefethen and J. A. C. Weideman, The exponentially convergent trapezoidal rule, SIAM Rev. 56 (2014), no. 3, 385–458. 10.1137/130932132Suche in Google Scholar

[24] J. A. C. Weideman and L. N. Trefethen, Parabolic and hyperbolic contours for computing the Bromwich integral, Math. Comp. 76 (2007), no. 259, 1341–1356. 10.1090/S0025-5718-07-01945-XSuche in Google Scholar

[25] W. Zhang and N. J. Higham, Matrix depot: An extensible test matrix collection for Julia, PeerJ Comput. Sci. 2 (2016), e58. 10.7717/peerj-cs.58Suche in Google Scholar

Received: 2024-01-17

Revised: 2024-05-27

Accepted: 2024-05-27

Published Online: 2024-10-26

This work is licensed under the Creative Commons Attribution 4.0 International License.

Artikel in diesem Heft

https://doi.org/10.1515/spma-2024-0013

Schlagwörter für diesen Artikel

matrix function; matrix exponential; numerical quadrature; double exponential formula

Creative Commons

BY 4.0