Electronic Transactions on Numerical Analysis.
Volume 54, pp. 558–580, 2021.
Copyright c 2021, Kent State University.
ISSN 1068–9613.
DOI: 10.1553/etna_vol54s558

ETNA
Kent State University and
Johann Radon Institute (RICAM)COMPUTING THE MATRIX FRACTIONAL POWER WITH THE DOUBLEEXPONENTIAL FORMULA∗
FUMINORI TATSUOKA†, TOMOHIRO SOGABE†, YUTO MIYATAKE‡, TOMOYA KEMMOCHI†, AND
SHAO-LIANG ZHANG†
Abstract. Two quadrature-based algorithms for computing the matrix fractional power Aα are presented in thispaper. These algorithms are based on the double exponential (DE) formula, which is well-known for its effectiveness incomputing improper integrals as well as in treating nearly arbitrary endpoint singularities. The DE formula transforms
a given integral into another integral that is suited for the trapezoidal rule; in this process, the integral interval istransformed into an infinite interval. Therefore, it is necessary to truncate the infinite interval to an appropriate
finite interval. In this paper, a truncation method, which is based on a truncation error analysis specialized to thecomputation of Aα , is proposed. Then, two algorithms are presented—one where Aα is computed with a fixed
number of abscissa points and one with Aα computed adaptively. Subsequently, the convergence rate of the DE
formula for Hermitian positive definite matrices is analyzed. The convergence rate analysis shows that the DE formula
converges faster than Gaussian quadrature when A is ill-conditioned and α is a non-unit fraction. Numerical results
show that our algorithms achieve the required accuracy and are faster than other algorithms in several situations.
Key words. matrix function, matrix fractional power, numerical quadrature, double exponential formula
AMS subject classifications. 65F60, 65D30

1. Introduction. We consider the computation of the matrix power Aα for a matrix
A ∈ Cn×n and a real number α ∈ R. The matrix exponential is defined asexp(X) = I + X + X 2 /2! + X 3 /3! + · · ·(X ∈ Cn×n ),and the matrix logarithm is defined as an inverse function of the matrix exponential, i.e., any
matrix Y ∈ Cn×n satisfying exp(Y ) = A is a logarithm of A. If all the eigenvalues of Y lie
in the strip {z ∈ C : |Im(z)| < π}, then Y is called the principal logarithm and is denoted bylog(A). Then, the principal matrix power Aα is defined as follows:Aα = exp(α log(A)).
It is known that Aα exists and is unique when all the eigenvalues of A lie in the set{z ∈ C : z ∈/ (−∞, 0]}; see, e.g., [11]. The matrix fractional power arises in several situations in computational science, e.g., fractional differential equations [2, 5, 16], lattice QCD
calculation [4], and the computation of the weighted matrix geometric mean [8].
In this study, we compute Aα using the double exponential (DE) formula [17] for the
integral representation:Z ∞sin(απ)α(1.1)
A(t1/α I + A)−1 dt(0 < α < 1);A =απ0∗ Received December 2, 2020. Accepted June 2, 2021. Published online on September 6, 2021. Recommended
by Dario Bini. This paper is an augmented version of our previous report (F. Tatsuoka, T. Sogabe, Y. Miyatake, and
S.-L. Zhang, A note on computing the matrix fractional power using the double exponential formula, Trans. Jpn.
Soc. Ind. Appl. Math., 28 (2018), pp. 142–161), which is in Japanese. There are two differences between thispaper and the previous report. One is that this paper presents the convergence rate analysis of the double exponential(DE) formula, while the previous report does not. The other is that this work employs a new transformationt(x) = exp(απ sinh(x)/2), while the previous one used the transformation t(x) = exp(α sinh(x)). Several
numerical results have shown that the transformation employed in this work is better than the previous one in terms of
convergence speed.
† Department of Applied Physics, Graduate School of Engineering, Nagoya University, Furo-cho, Chikusa-ku,
Nagoya 464-8603, Japan ({f-tatsuoka,sogabe,kemmochi,zhang}@na.nuap.nagoya-u.ac.jp).
‡ Cybermedia Center, Osaka University, 1-32 Machikaneyama, Toyonaka, Osaka 560-0043, Japan(miyatake@cas.cmc.osaka-u.ac.jp).

558

ETNA
Kent State University and
Johann Radon Institute (RICAM)COMPUTING Aα WITH THE DE FORMULA

559

see, e.g., [11, Eq. (1.4)]. Without loss of generality, we assume 0 < α < 1 in this work because
for any α ∈ R, we can compute Aα = Abαc Abαc−α , where bαc ∈ N and bαc − α ∈ [0, 1).
The DE formula transforms the integral (1.1) into another integral suited for the trapezoidal
rule. In particular, the transformed integrand decays double exponentially, and the transformedinterval is the infinite interval (−∞, ∞); see, Section 2. Hence, it is necessary to truncate theinfinite interval to a finite interval appropriately to utilize the DE formula.
In this paper, we first present an algorithm for computing Aα using the m-point DE
formula. To develop this algorithm, we analyze the truncation error of the DE formula for (1.1)
and propose a method for truncating the infinite interval to a finite interval within a given
tolerance. Then, we present an adaptive quadrature algorithm based on the first algorithm toachieve the required accuracy. This adaptive quadrature algorithm employs an a posteriori
error estimation technique discussed in [6]. Subsequently, we analyze the convergence rateof the DE formula for Hermitian positive definite (HPD) matrices and compare it with thatof Gaussian quadrature. The convergence rate analysis gives us an a priori estimate of thediscretization error similarly as that for Gaussian quadrature in [8]. Hence, for HPD matrices,
the number of required abscissa points can be estimated on the basis of this error estimation.
Several computational methods for Aα have been proposed. Examples are the Schur-Padéalgorithm [11, 12], the Schur logarithmic algorithm [13], and quadrature-based algorithms [2,6, 8, 9]. Among them, quadrature-based algorithms have two advantages. One advantage isthat these algorithms can compute Aα b (b ∈ Rn ) without computing Aα itself when A islarge and either sparse or structured. The other is that these algorithms can be parallelizedeasily in the sense that the integrand can be computed independently at each abscissa point.
See [18, Sect. 18] for more details.
Quadrature-based algorithms were developed based on two integral representations of
Aα . One representation utilizes the Cauchy integral [11, Eq. (1.1)], which has been thebasis for several algorithms [9]. These algorithms are specialized for the case where all
the eigenvalues of A lie on or near the positive real axis. Their convergence rate is better
than that of other quadrature-based algorithms, but they require complex arithmetic even
if A is a real matrix. To the best of our knowledge, for a general matrix A, no algorithms
based on the Cauchy integral are available because of the difficulty in selecting the integral
path. The other integral representation is (1.1). Unlike the Cauchy integral, the integral
representation (1.1) allows us to avoid selecting an integral path and complex arithmetic forreal matrices. In [2, 6, 8], the integral (1.1) was transformed into other integrals on [−1, 1],
and the transformed integrals were computed with Gaussian quadrature. For example, thetransformations t(u) = (1 + u)/(1 − u) and t(v) = (1 − v)α /(1 + v)α were considered.
There are two motivations for considering the DE formula. First, the DE formula may
converge faster than Gaussian quadrature. The DE formula is known for its effectiveness indealing with integrals on (half) infinite intervals as well as integrals with endpoint singularities;see, e.g., [14, 18]. Hence, the DE formula is often used when the convergence of Gaussian
quadrature is slow. Indeed, the convergence of Gaussian quadrature for (1.1) will be slow
when the condition number of A is large and α is a non-unit fraction. This slow convergenceis indicated by our numerical results and the analysis for HPD matrices in [8]. Second, theDE formula can be used as an efficient adaptive quadrature algorithm. One of the problemsof Gaussian quadrature is that it incurs high computational costs for an a posteriori error
estimate. For HPD matrices, one can use an a priori estimate technique that computes
λα for λ being the extreme eigenvalues of A. This estimation incurs lower costs than thecomputation of the integrand in (1.1); see, [8]. For general matrices, one can use an a posteriori
estimate technique that compares two approximations of Aα with different numbers of abscissa
points [6]. However, this estimation incurs high computational costs for Gaussian quadrature,

ETNA
Kent State University and
Johann Radon Institute (RICAM)

560

F. TATSUOKA, T. SOGABE, Y. MIYATAKE, T. KEMMOCHI, AND S.-L. ZHANG

because the integrand must be computed at all the abscissa points for a new approximation
with a different number of points. In contrast to Gaussian quadrature, the DE formula employs
the trapezoidal rule. The trapezoidal rule can reuse an approximation at a mesh of size h forcomputing a new approximation at a mesh of size h/2. Hence, the DE formula reduces thecomputational cost.
The rest of this paper is organized as follows. In the last of this section, a brief reviewof the DE formula is given. In Section 2, the truncation error of the DE formula for Aα
is analyzed, and the truncation method and the algorithms are presented. In Section 3, theconvergence rate of the DE formula is analyzed and compared with that of Gaussian quadrature.
Numerical results are shown in Section 4, and conclusions are given in Section 5.
Notation: Unless otherwise stated, throughout this paper, k · k denotes a consistentmatrix norm, e.g., the p-norm or the Frobenius norm, while k · k2 represents the 2-norm.
The symbol κ(A) indicates the condition number, i.e., κ(A) = kAk2 kA−1 k2 , and ρ(A) isthe spectral radius of A, i.e., |λmax | for λmax the largest absolute eigenvalue of A. The
inverse of sinh is denoted by asinh. The strip of width 2d about the real axis is denoted by
Dd := {z ∈ C : |Im(z)| < d}.
1.1. The DE formula. The DE formula exploits the fact that the trapezoidal rule for the
integrals of analytic functions on the real line converges exponentially; see, e.g., [14, 17].
Rb
The procedure based on the DE formula for a given integral a f (t) dt, where f is a scalarfunction, is as follows:
1. Apply a change of variables t = t(x) to the integral, where the transformed intervalis the infinite interval (−∞, ∞) and the transformed integrand t0 (x)f (t(x)) decays
double exponentially as x → ±∞.1
2. Truncate theR r infinite interval to a finite interval [l, r].
3. Compute l t0 (x)f (t(x)) dx by the trapezoidal rule.
If f is analytic on (a, b) and the change of variables is appropriate, then the transformed
integrand t0 (x)f (x) is analytic on (−∞, ∞) even if f (t) is not analytic at t = a, b. Thus, theDE formula can treat nearly arbitrary endpoint singularities. In addition, the transformation
with double exponential decay of the transformed integrand has a certain optimality property;
see, e.g., [14, 15]. Therefore, the DE formula can be used for computing integrals on the (half)
infinite interval.
For the integral (1.1), one can apply the double exponential formula. The integrand in (1.1)
is a matrix function of A in the sense of [10, Def. 1.2], and therefore we have the explicit
expressions of each element at hand. Considering the Jordan decomposition, each element∂kαof the integrand is 0 or of the form ∂t+ λ)−1 for k ≥ 0 and λ being the eigenvaluesk (tof A. Thus, with transformations such as t(x) = exp(απ sinh(x)/2), each element of theintegrand in (1.1) decays double exponentially as x → ∞, and one can compute integrals
simultaneously.
2. Algorithms based on the DE formula. This section presents two algorithms for computing Aα utilizing the DE formula. Here, the change of variables t(x) = exp(απ sinh(x)/2)
is considered for the DE formula. Thus, the algorithms are based on the integral representation
as follows:
Z ∞
α
A =
FDE (x) dx,−∞
1 More precisely, the expression “decays double exponentially” means that t0 (x)f (t(x)) is a function in the
Hardy space H∞ (Dd , ω), where ω is a function satisfying the conditions in [15, Thm. 3.2] and H∞ (Dd , ω) is afunction space defined in [15, Sect. 2].

ETNA
Kent State University and
Johann Radon Institute (RICAM)COMPUTING Aα WITH THE DE FORMULA

561

where
FDE (x) = t0 (x)F (t(x)),

sin(απ)A(t1/α I + A)−1 .
απ

F (t) =

The first algorithm computes Aα based on the m-point DE formula. To control the truncation
error, the algorithm selects a finite interval [l, r] satisfyingZ ∞Z r
FDE (x) dx −

(2.1)
−∞


2

≤

FDE (x) dx
l

2

Rr

for a given tolerance . Then, l FDE (x) dx is computed by using the m-point trapezoidal
Rrrule. The second algorithm adaptively computes l FDE (x) dx until the discretization error issmaller than or equal to /2. Hence, the total error of the second algorithm is smaller than orequal to .
The analysis of the truncation error for a given finite interval is provided in Section 2.1.
In Section 2.2, we propose a truncation method for (2.1). We explain the algorithms in
Section 2.3.
2.1. Truncation error analysis. The bound for the truncation error reads as follows:
L EMMA 2.1. Let A have no eigenvalues on the closed negative real axis, and α ∈ (0, 1).
For a given interval [l, r], let a = exp(απ sinh(l)/2) and b = exp(απ sinh(r)/2). Fora ≤ (2kA−1 k)−α and b ≥ (2kAk)α , we haveZ ∞Z r
FDE (x) dx −

−∞

(2.2)

FDE (x) dxl≤sin(απ)[α + (1 + α)kIk]
sin(απ)(3 − 2α)kAk 1−1/αa+b.
απ(1 + α)π(1 − α)(2 − α)Proof. From the triangle inequality, the truncation error is bounded by
Z ∞
−∞

Z ∞

Z l

Z r
FDE (x) dx ≤FDE (x) dx −

(2.3)

−∞

l

FDE (x) dx .

FDE (x) dx +r

We prove Lemma 2.1 by showing thatZ l

Z a
FDE (x) dx =

(2.4)
−∞

F (t) dt ≤sin(απ)[α + (1 + α)kIk]aαπ(1 + α)

F (t) dt ≤sin(απ)(3 − 2α)kAk 1−1/αb.
π(1 − α)(2 − α)

0

and
Z ∞

Z ∞
FDE (x) dx =

(2.5)
r

b

First, we show (2.4). For all t ∈ [0, a] with a ≤ (2kA−1 k)−α , we have t ≤ (2kA−1 k)−α ,and therefore, kt1/α A−1 k ≤ 1/2 (< 1). By applying the Neumann series expansion to(t1/α I + A)−1 , we get(t1/α I + A)−1 = A−1 [I − (−t1/α A−1 )]−1 = A−1∞
X(−1)k A−k tk/α .

k=0

ETNA
Kent State University and
Johann Radon Institute (RICAM)

562

F. TATSUOKA, T. SOGABE, Y. MIYATAKE, T. KEMMOCHI, AND S.-L. ZHANG

Therefore, the integral can be rewritten as follows:
Z aA(t1/α I + A)−1 dt = A

0

Z a

A−1

"∞
X

0

=

∞
X

k

(−1) A−k(−1)k A−k tk/α dt

k=0

Z a
t

k/α


dt =

0

k=0

= aI +

#

∞
X

∞
X

k

(−1) A−kk=0(−1)k A−k



k=1

1
ak/α+1k/α + 1
1ak/α+1 .
k/α + 1Thus, we get
Z asin(απ)AαπF (t) dt =0

Z a(t1/α I + A)−1 dt

0



∞
X
1
sin(απ)(−1)k A−kaI +ak/α+1απk/α + 1
k=1
"

#
∞
X
1
sin(απ)
−1 k
k/α+1
A
a kIk +a.
≤απk/α + 1
=

(2.6)

k=1

Furthermore, from a ≤ (2kA−1 k)−α , we obtain
∞
X

A

−1 k

k=1

(2.7)



1
ak/α+1k/α + 1



∞
X
k
αa
a1/α A−1
k+αk=1∞
X αa  1 k
α
≤
a.
=
1+α 2
1+α

=

k=1

We arrive at (2.4) by substituting (2.7) in (2.6).
Next, we show (2.5). The outline of this proof is almost the same as that of (2.4). For allt ∈ [b, ∞) with b ≥ (2kAk)α , we have kt−1/α Ak ≤ 1/2 (< 1). By applying the Neumann
series expansion to (t1/α I + A)−1 , we get(t1/α I + A)−1 = t−1/α [I − (−t−1/α A)]−1 = t−1/α∞X(−1)k t−k/α Ak .

k=0

Therefore, the integral can be rewritten as
Z ∞
(t

1/α

I + A)

#
Z ∞"
∞
X
−1/α
k −k/α kdt =t(−1) t
A dt

−1

b

b

=

∞ 
X

k

(−1) A

k=0

k=0
k

Z ∞
t

−(k+1)/α

b

dt =

∞
X

(−1)k Ak

k=0
∞

(2.8)



α
b1−(k+1)/αk+1−αX
α 1−1/αα=bI+(−1)k Akb1−(k+1)/α .
1−α
k+1−α
k=1

ETNA
Kent State University and
Johann Radon Institute (RICAM)COMPUTING Aα WITH THE DE FORMULA

563

Thus,
Z ∞

Z ∞
sin(απ)A(t1/α I + A)−1 dt
απ
b
"
#
∞
X
sin(απ)α 1−1/αα1−(k+1)/α
k k
=
b
b
A+A
(−1) A
απ
1−α
k+1−α
k=1
"
#
∞
sin(απ)1 1−1/α X
1
1−(k+1)/α
k
≤
kAk
b
+
b
.
kAk
π
1−α
k+1−α

F (t) dt =

(2.9)
b

k=1

α

Moreover, from b ≥ (2kAk) , we find(2.10)

∞
X

kAkk

k=1

∞
X
1
1
b1−(k+1)/α = b1−1/α
kb−1/α Akk
k+1−α
k+1−α
k=1
 k
∞
X
1
1 1−1/α1=b.
≤ b1−1/α
2−α 2
2−α
k=1

We obtain (2.5) by substituting (2.10) into (2.9). In conclusion, by combining (2.3), (2.4), and(2.5), we arrive at (2.2).
R EMARK 2.2. The technique of applying the Neumann series expansion to the integrandwas considered in [6, Thm. 2.1], and the following upper bound was derived:
Z ∞2b1−1/α(t1/α I + A)−1 dt ≤.
1/α − 1bHence, the relation
Z ∞F (t) dt ≤(2.11)b

2 sin(απ)kAk 1−1/αbπ(1 − α)holds. Note that the upper bound (2.5) is smaller than (2.11) because
2 sin(απ)kAk 1−1/α sin(απ)(3 − 2α)kAk 1−1/αb−bπ(1 − α)π(1 − α)(2 − α)
3 − 2αsin(απ)kAk 1−1/αsin(απ)kAk 1−1/α 1
b
2−
b
> 0.
=
=
π(1 − α)2−απ(1 − α)2−αThis difference originates from the separation of the sum into two terms in (2.8).
2.2. A truncation method based on the error analysis. On the basis of Lemma 2.1, we
propose a method to truncate the infinite interval to a finite interval within the given tolerance
below.
P ROPOSITION 2.3. For given  > 0, let
2 log(b)2 log(a),r = asinh,l = asinh
απ
απ
where




πα(1 + α), (2kA−1 k2 )−α ,4 sin(απ)(1 + 2α)()α/(α−1) π(1 − α)(2 − α)α, (2kAk2 ).
b = max
4 sin(απ)(3 − 2α) kAk2
a = min

ETNA
Kent State University and
Johann Radon Institute (RICAM)

564

F. TATSUOKA, T. SOGABE, Y. MIYATAKE, T. KEMMOCHI, AND S.-L. ZHANG

Then, the following holds:
Z ∞Z r
FDE (x) dx −FDE (x) dx(2.12)
−∞

l

≤
2


.
2

Proof. From the definitions of a and b, the inequalities a ≤ (2kA−1 k2 )−α andb ≥ (2kAk2 )α are valid. Therefore, the interval [l, r] satisfies the assumptions for Lemma 2.1.
Hence, it follows thatZ ∞Z r
FDE (x) dx −FDE (x) dx−∞
l2(2.13)
sin(απ)(1 + 2α)
sin(απ)(3 − 2α)kAk2 1−1/α≤a+b.
απ(1 + α)π(1 − α)(2 − α)From the definition of a and b, we have
(2.14)

sin(απ)(1 + 2α)a≤ ,απ(1 + α)
4sin(απ)(3 − 2α)kAk2 1−1/αb≤ .
π(1 − α)(2 − α)4

We obtain (2.12) by substituting (2.14) into (2.13).
2.3. Algorithms. We first explain the algorithm for computing Aα based on the mpoint DE formula.
This algorithm computes [l, r] that satisfies (2.12) by Proposition 2.3 andRrcomputes l FDE (x) dx by the m-point trapezoidal rule. The details of the algorithm aregiven in Algorithm 1.
Algorithm 1 m-point DE formula for computing Aα .
Input A, α ∈ (0, 1),  > 0, m.
1: l, r = GetInterval(A, α, )−1
2: Set F̃DE (x) = exp(απ sinh(x)/2) cosh(x) [exp(π sinh(x)/2)I + A]
3: h = (r − l)/(m − 1)
Pm−2
4: T = h[F̃DE (l) + F̃DE (r)]/2 + h
k=1 F̃DE (l + kh)
Output sin(απ)AT /2 ≈ Aα
5: function GetInterval(A, α, )6:Compute kAk2 , kA−1 k2
7:
a1 = [απ(1 + α)]/[4 sin(απ)(1 + 2α)], a2 = (2kA−1 k2 )−α8:a = min{a1 , a2 }9:
b1 = [π(1−α)(2−α)]α/(α−1) /[4 sin(απ)(3−2α) kAk2 ]α/(α−1) ,10:
b = max{b1 , b2 }11:
l = asinh(2 log(a)/απ),
12:return l, r
13: end function

b2 = (2kAk2 )αr = asinh(2 log(b)/απ)

When the tolerance  is sufficiently small, the accurate computation of kAk2 , kA−1 k2 instep 6 of Algorithm 1 is not necessary because the errors originating from these values havelittle effect on the truncation error. Assume that  is sufficiently small and [a, b] = [a1 , b1 ]. Let∆ ∈ R be such that the computed 2-norm of A is equal to kAk2 (1 + ∆). Then, the computed
value of b without rounding errors isα/(α−1)π(1 − α)(2 − α)b=b,
4 sin(απ)(3 − 2α) kAk2

ETNA
Kent State University and
Johann Radon Institute (RICAM)COMPUTING Aα WITH THE DE FORMULA

565

where b = /(1 + ∆). Hence, the upper bound of the truncation error isZ ∞Z r
FDE (x) dx −

−∞

l



1
1
1
FDE (x) dx ≤ ( + b ) ≤  1 +.
2
2
1 − |∆|For example, when ∆ = 10−2 , the upper bound of the truncation error is at most 1.005.
Therefore, the effect of these errors will be negligible. In practice, the effect of the error can
˜
be canceled by setting  of Algorithm 1 to  = ˜/(1 + 1/(1 − |∆|)),
where ˜ is the tolerance
˜
for the truncation error of the DE formula and ∆ is the tolerance for the relative error used in
the computation of kAk2 .
It may be desirable to apply bounds for the relative truncation error rather than for theabsolute truncation error; in other words, one may wish to select [l, r] that satisfies
RrAα − l FDE (x) dx 2
˜
≤
α
kA k2
2
for a given tolerance ˜. In this case, one can find such an interval by setting  = ρ(A)α ˜
as input for Algorithm 1. This is because ρ(A)α = ρ(Aα ) bounds kAα k2 from below. The
relation ρ(A)α = ρ(Aα ) can be verified from the fact that the eigenvalues of Aα are λαi ,where λi are the eigenvalues of A; see, e.g., [10, Thm. 1.13].
Before explaining the adaptive quadrature algorithm, we consider an a posteriori estimation of the discretization error. Let m be the number of abscissa points, h = (r − l)/(m − 1)be the mesh size, and T (h) be the trapezoidal rule with the mesh size h:m−2

T (h) :=

Xh
(FDE (l) + FDE (r)) + h
FDE (l + kh).
2
k=1

Then, T (h/2) can be computed as
 


m−1
h
1
h X
h
T
= T (h) +FDE l + (2k − 1)
.
2
2
2
2
k=1

The upper bound for the discretization error of T (h/2) can be obtained by applying thediscussionR in [6, p. 424] to the trapezoidal rule.
R r Assume that h is small enough to satisfyrkT (h) − l FDE (x) dxk ≤ µ and kT (h/2) − l FDE (x) dxk ≤ cµ, with 0 < c ≤ 1/2. If
kT (h/2) − T (h)k ≤ µ̃, thenZ r
T (h) −FDE (x) dx ≤ T (h) − Tl   Z rhh−
FDE (x) dx ,
+ T
2
2
l

that is, µ ≤ µ̃ + cµ, or equivalently, µ ≤ µ̃/(1 − c). Hence,
Z r(2.15)
l

 
 
h
h
FDE (x) dx − T≤ T (h) − T.
22

Now we state the adaptive quadrature algorithm. First, the algorithm computes[l, r], asRr
described in Algorithm 1. Then, the algorithm repeats the computation of l FDE (x) dx byincreasing the number of abscissa points until the discretization error is smaller than /2. The
details of the algorithm are given in Algorithm 2.

ETNA
Kent State University and
Johann Radon Institute (RICAM)

566

F. TATSUOKA, T. SOGABE, Y. MIYATAKE, T. KEMMOCHI, AND S.-L. ZHANG

Algorithm 2 Adaptive quadrature algorithm for Aα based on the DE formula.
Input A, α ∈ (0, 1),  > 0, m0 ≥ 2.
1: l, r = GetInterval(A, α, ). GetInterval is defined in Algorithm 1.
−1
2: Set F̃DE (x) = exp(απ sinh(x)/2) cosh(x) [exp(π sinh(x)/2)I + A]
3: s = −1
4: h0 = (r − l)/(m0 − 1)
Pm−2
5: T0 = h0 [F̃DE (l) + F̃DE (r)]/2 + h0
k=1 F̃DE (l + kh)
6: repeat
7:
s=s+1
8:
hs+1 = hs /2
Pms −1
9:
Ts+1 = Ts /2 + hs+1 k=1
F̃DE (l + (2k − 1)hs+1 )10:
ms+1 = 2ms − 1
11: while sin(απ)kATs+1 − ATs k/2 > /2
12: T = Ts+1
Output sin(απ)AT /2 ≈ Aα
R EMARK 2.4. The computational cost of Algorithm 1 for a dense matrix A is about(8m/3 + 2)n3 if the computational cost for computing the norms kAk2 and kA−1 k2 is of theorder O(n2 ). Similarly, the cost of Algorithm 2 is about (8ms+1 /3 + 2(s + 1))n3 . Further, thecost incurred in computing Aα b with Algorithm 1 is mcabscissa +cmul +cnorm , where cabscissais the computational cost of computing F̃DE (x)b, cmul is that of a matrix-vector multiplication,and cnorm is that of computing the norms kAk2 and kA−1 k2 . The computational cost of
Algorithm 2 is ms+1 cabscissa + (s + 1)cmul + cnorm . Because one can use low-accuracy
norms, the total computational cost will be almost proportional to the total number of abscissa
points.
3. Convergence of the DE formula for HPD matrices. When A is an HPD matrix, theanalysis of a quadrature formula for Aα can be reduced to that of the quadrature formula forthe scalar fractional power. For the DE formula, the following relation holds:
Z ∞FDE (x) dx − h
−∞

∞
X
k=−∞

FDE (kh)2

Z ∞fDE (x, λ) dx − h= maxλ∈Λ(A)−∞∞X

fDE (kh, λ) ,k=−∞where h is the mesh size, Λ(A) is the spectrum of A, and fDE is the scalar counterpart ofFDE (x). Further, fDE can be expressed as(3.1)

fDE (x, λ) =sin(απ)λ exp(απ sinh(x)/2) cosh(x).
2exp(π sinh(x)/2) + λSection 3.1
R ∞ presents our analysis of the discretization error of the trapezoidal rule for the
integral −∞ fDE dx, i.e., the DE formula for λα . The error of the DE formula for Aα
is discussed in Section 3.2, and the DE formula is compared with Gaussian quadrature in
Section 3.3.
3.1. Discretization error of the DE formula for λα . First, we recall the following
upper bound for the discretization error of the trapezoidal rule for an integral of an analytic
function on the real axis:

ETNA
Kent State University and
Johann Radon Institute (RICAM)COMPUTING Aα WITH THE DE FORMULA

567

T HEOREM 3.1 (see, [18, Thm. 5.1]). Suppose that g is analytic in the strip
Dd = {z ∈ C : |Im(z)| < d} for some d > 0. Further, suppose that g(z) → 0 uniformly
as |z| → ∞ in the strip, and for some c, it satisfies
Z ∞
(3.2)
|g(x + iy)| dx ≤ c−∞for all y ∈ (−d, d). Then, for any h > 0, the sum h
h

∞
X
k=−∞

P∞

k=−∞ g(kh) exists and satisfies

Z ∞g(kh) −

g(x) dx ≤−∞2c.
exp(2πd/h) − 1

To apply Theorem 3.1 to fDE , we need to identify the strip in which fDE is analytic. We
identify such a region using the following proposition:
P ROPOSITION 3.2. For λ > 0, fDE (z, λ) defined in (3.1) is analytic in the strip Dd0 (λ) ,
where
s

p
2 + 5π 2 /4 −2 + 5π 2 /4]2 − π 4
(log
λ)
[(log
λ)
.
d0 (λ) = arcsin (3.3)
π 2 /2
Proof. We prove Proposition 3.2 by calculating the imaginary part of a pole of fDE closest
to the real axis among all the poles of fDE . Because exp(απ sinh(z)/2) cosh(z) is analytic
in C, we consider only the poles of 1/[exp(π sinh(z))/2) + λ].
In fact, the poles are calculated exactly. Consider the following equation:


2
(3.4)
exp
sinh z = −λ.
πLet x and y be the real and imaginary parts of z in (3.4), respectively. Then, (3.4) can be
rewritten as
π

π

π
(3.5)
exp
sinh z = exp
sinh x cos y + i cosh x sin y ,
2
2
2
(3.6)
−λ = exp(log λ + ikπ)(k = ±1, ±3, ±5 . . . ).
From the right-hand side of (3.5) and that of (3.6), it follows that
(3.7)

π sinh(x) cos(y)/2 = log λ,(3.8)

π cosh(x) sin(y)/2 = kπ.

To solve (3.7), we consider two cases: either λ is 1 or not. In case λ = 1, equation (3.7) leads
to cos y = 0. Therefore, the minimum absolute solution y is ±π/2.
In case λ 6= 1, cos y 6= 0 because log λ 6= 0. After squaring the left and right of both
equations of (3.7), we obtain 2
(log λ)2π+(3.9)
sin2 y = k 2 π 2
4
1 − sin2 yby substituting the relations cosh2 x = 1 + sinh2 x and cos2 y = 1 − sin2 y(6= 0). Equation (3.9) can be solved easily because (3.9) is a quadratic equation in sin2 y. Thus, theimaginary parts of the solutions are
s

p
2
2
2
2
2
2
2
2
2
2
4
(log λ) + π /4 + k π − [(log λ) + π /4 + k π ] − k π .
yk = arcsin 
π 2 /2

ETNA
Kent State University and
Johann Radon Institute (RICAM)

568

F. TATSUOKA, T. SOGABE, Y. MIYATAKE, T. KEMMOCHI, AND S.-L. ZHANG

d0( )

/2
/4
0
10 16 10 81

108 1016

F IG . 3.1. Distance d0 (λ) (in (3.3)) of a strip where fDE (z, λ) is analytic. Note that the horizontal axis islogarithmic.

We can verify that |y1 | = |y−1 | ≤ |yk | for all |k| > 1. Therefore, the imaginary parts of thepoles closest to the real axis are ±y1 .
Because y1 = π/2 when λ = 1, we can combine the two cases for λ. In conclusion,fDE (z, λ) is analytic in Dd0 (λ) , where d0 (λ) = y1 .
To apply Theorem 3.1 to fDE , we make sure that fDE satisfies the condition (3.2). Let δ
be a sufficiently small positive number. Then, fDE is continuous and bounded in Dd0 (λ)−δ . In
addition, by using the results in Appendix A, we can show that there exist constants cl , cr > 0
and l, r ∈ R such that for all |y| ≤ d0 (λ) − δ (< π/2),
απ cos(d0 (λ) − δ)sinh(|x|) cosh x|fDE (x + iy)| < cl exp −2
(1 − α)π cos(d0 (λ) − δ)|fDE (x + iy)| < cr exp −sinh x cosh x
2(x < l),(x > r).

Therefore, fDE satisfies condition (3.2).
To conclude, the discretization error of the DE formula for λα is estimated as(3.10)

λα − h

∞
X
k=−∞




2πd0 (λ)
fDE (kh, λ) ≤ O exp −.
h

The estimate (3.10) has a limitation in the sense that it only yields the convergence rate. The
error of the DE formula goes to 0 as λ → 0 since fDE (x, λ) → 0. However, (3.10) does notindicate this because d0 (λ) → 0 as λ → 0.
3.2. Convergence of the DE formula for Aα . First, let us consider the eigenvaluesof A that determine the convergence rate of the DE formula for Aα . Since the error of the
DE formula for λα is O(exp(−2πd0 (λ)/h)), its convergence is slow when | log(λ)| is large,i.e., λ is far away from 1. This is because d0 (λ) decreases monotonically when λ > 1 andincreases monotonically when 0 < λ < 1. In addition, the convergence rate of the DE formulafor λα is approximately equal to that of the DE formula for (1/λ)α because d0 (λ) = d0 (1/λ).
This relation d0 (λ) = d0 (1/λ) is true because λ appears only in the form of (log λ)2 ind0 (λ). These properties of d0 (λ) can be verified in Figure 3.1, which displays d0 (λ) forλ ∈ [10−16 , 1016 ]. Therefore, for sufficiently small h, the discretization error depends on the

ETNA
Kent State University and
Johann Radon Institute (RICAM)COMPUTING Aα WITH THE DE FORMULA

569

maximum and the minimum eigenvalues of A, say λmax , λmin , and the error can be written as∞
XAα − hFDE (kh)

k=−∞

=

max

λ∈{λmax ,λmin }λα − h2∞
X

fDE (kh, λ) .

k=−∞

For h that is not sufficiently small, the discretization error may depend on eigenvalues otherthan the extreme eigenvalues. This is because the error of the DE formula depends onthe coefficient multiplying the exponential factor in the discretization error as well as the
exponential factor.
Now, we assume that λmax λmin = 1 because this assumption minimizes the value
max{| log(λmax )|, | log(λmin )|},
which determines the convergence rate of the DE formula. Note that the previous condition
λmax λmin = 1 can be satisfied without loss of generality because we can computepAα = (λmax λmin )α/2 (A/ λmax λmin )α ,√where the product of the maximum and the minimum eigenvalues of A/ λmax λmin is 1.
Under this assumption, we have
∞
X



2πd0 (λmax )A −hFDE (kh) ≤ O exp −
h
k=−∞
2
!!
p
2πd0 ( κ(A))= O exp −hα(3.11)



because d0 (λmax ) = d0 (1/λmax ) = d0 (λmin ).
Finally, to compare the DE formula with Gaussian quadrature, let us rewrite the error (3.11) in terms of the number of abscissa points m. Assume that we have a finite interval[l, r] for which the truncation error is smaller than the discretization error. Then, the total error
can be rewritten asm−2XhA − (FDE (l) + FDE (r)) − h
FDE (l + kh)
2
k=1
2
p
 !!
2πd0
κ(A)
≤ O exp −
m
,
r−l
α

(3.12)

where h = (r − l)/(m − 1). The error (3.12) can be estimated by computing λα
max using the
m-point DE formula.
3.3. Comparison of the convergence speed. In this subsection, we compare the convergence speeds of three quadrature formulas. The first one, denoted by DE, is the m-point
DE formula (Algorithm 1). The second one, denoted by GJ1, represents the application of theGauss-Jacobi (GJ) quadrature to(3.13)

2 sin(απ)
A =
A
απ
α

Z 1
−1

h
i−1
du,
(1 − u)1/α−2 (1 + u)1/α I + (1 − u)1/α A

ETNAKent State University and
Johann Radon Institute (RICAM)

570

F. TATSUOKA, T. SOGABE, Y. MIYATAKE, T. KEMMOCHI, AND S.-L. ZHANG

which is obtained by applying t(u) = (1 + u)/(1 − u) to (1.1). Integral representations similarto (3.13) are considered in [6, 8]. The third one, denoted by GJ2, represents the application of
the GJ quadrature to(3.14)

Aα =2 sin(απ)
AπZ 1(1 − v)α−1 (1 + v)−α [(1 − v)I + (1 + v)A]

−1

dv,

−1

which is obtained by applying t(v) = (1 − v)α /(1 + v)α to (1.1). Integral representations
similar to (3.14) are considered in [2, 8].
In [8], under the assumption that λmax λmin = 1, the error of the GJ2 method is estimated2
to be



 
1 + [κ(A)]1/4
O exp −2 log(3.15)
m.
|1 − [κ(A)]1/4 |The error of GJ1 when α is a unit fraction is estimated [8] as! !!
p
1 + [κ(A)]α/2 + 2[κ(A)]α/2 (1 − cos(απ))
p(3.16)m.
O exp −2 log
1 + [κ(A)]α + 2[κ(A)]α/2 cos(απ)
When α is a non-unit fraction, the integrand in (3.13) is not analytic at u = ±1, and the GJ
quadrature for (3.13) may not converge exponentially.
Based on the above discussion, we can represent the error of the three quadrature formulasin the form O(exp(−φm)), where φ is a constant depending on the quadrature formulas, κ(A),and α. Figure 3.2 displays the values of φ for α = 0.1, 0.2, . . . , 0.9 and κ(A) ∈ [1, 1016 ]. For
DE, a finite interval [l, r] is obtained by using the subroutine GetInterval in Algorithm 1.
The parameter  is set to  = 2−53 κ(A)α/2 so that the relative error is smaller than 2−53 ≈1.1 × 10−16 . For GJ1, the speed φ is only considered when α is a unit fraction because of
the analyticity of the integrand. Note that the vertical axes in Figure 3.2 show the values φ.
Hence, the convergence is fast for large φ.
Figure 3.2 indicates that the convergence speed of DE is higher than that of GJ2 whenκ(A) is larger than about 104 . Conversely, the convergence speed of DE is lower than that of
GJ2 for small κ(A). When α is a unit fraction and κ(A) is large, the convergence speed of
GJ1 is higher than that of each of the other algorithms except for α = 0.5. From Figure 3.2,
we can select the fastest-converging quadrature formula. For example, when we compute Aα
where α = 0.5 and κ(A) = 1010 , DE is the fastest to converge among the three.
4. Numerical experiments. The numerical experiments were carried out by using Julia 1.5.1 on a Core-i7 (3.6 GHz) CPU with 16 GB RAM. The IEEE double-precision arithmetic
is used unless otherwise stated. For arbitrary precision arithmetic, the programs use the
BigFloat data type of Julia, which gives roughly 77 significant decimal digits. Abscissas
and weights in the GJ quadrature are computed with FastGaussQuadrature.jl.3 The
test matrices are listed in Table 4.1. The programs for these experiments are available on
Github.4
To avoid computing extremely large or extremely small values, test matrices are scaled
before the computation of Aα . Specifically, Aα = c−α (cA)α is computed with the value
2 Generally, the error of the m-point Gaussian quadrature is given by O(τ −2m ) for some constant τ . In thispaper, we use the representation O(exp(−2 log(τ )m)) because we want to compare the convergence of Gaussian
quadrature with that of the DE formula.
3 https://github.com/JuliaApproximation/FastGaussQuadrature.jl
4 https://github.com/f-ttok/article-powmde/

ETNA
Kent State University and
Johann Radon Institute (RICAM)COMPUTING Aα WITH THE DE FORMULA= 0.1 = 1/10
3.0 × 104
8.4 × 104

100
10 1
10 2
10 3

DE
GJ1
GJ2

100
100
10 1
10 2
10 3

108 1016= 0.4
1.1 × 104

DE
GJ2

100
100
10 1
10 2
10 3

108 1016= 0.7
1.1 × 104

DE
GJ2

100

108

1016

= 0.2 = 1/5
2.1 × 103

100
10 1
10 2
10 3

DE
GJ1
GJ2

100
100
10 1
10 2
10 3

108 1016= 0.5 = 1/2
7.0 × 101 9
1.9 × 10
DE
GJ1
GJ2

100
100
10 1
10 2
10 3

108 1016= 0.8
1.3 × 104

DE
GJ2

100

108

1016

571

= 0.3
1.2 × 104

100
10 1
10 2
10 3

DE
GJ2

100
100
10 1
10 2
10 3

108 1016= 0.6
1.0 × 104

DE
GJ2

100
100
10 1
10 2
10 3

108 1016= 0.9
1.8 × 104

DE
GJ2

100

108

1016

F IG . 3.2. Convergence speed of the following quadrature formulas: the DE formula (DE), the GJ quadraturefor (3.13) (GJ1), and the GJ quadrature for (3.14) (GJ2). The upper bound of the error for DE, GJ1, and GJ2 areshown in (3.12), (3.16), and (3.15), respectively. The horizontal axes are the condition number κ(A). The verticalaxes are the convergence speed, i.e., the constant φ in the error of the m-point quadrature formula exp(−φm).
Therefore, the convergence of a quadrature formula is fast if the values are large. The condition numbers are indicatedfor the points where the fastest quadrature formula is changed.

√
c = 1/ σmax σmin , where σmax and σmin are the maximum and the minimum singular valuesof A, respectively. In this section, we denote the scaled matrix cA by Ã for convenience.
When A is an HPD matrix, the scaling parameter c can be chosen specialized to GJ2 so thatthe error of GJ2 is small [3]; see, Appendix B. We denote this algorithm GJ2pre because we
can regard selecting c appropriately as a preconditioning to make good use of GJ2.
4.1. Test 1: Selected intervals in our algorithms. In this test, we verify that the subroutine GetInterval truncates the infinite interval appropriately. We consider the computationof the square root of test matrices according to the following procedure. First, we compute the reference solution A1/2 by using the Denman-Beavers (DB) iteration [10, p. 148]
with arbitrary precision. Next, we compute the finite interval in double precision. We set√
 = (cρ(A))α 10−7 , (cρ(A))α 10−14 , where c = 1/ σmax σmin , so that the relative error in

ETNA
Kent State University and
Johann Radon Institute (RICAM)

572

F. TATSUOKA, T. SOGABE, Y. MIYATAKE, T. KEMMOCHI, AND S.-L. ZHANGTABLE 4.1
Test matrices for the experiments. All matrices are real. The term SPD stands for symmetric positive definite.
For pores_1, we computed (−A)α because all the eigenvalues of pores_1 lie in the left half plane. For cell1,
TSOPF_RS_b9_c6, and circuit_3 we computed (A + 10−8 I)α , (A + 40.35I)α , and (A + 3I)α , respectively,
because they have real negative eigenvalues.

Test
1,2
3

4

5

Matrixex5 [7]pores_1 [7]nos4 [7]bcsstk04 [7]lund_b [7]SPD_well
SPD_ill
NS_well
NS_ills2rmt3m1 [7]
fv3 [7]poisson200
cell1 [7]TSOPF_RS_b9_c6 [7]circuit_3 [7]

n
27
30
100
132
147
100
100
100
100
5489
9801
40000
7055
7224
12127

Condition number
6.6 × 107
1.8 × 106
1.6 × 103
2.3 × 106
3.0 × 104
1.0 × 102
1.0 × 107
1.0 × 102
1.0 × 107
2.5 × 108
2.0 × 103
1.6 × 104
8.5 × 108
5.6 × 104
1.1 × 102

Properties
SPD
Nonsymmetric
SPD
SPD
SPD
SPD
SPD
NonsymmetricNonsymmetric
SPD
SPD
SPD
NonsymmetricNonsymmetricNonsymmetric

TABLE 4.2
Right end values r of the finite intervals [l, r] selected by Algorithm 1.

Matrixex5
ex5
pores_1
pores_1


10−7
10−14
10−7
10−14

svdvals
4.0189456993
4.5713980347
3.9825518994
4.5506094014

Arpack
4.0189454235
4.5713979514
3.9825518993
4.5506093993

Arpackmod
4.0501871836
4.5895014861
4.0149323090
4.5690896458

the 2-norm is smaller than 10−7 or 10−14 . The norms kÃk2 and kÃ−1 k2 are computed inthree ways. First, the computation is performed using the function svdvals in Julia, which
computes all the singular values of a matrix. For the second method, the computation is
performed to achieve a three-digit accuracy by using Arpack.jl,5 which computes the
eigenvalues and singular values by the Arnoldi method. For the third method, the computationis performed using Arpack.jl as in the second method, but we set  = ˜/(1+1/(1−10−3 ))for ˜ = 10−7 , 10−4 . The third method is denoted as Arpackmod. Then, we compute A1/2
by using the m-point DE formula with arbitrary precision to avoid the error caused by a matrix
inversion. For reference, we also compute A1/2 with the wider finite interval [−6, 6].
As a first result, the right end values r of the finite intervals [l, r] computed by
GetInterval are listed in Table 4.2. The data in Table 4.2 show that there is only aslight difference between the interval with accurate norms (svdvals) and that with roughnorms (Arpack). If we consider the error of the norms (Arpackmod), GetInterval
provides a slightly wider interval than that of svdvals. Hence, the error of the norms can beneglected.

5 https://github.com/JuliaLinearAlgebra/Arpack.jl

ETNA
Kent State University and
Johann Radon Institute (RICAM)COMPUTING Aα WITH THE DE FORMULA

ex5 (svdvals)573

ex5 (Arpack)ex5 (Arpackmod)

100

100

100

10 7

10 7

10 7

10 14

10 14

10 14

0

100

200

0

pores_1 (svdvals)

100

200

pores_1 (Arpack)

0

100

100

100

10 7

10 7

10 7

10 14

10 14

10 14

0

100

200
= 10 7

0100= 10 14

200

100

200

pores_1 (Arpackmod)

0

100

200

[l, r] = [ 6, 6]

F IG . 4.1. Convergence histories of the DE formula for A1/2 in Test 1. The vertical axes show the relative
error in the 2-norm, and the horizontal axis shows the number of abscissa points. The labels in each graph indicate
the test matrix and the method for computing the norms kÃk2 , kÃ−1 k2 . The notation svdvals indicates thatthe parameters are computed by using svdvals, Arpack indicates that the norms are computed approximately
by using Arpack.jl, and Arpackmod indicates that the norms are computed by using Arpack.jl, but  is set
so that the error of norms does not affect the estimate of the truncation error. The values  in the legend represent
the tolerance of the relative error, and “[l, r] = [−6, 6]” indicates that the approximation is computed with a wide
interval [−6, 6].

Next, in Figure 4.1, we present the convergence histories of the DE formula. Figure 4.1
indicates that the approximations achieved the required accuracy in all cases. Therefore, thetruncation error of Algorithm 1 is smaller than the given tolerance. As expected from theresults listed in Table 4.2, the computational results are accurate regardless of the accuracyof the parameters. In addition, the DE formula with selected finite intervals converges fasterthan the DE formula with the wide interval. In conclusion, GetInterval determines finite
intervals appropriately, and we can use the roughly computed parameters.
4.2. Test 2: Accuracy of Algorithm 2. In this test, we investigate the accuracy of
Algorithm 2 by computing Aα (α = 0.2, 0.5, 0.8). The reference solution Aα is computed byusing the DB iteration and the inverse Newton method [10, Alg. 7.14] with arbitrary precision.
The reference solution A0.8 is computed via A0.8 = (A1/5 )4 . Then, we compute Ãα by using
Algorithm 2 with arbitrary precision. In Algorithm 2, the parameter  is set so that the relative
error in the 2-norm is smaller than 10−7 and 10−14 .
The relative error and its estimate in Algorithm 2 are given in Figure 4.2. Figure 4.2
indicates that Algorithm 2 achieves the required accuracy for all cases. The behavior of theerrors shows that Algorithm 2 controls the truncation error in (2.12), and the behavior of theestimates shows that Algorithm 2 controls the discretization error in (2.15).
4.3. Test 3: Verifying the convergence speed of the DE formula. In this test, weverify the convergence speed of the DE formula for the fractional power of symmetric positive
definite (SPD) matrices. We compute Ã1/2 by using the DE formula (Algorithm 1) with

ETNA
Kent State University and
Johann Radon Institute (RICAM)

574

F. TATSUOKA, T. SOGABE, Y. MIYATAKE, T. KEMMOCHI, AND S.-L. ZHANG

ex5, = 0.2

ex5, = 0.5

ex5, = 0.8

100

100

100

10 7

10 7

10 7

10 14

10 14

10 14

0

200 400
pores_1, = 0.2

0200 400
pores_1, = 0.5

0200 400
pores_1, = 0.8

100

100

100

10 7

10 7

10 7

10 14

10 14

10 14

400
0
200
estimate ( = 10 14)
error ( = 10 14)

400
0
200
estimate ( = 10 7)
error ( = 10 7)

0

200

400

F IG . 4.2. Relative error and its estimate in Algorithm 1. The vertical axes show the relative error in the 2-norm,
and the horizontal axes show the number of abscissa points ms . The values  represent the tolerance for the relative
error.
1/2

 = 2−53 (cρ(Ã))α in double-precision arithmetic. Subsequently, we compute λ̃max , i.e., thesquare root of the maximum eigenvalue of Ã, by using the DE formula with an integral interval
identical to that of the DE formula for Ã1/2 . The convergence histories of the DE formula arepresented in Figure 4.3. In addition, we illustrate the estimates of the convergence rate (3.12)
in Figure 4.3.
Figure 4.3 shows that the convergence of the DE formula for Ã1/2 is almost equal to that
1/2of the DE formula for λ̃max . Furthermore, the convergence rates of the quadrature formulas
are approximately equal to the estimates. These results support our convergence rate analysis,
and hence, we can predict the convergence of the DE formula for Ãα from the behavior of the
DE formula for λ̃αmax .
4.4. Test 4: Comparison of the DE formula with Gauss-Jacobi quadrature. We
compare the convergence of the DE formula with that of GJ quadrature. For this test, we
generate four test matrices:
1. We generate two SPD matrices (SPD_well and SPD_ill).
1-1. We generate an orthogonal matrix Q by taking the orthogonal vector of the QR
decomposition of a random 100 × 100 matrix.
1-2. We generate a diagonal matrix D = diag(d1 , . . . , dn ) whose diagonal elements
are derived from the geometric sequence {di }i=1,...,100 , where d1 = κ−1/2 andd100 = κ1/2 for κ = 102 .
1-3. ASPD_well = QDQ>
1-4. We repeat Step 1–2 and 1–3 with the setting κ = 107 and generate ASPD_ill .

ETNA
Kent State University and
Johann Radon Institute (RICAM)COMPUTING Aα WITH THE DE FORMULA

lund_b

100
10 8
10 16

100

bcsstk04

100

10 16
200 1/2
0
max

nos4

100

10 8
0

575

10 8
100

A1/2

200

10 16

0

estimate

100200

F IG . 4.3. Convergence histories of the DE formula for the square root of the test matrices and its maximum
eigenvalues. The vertical axes show the relative error in the 2-norm, and the horizontal axes show the number ofabscissa points.

2. We generate two nonsymmetric matrices (NS_well and NS_ill).
2-1. We generate a random 100 × 100 matrix R.
2-2. We generate ANS_well = exp(cR) with c = cwell ≈ 8.633 × 10−2 so that
κ(ANS_well ) ≈ 102 . The matrix exponential is computed by using the exp
function in Julia.
2-3. We scale ANS_well so that the product of the maximum singular value and theminimum singular value becomes 1.
2-4. We repeat Step 2–2 and 2–3 and generate ANS_ill with the setting c = cill ≈3.029 × 10−1 so that κ(ANS_ill ) ≈ 107 .
Then, we compute Aα for α = 0.2, 0.5, 0.8. The reference solutions are computed as in Test 1.
The convergence histories are presented in Figure 4.4. Figure 4.4 shows that the DE formula
works well in all cases. In addition, when α is a non-unit fraction (α = 0.8) and κ(A) is large(SPD_ill and NS_ill), DE is the fastest to converge among the four or three quadratureformulas.
4.5. Test 5: Computational time for Aα b. In this test, we measure the computationaltime of quadrature-based algorithms for computing Aα b. The test matrix poisson200 is
generated by A = L200 ⊗ I200 + I200 ⊗ L200 , where L200 = tridiag(−1, 2, −1) ∈ R200×200 ,I200 is the identity matrix of size 200, and ⊗ is the Kronecker product. The vector b is set to arandom vector with norm 1. The test procedure is as follows:
SPD matrices: For SPD matrices, we consider five quadrature-based algorithms: GJ1, GJ2,GJ2pre, DE, and the algorithm based on the Cauchy integral presented in [9, Sect. 3],
which is denoted by Cauchy. The algorithm Cauchy requires the computationof the complete elliptic integral of the first kind and the Jacobi elliptic functions.
We implement them by using Elliptic.jl6 for real arguments. For complex
arguments, we use the relations in [1, Eqs. (16.21.1)–(16.21.4)].
First, we estimate the number of abscissa points m such that the absolute errorkAα b − xk2 is smaller than or equal to 10−6 . For GJ1, GJ2, and DE, we canpredict the convergence of these algorithms by the scalar convergence of the extreme
eigenvalue of Ã. Hence, we count the number of abscissa points in the scalar caseand use the same number for the matrix case. For GJ2pre, we use the upper boundfor the error, which can be computed by the extreme eigenvalues of A, derived in [3].
The extreme eigenvalues are computed with three-digit accuracy using Arpack.jl.
6 https://github.com/nolta/Elliptic.jl

ETNA
Kent State University and
Johann Radon Institute (RICAM)

576

F. TATSUOKA, T. SOGABE, Y. MIYATAKE, T. KEMMOCHI, AND S.-L. ZHANG

100

SPD_well,

= 0.2

10 8
10 16
100

100

0

100

SPD_ill,

200
= 0.2

100

0

100

NS_well,

200
= 0.2

10 16
100

10 16
100

0

100

NS_ill,

200
= 0.2

10 16
100

0

100

SPD_ill,

200
= 0.5

100

200
DE

10 16

SPD_well,

= 0.8

10 16
100

0

100

SPD_ill,

200
= 0.8

10 8
0

100

NS_well,

200
= 0.5

10 16
100

0

100

NS_well,

200
= 0.8

10 8
0

100

NS_ill,

200
= 0.5

10 8
0

100
10 8

10 8

10 8
10 16

= 0.5

10 8

10 8
10 16

SPD_well,

10 8

10 8
10 16

100

10 16
100

0

100

NS_ill,

200
= 0.8

10 8
0

GJ1

100

200
GJ2

10 16

0

100

200

GJ2preF IG . 4.4. Convergence histories of the four or three quadrature formulas: The m-point DE formula (DE), theGJ quadrature for (3.13) (GJ1), the GJ quadrature for (3.14) (GJ2), and preconditioned GJ2 for HPD matrices(GJ2pre). The vertical axes show the relative error in the Frobenius norm, and the horizontal axes show the numberof abscissa points.

For Cauchy, the number of abscissa points is estimated by applying Cauchy to a
2 × 2 matrix whose eigenvalues are the extreme eigenvalues of Ã as in [8]. When theestimated number of the abscissa points is greater than 1000, we set this number to1000. Then, we compute Aα b. The linear systems in the integrand are solved using asparse direct linear solver in SuiteSparse.jl.7

7 https://github.com/JuliaLinearAlgebra/SuiteSparse.jl

ETNA
Kent State University and
Johann Radon Institute (RICAM)COMPUTING Aα WITH THE DE FORMULA

577

Nonsymmetric matrices: For general matrices, we consider three adaptive quadrature algorithms, namely, Algorithm 2 (denoted by DE) and two algorithms based on the GJ
quadrature for (3.13) and (3.14) (denoted by GJ1 and GJ2, respectively). Algorithms
GJ1 and GJ2 compute approximations by using the (2s m0 )-point GJ quadrature fors = 0, 1, 2, . . . As in [6, Alg. 5.2], the error of the GJ quadrature is estimated by
kx2m − Aα bk2 ≤ kx2m − xm k2 ,where xm is the approximation of the m-point GJ quadrature. In contrast to thealgorithm of [6, Alg. 5.2], our algorithms do not compute the Schur decomposition
and the matrix square root because A is sparse.
As for SPD matrices, we set the tolerance so that kAα b − xk2 ≤ 10−6 . For all
algorithms, we set the initial number of abscissa points m0 to 8. The extreme singularvalues for scaling the matrices are computed via Arpack.jl, and the linear systems
are solved using SuiteSparse.jl. The algorithms are stopped when the numberof the evaluations of the integrand exceeds 1000.
Table 4.3 lists the computational times and the number of evaluations of the integrandof the quadrature-based algorithms for SPD matrices. This computational time includes
the time for computing the extreme singular values and the extreme eigenvalues for scaling
the test matrices. When α is a non-unit fraction (α = 0.8) and A is an ill-conditioned
matrix (s2rmt3m1), DE is the fastest. Otherwise, GJ1 or GJ2pre are the fastest. Although
Cauchy requires a small number of abscissa points, it is not the fastest in this test because of
the use of complex arithmetic. Table 4.4 lists the times of the algorithms for nonsymmetricmatrices. For well-conditioned matrix (circuit_3), GJ2 is the fastest, otherwise, DE isthe fastest. The algorithm DE seems to be more robust than GJ1 and GJ2 because of theefficiency of the discretization error estimation.
TABLE 4.3
Comparison of quadrature-based algorithms for SPD matrices in terms of CPU time (in seconds) and thenumber of evaluations of the integrand (in parentheses). If the number of evaluations exceeds 1000, the corresponding
results are underlined. The results of the fastest algorithm among the five algorithms are written in bold font.
Matrixs2rmt3m1
s2rmt3m1
fv3
fv3
poisson200
poisson200

α
0.2
0.8
0.2
0.8
0.2
0.8

GJ1
0.59
14.37
0.30
2.10
1.43
14.34

(39)
(1000)
(23)
(193)
(25)
(352)

GJ2
7.20
(511)
10.55
(731)
0.30
(24)
0.32
(26)
1.82
(41)
1.97
(45)

GJ2pre
6.38
(451)
6.98
(485)
0.28
(23)
0.26
(21)
1.69
(38)
1.55
(34)

DE
1.07
1.52
0.36
0.38
1.59
1.65

(73)
(95)
(30)
(29)
(33)
(33)

Cauchy
1.30
(21)
1.90
(31)
0.41
(11)
0.38
(10)
2.01
(13)
1.91
(12)

TABLE 4.4
Comparison of quadrature-based algorithms for nonsymmetric matrices in terms of CPU time (in seconds)
and the number of evaluations of the integrand (in parentheses). If the number of evaluations exceeds 1000, thecorresponding results are underlined. The results of the fastest algorithm among the three algorithms are written inbold font.

Matrix
cell1
cell1
TSOPF_RS_b9_c6
TSOPF_RS_b9_c6
circuit_3
circuit_3

α
0.2
0.8
0.2
0.8
0.2
0.8

GJ1
4.00
(248)
15.88 (1016)
0.93
(120)
7.82 (1016)
2.36
(120)
2.21
(120)

GJ2
15.86 (1016)
15.87 (1016)
1.93
(248)
1.92
(248)
0.81
(24)
0.91
(24)

3.55
3.60
0.88
0.87
1.36
1.33

DE
(225)
(225)
(113)
(113)
(57)
(57)

ETNA
Kent State University and
Johann Radon Institute (RICAM)

578

F. TATSUOKA, T. SOGABE, Y. MIYATAKE, T. KEMMOCHI, AND S.-L. ZHANG

5. Conclusion. In this work, we considered the DE formula to compute Aα . To utilize
the DE formula, we proposed a method of truncating the infinite interval based on the analysisof the truncation error specialized for Aα . Based on the truncation error analysis, we presented
two algorithms, including an adaptive quadrature algorithm. We analyzed the convergence
rate of the DE formula for HPD matrices and found that when α is a non-unit fraction andκ(A) is large, the DE formula converges faster than the GJ quadrature.
We performed five numerical tests. The results of the first two tests show that ouralgorithms achieve the required accuracy. The result of the third test shows that our estimateof the convergence rate is appropriate. The results of the fourth and fifth tests show that ouralgorithms work well, especially when the algorithms based on the GJ quadrature converge
slowly.
In the future, we plan to consider the practical performance of our algorithm with parallel
computing when applied to large practical problems.
Acknowledgments. We would like to thank Prof. Ken’ichiro Tanaka for his valuablecomments on the convergence rate analysis of the DE formula. We are grateful to the twoanonymous referees for the careful reading and the comments that substantially enhanced
the quality of the manuscript. This work was supported by JSPS KAKENHI Grant Numbers
18J22501 and 18H05392. The authors would like to thank Enago (www.enago.jp) for the
English language review.
Appendix A. On the behavior of |fDE (x + iy, λ)| as x → ±∞. In this section, we
consider the asymptotic behavior of fDE defined in (3.1). We show that there exist constants
cl , cr , l, and r such that for any y satisfying |y| < d0 (λ) − δ (< π/2),
απ cos(d0 (λ) − δ)(A.1) |fDE (x + iy)| < cl exp −sinh(|x|) cosh x(x < l),2
(1 − α)π cos(d0 (λ) − δ)(A.2) |fDE (x + iy)| < cr exp −sinh x cosh x (x > r).
2
p
From | cosh(x + iy)| = cosh2 x + cos2 y − 1, it follows that
(A.3)

| cosh(x + iy)| ≤ cosh x.

Besides, for x < 0,
 απ

 απ

(A.4)
exp
sinh(x + iy) ≤ exp −sinh(|x|) cos(d0 (λ) − δ) .
2
2
Because fDE is bounded in Dd0 (λ)−δ , there exists a constant c1 that satisfies
π

(A.5)
exp
sinh(x + iy) + λ > c1 .
2
By combining (A.3), (A.4), (A.5), we get (A.1).
Next, it is true that


 απ
 απ
(A.6)
exp
sinh(x + iy) = exp
sinh(x) cos(y) .
22In addition, there exists a constant c2 such that for all x > 0 satisfying the condition
exp(π sinh x cos(d0 (λ) − δ)/2) > λ, it holds that
π

π

π

exp
sinh(x + iy) + λ ≥ exp
sinh x cos y − λ ≥ c2 exp
sinh x cos y .
2
2
2

ETNA
Kent State University and
Johann Radon Institute (RICAM)COMPUTING Aα WITH THE DE FORMULA

579

Hence,(A.7)

| exp(απ sinh(x + iy)/2)|
1(1 − α)π≤exp −sinh(x) cos(y)| exp(π sinh(x + iy)/2) + λ|
c2
2


1
(1 − α)π cos(d0 (λ) − δ)≤exp −sinh(x) .
c2
2

By combining (A.3), (A.6), and (A.7), we get (A.2).
Appendix B. A preconditioning technique of GJ2 for HPD matrices. The study [3]considers the computation of fractional powers of a positive self-adjoint operator L by using
the GJ quadrature for the integral
(B.1)

L

−α

2 sin(απ)τ 1−α=π

Z 1(1 − t)−αα−2(1 + t)

−1




1−t
τ
I + L dt,
1+t

where α ∈ (0, 1), I is the identity operator, and τ > 0 is a parameter [3, Eq. (2)]. When
L is an HPD matrix, we can rewrite (B.1) as (3.14), which is the integral used in GJ2, by
substituting L−1 = A and τ = 1.
One of the contributions of [3] is a method for selecting τ to reduce the error. The
proposed selection is(τ − (m) (m < m̄),τ=
τ + (m) (m ≥ m̄),where m is the number of abscissa points,
1/2 1/4
αµmax
2m̄ = √ log(e κ)
κ
κ=
,
µmin
2 2



 α 24em2
exp 2W,
τ − (m) = µmin
2em
α2

2
v
!2
u
1/2
1/2
u
 αµmax log(κ) t αµmax log(κ)
τ + (m) = −++ (µmax µmin )1/2  ,8m
8m
W is the Lambert W-function, and µmax and µmin are the maximum and the minimum
eigenvalue of L, respectively. One can apply the preconditioning technique to GJ2, which
computes Aα for α ∈ (0, 1), by substituting µmax = 1/λmin , µmin = 1/λmax , where
λmax and λmin are the maximum and the minimum eigenvalue of A, and then computing
Aα = τ −α (τ A)α .
REFERENCES[1] M. A BRAMOWITZ AND I. A. S TEGUN, Handbook of Mathematical Functions with Formulas, Graphs,
and Mathematical Tables, vol. 55 of National Bureau of Standards Applied Mathematics Series, U.S.
Government Printing Office, Washington, D.C., 1964.
[2] L. ACETO AND P. N OVATI, Rational approximation to the fractional Laplacian operator in reaction-diffusion
problems, SIAM J. Sci. Comput., 39 (2017), pp. A214–A228.
[3], Rational approximations to fractional powers of self-adjoint positive operators, Numer. Math., 143(2019), pp. 1–16.

ETNA
Kent State University and
Johann Radon Institute (RICAM)

580

F. TATSUOKA, T. SOGABE, Y. MIYATAKE, T. KEMMOCHI, AND S.-L. ZHANG[4] S. AOKI , M. F UKUGITA , S. H ASHIMOTO , K. I. I SHIKAWA , N. I SHIZUKA , Y. I WASAKI , K. K ANAYA ,T. K ANEKO , Y. K URAMASHI , M. O KAWA , N. T SUTSUI , A. U KAWA , N. YAMADA , AND T. YOSHI Ã L’,
An exact algorithm for any-flavor lattice QCD with KogutâĂŞSusskind fermion, Comput. Phys. Commun.,155 (2003), pp. 183–208.
[5] K. B URRAGE , N. H ALE , AND D. K AY, An efficient implicit FEM scheme for fractional-in-space reactiondiffusion equations, SIAM J. Sci. Comput., 34 (2012), pp. A2145–A2172.
[6] J. R. C ARDOSO, Computation of the matrix pth root and its Fréchet derivative by integrals, Electron. Trans.
Numer. Anal., 39 (2012), pp. 414–436.
http://etna.ricam.oeaw.ac.at/vol.39.2012/pp414-436.dir/pp414-436.pdf[7] T. A. DAVIS AND Y. H U, The University of Florida sparse matrix collection, ACM Trans. Math. Software, 38(2011), Art. 1, 25 pages.
[8] M. FASI AND B. I ANNAZZO, Computing the weighted geometric mean of two large-scale matrices and its
inverse times a vector, SIAM J. Matrix Anal. Appl., 39 (2018), pp. 178–203.
[9] N. H ALE , N. J. H IGHAM , AND L. N. T REFETHEN, Computing Aα , log(A), and related matrix functions bycontour integrals, SIAM J. Numer. Anal., 46 (2008), pp. 2505–2523.
[10] N. J. H IGHAM, Functions of Matrices, SIAM, Philadelphia, 2008.
[11] N. J. H IGHAM AND L. L IN, A Schur-Padé algorithm for fractional powers of a matrix, SIAM J. Matrix Anal.
Appl., 32 (2011), pp. 1056–1078.
, An improved Schur-Padé algorithm for fractional powers of a matrix and their Fréchet derivatives,[12]
SIAM J. Matrix Anal. Appl., 34 (2013), pp. 1341–1360.
[13] B. I ANNAZZO AND C. M ANASSE, A Schur logarithmic algorithm for fractional powers of matrices, SIAM J.
Matrix Anal. Appl., 34 (2013), pp. 794–813.
[14] M. M ORI, Discovery of the double exponential transformation and its developments, Publ. Res. Inst. Math.
Sci., 41 (2005), pp. 897–935.
[15] M. S UGIHARA, Optimality of the double exponential formula—functional analysis approach, Numer. Math.,
75 (1997), pp. 379–395.
[16] B. J. S ZEKERES AND F. I ZSÁK, Finite difference approximation of space-fractional diffusion problems: thematrix transformation method, Comput. Math. Appl., 73 (2017), pp. 261–269.
[17] H. TAKAHASI AND M. M ORI, Double exponential formulas for numerical integration, Publ. Res. Inst. Math.
Sci., 9 (1973/74), pp. 721–741.
[18] L. N. T REFETHEN AND J. A. C. W EIDEMAN, The exponentially convergent trapezoidal rule, SIAM Rev., 56(2014), pp. 385–458.