Algorithms for the computation of the matrix logarithm based on the double exponential formula

TL;DR

This paper introduces a novel approach for computing the matrix logarithm using the double exponential numerical quadrature, improving efficiency and accuracy especially for ill-conditioned matrices.

Contribution

It develops a method to select finite integration intervals and presents two algorithms, including one with error control, for matrix logarithm computation using the double exponential formula.

Findings

Enhanced convergence for ill-conditioned matrices

Effective finite interval selection method

Algorithms with error control for improved accuracy

Abstract

We consider the computation of the matrix logarithm by using numerical quadrature. The efficiency of numerical quadrature depends on the integrand and the choice of quadrature formula. The Gauss--Legendre quadrature has been conventionally employed; however, the convergence could be slow for ill-conditioned matrices. This effect may stem from the rapid change of the integrand values. To avoid such situations, we focus on the double exponential formula, which has been developed to address integrands with endpoint singularity. In order to utilize the double exponential formula, we must determine a suitable finite integration interval, which provides the required accuracy and efficiency. In this paper, we present a method for selecting a suitable finite interval based on an error analysis as well as two algorithms, and one of these algorithms addresses error control.

Figures

Tables

Table 1. Table 1: Test matrices. The condition number of A 𝐴 A is denoted by κ 2 ​ ( A ) = ‖ A ‖ 2 ​ ‖ A − 1 ‖ 2 subscript 𝜅 2 𝐴 subscript norm 𝐴 2 subscript norm superscript 𝐴 1 2 \kappa_{2}(A)=\|A\|_{2}\|A^{-1}\|_{2} .

Matrix	Size	${\kappa }_2\left(A\right)$	Structure
SPD 1	50	$1.0\times {10}^1$	SPD
SPD 2	50	$1.0\times {10}^4$	SPD
SPD 3	50	$1.0\times {10}^7$	SPD
parter [17]	10	$2.4\times {10}^0$	Nonsymmetric
frank [17]	10	$2.9\times {10}^7$	Nonsymmetric
vand [17]	10	$3.1\times {10}^{12}$	Nonsymmetric
bcsstk02 [4]	66	$4.3\times {10}^3$	SPD
bcsstk03 [4]	112	$6.8\times {10}^6$	SPD
ck104 [4]	104	$5.5\times {10}^3$	Nonsymmetric

Table 2. Table 2: Comparison between Algorithm 2 (based on the DE formula) and Algorithm 3 (based on the GL quadrature), in terms of the number of integrand evaluations (in bold) and the relative error when the algorithm stopped (in parentheses). The notation “ – (–)” means that the algorithms did not stop before the number of integrand evaluations reached the limit.

	Algorithm 2 (DE)		Algorithm 3 (GL)
$\zeta$	${10}^{-8}$	${10}^{-11}$	${10}^{-8}$	${10}^{-11}$
SPD 1	61 ($2.2\times {10}^{-9}$)	61 ($2.7\times {10}^{-12}$)	48 ($4.6\times {10}^{-16}$)	112 ($5.7\times {10}^{-16}$)
SPD 2	121 ($6.7\times {10}^{-10}$)	241 ($6.4\times {10}^{-13}$)	1008 ($1.8\times {10}^{-15}$)	1008 ($1.8\times {10}^{-15}$)
SPD 3	241 ($3.0\times {10}^{-10}$)	481 ($4.9\times {10}^{-13}$)	– (–)	– (–)
parter	61 ($2.6\times {10}^{-9}$)	121 ($2.3\times {10}^{-12}$)	112 ($3.3\times {10}^{-16}$)	112 ($3.3\times {10}^{-16}$)
frank	481 ($1.0\times {10}^{-12}$)	1921 ($2.1\times {10}^{-13}$)	496 ($1.5\times {10}^{-11}$)	– (–)
vand	– (–)	– (–)	– (–)	– (–)
bcsstk02	121 ($2.8\times {10}^{-9}$)	121 ($3.1\times {10}^{-12}$)	496 ($1.7\times {10}^{-15}$)	1008 ($1.0\times {10}^{-15}$)
bcsstk03	241 ($1.4\times {10}^{-9}$)	241 ($1.5\times {10}^{-12}$)	– (–)	– (–)
ck104	121 ($6.7\times {10}^{-10}$)	121 ($7.4\times {10}^{-13}$)	496 ($2.2\times {10}^{-15}$)	496 ($2.2\times {10}^{-15}$)

Equations

$$\exp (X)=A,$$

$$\displaystyle\log(A)=(A-I)\int_{0}^{1}\bigl{[}t(A-I)+I\bigr{]}^{-1}\,\mathrm{d}t.$$

$$\log (A)=(A-I){\int }_{-1}^1{\left[(1+u)(A-I)+2I\right]}^{-1}\mathrm{d}u.$$

$$\log (A)=(A-I){\int }_{-\infty }^{\infty }{F}_{\mathrm{DE}}(x)\mathrm{d}x,$$

$$F_{\mathrm{DE}}(x):=\cosh (x){\mathrm{sech}}^2(\sinh (x)){\left[(1+\tanh (\sinh (x)))(A-I)+2I\right]}^{-1}.$$

$$\Vert \log (A)-(A-I){\int }_l^r{F}_{\mathrm{DE}}(x)\mathrm{d}x\Vert .$$

$$\Vert \log (A)-(A-I){\int }_l^r{F}_{\mathrm{DE}}(x)\mathrm{d}x\Vert =\Vert (A-I)\left[{\int }_{-\infty }^l{F}_{\mathrm{DE}}(x)\mathrm{d}x+{\int }_r^{\infty }{F}_{\mathrm{DE}}(x)\mathrm{d}x\right]\Vert .$$

$$\Vert (A-I)\left[{\int }_{-\infty }^l{F}_{\mathrm{DE}}(x)\mathrm{d}x+{\int }_r^{\infty }{F}_{\mathrm{DE}}(x)\mathrm{d}x\right]\Vert$$

$$\le \Vert (A-I){\int }_{-\infty }^l{F}_{\mathrm{DE}}(x)\mathrm{d}x\Vert +\Vert (A-I){\int }_r^{\infty }{F}_{\mathrm{DE}}(x)\mathrm{d}x\Vert .$$

$$\Vert (A-I){\int }_{-\infty }^l{F}_{\mathrm{DE}}(x)\mathrm{d}x\Vert =\Vert (A-I){\int }_0^{a}F(t)\mathrm{d}t\Vert ,$$

$$F(t):=[t(A-I)+I{]}^{-1},a:=\frac{\tanh (\sinh (l))+1}{2}.$$

$$\Vert (A-I){\int }_0^{a}F(t)\mathrm{d}t\Vert \le \frac{3a}{2}\Vert A-I\Vert .$$

$$\Vert t(I-A)\Vert \le \frac{1}{2}(<1).$$

$$F(t)=[t(A-I)+I{]}^{-1}=\sum _{k=0}^{\infty }{t}^k(I-A{)}^k.$$

$$(A-I){\int }_0^{a}F(t)\mathrm{d}t$$

$$=(A-I)\sum _{k=0}^{\infty }\frac{a^{k+1}}{k+1}(I-A{)}^k$$

$$=a(A-I)+a(A-I)\sum _{k=1}^{\infty }\left[\frac{a^k}{k+1}(I-A{)}^k\right].$$

$$\Vert (A-I){\int }_0^{a}F(t)\mathrm{d}t\Vert$$

$$\le a\Vert A-I\Vert +a\Vert A-I\Vert \sum _{k=1}^{\infty }\left[\frac{1}{2}{\left(\frac{1}{2}\right)}^k\right]$$

$$=\frac{3a}{2}\Vert A-I\Vert .$$

$$\Vert (A-I){\int }_r^{\infty }{F}_{\mathrm{DE}}(x)\mathrm{d}x\Vert =\Vert (A-I){\int }_b^{1}F(t)\mathrm{d}t\Vert ,$$

$$\Vert (A-I){\int }_b^{1}F(t)\mathrm{d}t\Vert \le \left(-\log (b)+\frac{1-b}{2b}\right)\Vert A-I\Vert \Vert A^{-1}\Vert .$$

$$\Vert \frac{t-1}{t}{A}^{-1}\Vert \le \frac{1}{2}(<1).$$

$$F(t)$$

$$=\frac{1}{t}{A}^{-1}\sum _{k=0}^{\infty }{\left(\frac{t-1}{t}\right)}^k{A}^{-k}.$$

$$(A-I){\int }_b^{1}F(t)\mathrm{d}t$$

$$=(A-I)A^{-1}\left[-\log (b)I+\sum _{k=1}^{\infty }{A}^{-k}{\int }_b^1\frac{1}{t}{\left(\frac{t-1}{t}\right)}^k\mathrm{d}t\right].$$

$$\Vert (A-I){\int }_b^{1}F(t)\mathrm{d}t\Vert$$

$$\sum _{k=1}^{\infty }\left[{\Vert A^{-1}\Vert }^k{\int }_b^1\frac{1}{t}{\left(\frac{1-t}{t}\right)}^k\mathrm{d}t\right]$$

$$=\frac{1-b}{b}\sum _{k=1}^{\infty }\left[\frac{1}{k+1}{\Vert \frac{1-b}{b}{A}^{-1}\Vert }^k\right]$$

$$\le \frac{1-b}{b}\sum _{k=1}^{\infty }\left[\frac{1}{2}{\left(\frac{1}{2}\right)}^k\right]$$

$$=\frac{1-b}{2b}.$$

$$\Vert \log (A)-(A-I){\int }_l^r{F}_{\mathrm{DE}}(x)\mathrm{d}x\Vert \le \frac{3}{2}a\Vert A-I\Vert +\left(-\log (b)+\frac{1-b}{2b}\right)\Vert A-I\Vert \Vert A^{-1}\Vert .$$

$$\frac{\Vert \log (A)-(A-I){\int }_l^r{F}_{\mathrm{DE}}(x)\mathrm{d}x\Vert }{\Vert \log (A)\Vert }\le ϵ.$$

$$\Vert \log (A)\Vert \ge |\log (\rho (A))|,$$

Full text

1 Introduction

A logarithm of $A\in\mathbb{R}^{n\times n}$ is defined as any matrix $X$ such that

[TABLE]

where $\exp(X):=I+X+\frac{1}{2!}X^{2}+\frac{1}{3!}X^{3}+\cdots$ [10, p. 269]. If all eigenvalues of $A$ lie in the set $\mathbb{C}\setminus(-\infty,0]$, there is a unique logarithm of $A$ whose eigenvalues all lie in the strip $\{z\in\mathbb{C}:-\pi<\mathrm{Im}(z)<\pi\}$ [10, Thm. 1.31]. This logarithm is called the principal logarithm of $A$, and denoted by $\log(A)$. Throughout this paper, we assume that all eigenvalues of $A$ lie in the set $\mathbb{C}\setminus(-\infty,0]$, and we consider the principal logarithm of $A$.

The matrix logarithm is utilized in many fields of research, such as quantum mechanics [16], quantum chemistry [9], biomolecular dynamics [11], buckling simulation [14], and machine learning [8, 12, 6, 7]. The computational methods include the inverse scaling and squaring (ISS) algorithm [1], an algorithm based on the arithmetic-geometric mean (AGM) iteration [3], and numerical quadrature. In this paper, we focus on numerical quadrature, which employs the following integral representation (see e.g. [10, Thm. 11.1]):

[TABLE]

First, numerical quadrature can make use of the sparseness of $A$ if $A$ is sparse. It is useful when computing the multiplication of the matrix logarithm and a vector $\log(A)\bm{b}~{}(\bm{b}\in\mathbb{R}^{n})$, which appears in applications such as [8, 6, 7]111 In the three applications, the computation of $\log(\det A)$ is performed based on some connections between $\log(A)\bm{b}$ and $\log(\det A)$. For more details, see, e.g., [13]. , without computing and storing dense matrices. Conversely, the ISS algorithm and the algorithm based on the AGM iteration include the computation of the matrix square root, which means that the calculation involves dense matrices even if $A$ is sparse. The second reason is that numerical quadrature is potentially more favorable for parallel computers because of independent computation of the integrand on each abscissa.

Because the integrand in (2) includes matrix inversion, the computational cost of numerical quadrature depends on the number of evaluations of the integrand. Although numerical quadrature is suitable for parallelization, the quadrature formula should be selected carefully to reduce the computational cost and save computational resources.

The method conventionally used to compute (2) is the Gauss–Legendre (GL) quadrature. If the spectral radius of $A-I$ is smaller than 1, the GL quadrature, which can be regarded as a rational approximation of $\log(A)$, coincides with the Padé approximants of $\log(A)$ at $I$ [5, Thm. 4.3]. Therefore, it is natural to use the GL quadrature to reduce the number of abscissas when $A$ is close to $I$. However, the convergence of the GL quadrature becomes slow when $A$ is not close to $I$. For example, the convergence in our experiments became slow when $A$ was ill-conditioned which may be explained by rapid changes in the integrand value when it is closer to the endpoint of the interval.

In this paper, we consider the double exponential (DE) formula [15], which can be used to compute integrals with singularities at one (or both) of the endpoints. For this reason, the DE formula may be useful in scenarios in which the GL quadrature does not perform well. However, when using the DE formula, a finite interval needs to be selected because the integrand in (2) is transformed into a corresponding function on the infinite interval. This selection is important, i.e., if the finite interval is too narrow, the accuracy of the computational result becomes low, but if it is too wide, the convergence of the DE formula becomes slow.

By performing an error analysis, we provide a method of selecting the appropriate finite interval, as well as two algorithms for the computation of $\log(A)$ based on the $m$-point DE formula.

The remainder of this paper is organized as follows: in Section 2, we present an error analysis and propose two algorithms; in Section 3, we show the results of numerical experiments; in Section 4, we conclude the study.

Notation: Unless otherwise stated, $\|\cdot\|$ denotes a consistent matrix norm, e.g., the $p$-norm or the Frobenius norm, and $\|\cdot\|_{2}$ and $\|\cdot\|_{F}$ denote the 2-norm and the Frobenius norm, respectively. The inverse functions of $\sinh$ and $\tanh$ are referred to as $\mathrm{asinh}$ and $\mathrm{atanh}$, respectively.

2 Algorithms for the computation of $\log(A)$ based on the DE formula

In this section we propose a method of selecting a finite interval for the DE formula by estimating the interval truncation error and present two algorithms. Before considering the truncation error, let us apply the following transformations to (2). By substituting $u=2t-1$ in (2), we obtain

[TABLE]

Then, by applying the DE transformation $u=\tanh(\sinh(x))$, it follows that

[TABLE]

where

[TABLE]

The matrix $(1+\tanh(\sinh(x)))(A-I)+2I$ in the integrand $F_{\mathrm{DE}}$ is nonsingular for any $x\in$ $\mathbb{R}$. In Subsection 2.1, we derive an upper bound on the error between the integral in (4) and the same integral defined in the finite interval $[l,r]$,

[TABLE]

In Subsection 2.2, we propose a method of selecting the interval $[l,r]$ so that the relative truncation error is small than or approximately equal to a given tolerance $\epsilon$. Our algorithms are described in Subsection 2.3.

2.1 Estimation of the error from the interval truncation

The error that stems from the interval truncation (6) can be rewritten as

[TABLE]

Using the triangle inequality in the right-hand side (RHS) of (7), it holds that

[TABLE]

By estimating the RHS of (8), we obtain an upper bound on (6).

Initially, we focus on the first term which is on the RHS of (8). To avoid cumbersome notation, instead of $F_{\mathrm{DE}}$, which including hyperbolic functions, we consider the integrand of (2). By applying the transformation $t=[\tanh(\sinh(x))+1]/2$, we have

[TABLE]

where,

[TABLE]

The following lemma shows an upper bound on the RHS of (10), if $a$ is small enough to warrant the use of the Neumann series expansion of $F(t)$.

Lemma 2.1**.**

Suppose that $A\neq I$. Then, for a\in\bigl{(}0,1/(2\|A-I\|)\bigr{]}, we have

[TABLE]

Proof**.**

For all $t\in[0,a]$ where a\in\bigl{(}0,1/(2\|A-I\|)\bigr{]}, it holds that

[TABLE]

By applying the Neumann series expansion to $F(t)$ we get the following:

[TABLE]

Therefore, the integral of $F$ can be rewritten as follows:

[TABLE]

By using triangle inequality and consistency of the norm we get the following:

[TABLE]

$\square$

The calculation to estimate the second term on the RHS of (8) is similar to that of the first term. By applying the transformation $t=[\tanh(\sinh(x))+1]/2$, we get the following:

[TABLE]

where $b=[\tanh(\sinh(r))+1]/2$.

The following lemma shows an upper bound on (21), if $b$ is close enough to 1 to warrant the use of the Neumann series expansion of $F(t)$.

Lemma 2.2**.**

For b\in\bigl{[}2\|A^{-1}\|/(2\|A^{-1}\|+1),1\bigr{)}, we have

[TABLE]

Proof**.**

The outline of this proof is similar to that of Lemma 2.1. For all $t\in[b,1]$ where b\in\bigl{[}2\|A^{-1}\|/(2\|A^{-1}\|+1),1\bigr{)}, it holds that

[TABLE]

By applying the Neumann series expansion to $F(t)$, we get

[TABLE]

Therefore, the integral of $F$ can be rewritten as follows:

[TABLE]

By using the triangle inequality and consistency of the norm, it follows

[TABLE]

In (28), it holds that $|-\log(b)|=-\log(b)$ because b\in\bigl{[}2\|A^{-1}\|/(2\|A^{-1}\|+1),1\bigr{)}, and $|(t-1)/t|=(1-t)/t$ because $t\in[b,1]$. For the second term in the bracket on the RHS of (28), we have:

[TABLE]

By substituting (32) in (28), we obtain the inequality (22). $\square$

In the final part of this Subsection, we estimate the upper bound on (6).

Proposition 2.1**.**

Suppose that $A\neq I$. For a given interval $[l,r]$, let $a=[\tanh(\sinh(l))+1]/2$ and $b=[\tanh(\sinh(r))+1]/2$. Then, if $a\leq 1/(2\|A-I\|)$ and $b\geq 2\|A^{-1}\|/(2\|A^{-1}\|+1)$, it holds that

[TABLE]

Proof**.**

By combining (7) and (8), and as well as substituting (12) and (22) in (8), we get (33). $\square$

2.2 Setting the integration interval

To develop algorithms for computing $\log(A)$ based on the DE formula, we need to determine the appropriate finite integration interval, $[l,r]$ in advance. The finite interval should be ideally set so that the relative error is guaranteed to be smaller than or equal to a given tolerance, $\epsilon>0$, i.e.,

[TABLE]

To accomplish this, a lower bound on $\|\log(A)\|$ must be estimated. The following lemma shows a lower bound in terms of the spectral radius of $A$.

Lemma 2.3**.**

Let $\rho(\cdot)$ be the spectral radius, i.e., the largest absolute value of eigenvalues. Then, the following two inequalities hold:

[TABLE]

Proof**.**

By using the consistency of the norm and the fact any eigenvalue of $\log(A)$ is equal to $\log(\lambda)$, for some $\lambda$ being an eigenvalue of $A$, we have the lower bound in (35) as $\|\log(A)\|\geq\rho(\log(A))\geq|\log(\rho(A))|$. The lower bound in (36) can be obtained in a similar way. $\square$

In the following proposition, we show how to set a finite interval such that the relative truncation error in 2-norm is smaller than or equal to a given tolerance, $\epsilon>0$.

Proposition 2.2**.**

Suppose that $A\neq I$. Let $\theta$ be a lower bound on $\|\log(A)\|_{2}$, the tolerance $\epsilon>0$ satisfy

[TABLE]

and $s$ be a real positive solution of the equation

[TABLE]

If

[TABLE]

where

[TABLE]

then, it holds that

[TABLE]

Proof**.**

From the definition of $a$ and $b$, it is true that $a\leq 1/(2\|A-I\|_{2})$ and $b\geq 2\|A^{-1}\|_{2}/(2\|A^{-1}\|_{2}+1)$. In addition, from

[TABLE]

it follows that $a<b$, and $l<r$. Therefore, we can choose a finite interval $[l,r]$ that satisfies the assumptions of Proposition 2.1. By dividing the inequality (33) by $\|\log(A)\|_{2}\geq\theta$, it follows:

[TABLE]

From the definition of $a$ and $b$, it holds that $a\leq\theta\epsilon/(3\|A-I\|_{2})$ and $b\geq s$. Therefore,

[TABLE]

We obtain (41) by substituting (46) in (45). $\square$

2.3 Algorithms

In this subsection we establish two algorithms based on the results in subsections 2.1 and 2.2. One of the algorithms is designed to compute $\log(A)$ using the $m$-point DE formula on a finite interval with an interval truncation error smaller than or approximately equal to a given tolerance, $\epsilon>0$. The other algorithm is an adaptive quadrature algorithm designed to compute $\log(A)$ by automatically adding abscissas until the error is smaller than or approximately equal to a given tolerance $\zeta>0$.

If the tolerance $\epsilon$ given in Proposition 2.2 is sufficiently small, a linear approximation to the nonlinear equation (38) can be used to determine an appropriate interval. We describe our calculation in detail below.

Suppose that $\epsilon$ is sufficiently small and the solution $s$ of (38) is approximately equal to 1. Then, because $3(1-s)/2$ is the first-order Taylor approximation to $-\log(s)+(1-s)/2s$ at $s=1$, the solution $s$ can be approximated by using the solution $\tilde{s}$ of the following equation:

[TABLE]

The solution of (47) is given by

[TABLE]

Under the assumptions of Proposition 2.2, by choosing $\tilde{b}$ as

[TABLE]

instead of (40), and setting $\tilde{r}=\mathrm{asinh}(\mathrm{atanh}(2\tilde{b}-1))$, the interval truncation will be smaller than or approximately equal to $\epsilon$. The summary of the first algorithm, which computes $\log(A)$ using the $m$-point DE formula whose interval truncation error is smaller than or approximately equal to $\epsilon$ is as shown in Algorithm 1.

When $\epsilon$ is sufficiently small, an accurate computation of $\|I-A\|_{2}$, $\|A^{-1}\|_{2}$ and $\rho(A)$ at Step 2 of Algorithm 1 may not be required because the errors that stem from of these values have little effect on the accuracy of $\log(A)$. We give more detail in the following paragraph.

Assume that $\epsilon$ in Algorithm 1 is sufficiently small, and that $a$ and $b$ in Step 9 are chosen as

[TABLE]

where $\theta$ is a lower bound on $\|\log(A)\|_{2}$. Let $\Delta_{1}$ and $\Delta_{2}$ be the relative errors of $\theta/\|A-I\|_{2}$ and $1/\|A^{-1}\|_{2}$ respectively. Then, the computational result of $a$ is equal to

[TABLE]

and that of $b$ is equal to

[TABLE]

Therefore, the upper bound on the truncation error $\epsilon$ computed by considering the relative errors $\Delta_{1},~{}\Delta_{2}$ is almost equal to the upper bound on the truncation error when the tolerance is set as $\epsilon(1+\Delta_{1}+\Delta_{2}+\Delta_{1}\Delta_{2})$. For example, when $\Delta_{1},\Delta_{2}=10^{-2}$, $\epsilon(1+\Delta_{1}+\Delta_{2}+\Delta_{1}\Delta_{2})\approx 1.02\epsilon$, which means that the upper bound on the truncation error changes by approximately 2%. If $\epsilon$ is sufficiently small, the effect of these errors will be negligible.

At Step 3, a lower bound on $\|\log(A)\|_{2}$ is computed based on (35). By setting $\theta=\max\{|\log(\rho(A))|,$ $|\log(\rho(A^{-1}))|\}$, a tighter lower bound can be obtained. In particular, when $A$ is positive definite, $|\log(\rho(A^{-1}))|$ can be obtained without additional computation because $\rho(A^{-1})=\|A^{-1}\|_{2}$ is already computed in Step 2.

The computational cost of Algorithm 1 for dense $A$ is $(2m+2)n^{3}+\mathcal{O}(n^{2})$. When $A$ is sparse, evaluating $\log(A)\bm{b}$ using Algorithm 1 has computational cost $mc_{\mathrm{abscissa}}+c_{\mathrm{mul}}+c_{\mathrm{param}}$, where $c_{\mathrm{abscissa}}$ is the computational cost of computing $F_{\mathrm{DE}}(x)\bm{b}$, $c_{\mathrm{mul}}$ is the computational cost of a matrix-vector multiplication, and $c_{\mathrm{param}}$ is the computational cost of computing parameters $\rho(A)$, $\|A-I\|_{2}$, and $\|A^{-1}\|_{2}$. If the parameters are computed approximately and $F_{\mathrm{DE}}(x)\bm{b}$ is computed accurately, then $c_{\mathrm{param}}$ will be smaller than $c_{\mathrm{abscissa}}$. Therefore, the computational cost will largely depend on $m$.

Once an appropriate finite interval is obtained, the accuracy of the DE formula can be improved with the following procedure. Let $m$ be the number of abscissas, $h=(r-l)/(m-1)$ be the mesh size, and $T(h)$ be the trapezoidal rule for the mesh size $h$:

[TABLE]

Then, $T(h/2)$ can be computed using $T(h)$:

[TABLE]

In addition, we can apply the following estimation of the trapezoidal error for a sufficiently small value of $h$ using [2, Eq. (4.3)]:

[TABLE]

Our adaptive quadrature algorithm which is based on (55) is presented as Algorithm 2.

The computational cost of Algorithm 2 for dense $A$ is $(2m_{k+1}+2)n^{3}+\mathcal{O}(n^{2})=[2^{k+1}(m_{0}-1)+4]n^{3}+\mathcal{O}(n^{2})$, and the computational cost of $\log(A)\bm{b}$ with sparse $A$ is $[2^{k}(m_{0}-1)+1]c_{\mathrm{abscissa}}+c_{\mathrm{mul}}+c_{\mathrm{param}}$.

3 Numerical experiments

The numerical experiments were carried out using Julia 1.0 on a Core-i7 (3.4GHz) CPU with 16GB RAM. We used the IEEE double precision arithmetic. We computed abscissas and weights in the GL quadrature with QuadGK.jl (https://github.com/JuliaMath/QuadGK.jl).

3.1 Experiment 1: Checking the convergence

In this experiment, we checked the convergence of the GL quadrature and the DE formula. Our test matrices are presented in Table 1.

We generated the first three matrices in Table 1 by using the following procedure:

We generated an orthogonal matrix $Q$ by QR decomposition of a random $50\times 50$ matrix. 2. 2.

We generated a diagonal matrix whose diagonal elements were from the geometric sequence: $\{d_{i}\}_{i=1,\ldots,50}$ where $d_{1}=\kappa^{-1/2}$ and $d_{50}=\kappa^{1/2}$ for $\kappa=10^{1}$. 3. 3.

$A=QDQ^{\top}$. 4. 4.

We repeated Step 2 and Step 3 by setting $\kappa$ as $10^{4}$ and $10^{7}$.

The experimental procedure is as follows:

We scaled the test matrices as $\tilde{A}=(10/\rho(A))A$ because some matrices had values that were too large to use in computation. 2. 2.

We computed the reference $\log(\tilde{A})$ with the arbitrary precision arithmetic and rounded the result to double precision. We implemented the ISS algorithm [10, Alg. 11.10] with the BigFloat type of Julia. 3. 3.

We computed $\log(\tilde{A})$ using Algorithm 1, where $\epsilon=2^{-53}\approx 1.1\times 10^{-16}$. If the test matrix was symmetric positive definite, we set $\theta=\max\{|\log(\rho(\tilde{A}))|,|\log(\rho(\tilde{A}^{-1}))|\}$ as stated in Subsection 2.3. We computed $\rho(\tilde{A})$ using the eigvals function of Julia, which computes all eigenvalues of $\tilde{A}$. Similarly, we computed $\|I-\tilde{A}\|_{2}$ and $\|\tilde{A}^{-1}\|_{2}$ using the svdvals function, which computes all singular values of $\tilde{A}$ 222 If $A$ is large, using eigvals and svdvals may be inefficient. Instead, for Julia, a package Arapack.jl (https://github.com/JuliaLinearAlgebra/Arpack.jl) is available, which can compute a part of eigenvalues and singular values based on the Lanczos (or the Arnoldi) process with the desired accuracy. We present some numerical results, for which $\rho(\tilde{A})$, $\|\tilde{A}-I\|_{2}$, and $\|\tilde{A}^{-1}\|_{2}$ are computed with low accuracy, as shown in Appendix A. . 4. 4.

We computed $\log(\tilde{A})$ by applying the GL quadrature to (3).

Figure 1 shows the convergence histories of the DE formula and the GL quadrature for each matrix. Several observations can be made:

• •

The accuracy of the DE formula is almost the same as that of the GL quadrature, and the accuracies of the DE formula and the GL quadrature depend on the condition number of test matrices.

• •

For well-conditioned matrices, such as SPD 1 and parter matrix, the GL quadrature converged faster than the DE formula. Conversely, for the ill-conditioned matrices, such as SPD 3 and vand matrix, the DE formula converged faster than the GL quadrature.

The above observations suggest that Algorithm 1 selects an appropriate interval and the DE formula is suitable for ill-conditioned matrices.

3.2 Experiment 2: Checking Algorithm 2

In this experiment, we check the performance of Algorithm 2 by using the same matrices that were used in Experiment 1 (see Subsection 3.1). We compared Algorithm 2 with Algorithm 3, which is based on the GL quadrature (see Appendix B).

We conducted the experiment using the following procedure:

We computed $\log(\tilde{A})$ by using Algorithm 2. We set $m_{0}=16$ and $\zeta=\epsilon\in\{10^{-8},10^{-11}\}$. In Step 2 of Algorithm 2, which calls Algorithm 1, we computed the spectral radius and the 2-norm of matrices using eigvals and svdvals functions, as is done in Experiment 1. We stopped the computation once the number of integrand evaluations reached 1921. 2. 2.

We computed $\log(\tilde{A})$ by using Algorithm 3. We set $m_{0}=16$ and $\zeta\in\{10^{-8},10^{-11}\}$. In Step 2 of Algorithm 3, the lower bound $\theta$ was computed using (35) and (36) in the same way as was done for the DE formula. If the number of integrand evaluations was more than 2032, we stopped the computation.

The number of integrand evaluations and the corresponding relative error when the two algorithms stopped are shown in Table 2.

Several observations can be made:

• •

Algorithm 2 successfully computed the logarithm with the desired accuracy within 2000 integrand evaluations for all test matrices, except vand matrix, whereas Algorithm 3 could not succeed for three matrices. For vand matrix, the stopping criterion $\zeta$ may be too strict because the convergence history of vand matrix (as shown in Figure 1) hardly reach the value of $10^{-8}$. Our future studies will focus on the method for selecting a suitable stopping criterion.

• •

Even if the condition number of $A$ is small as is the case for SPD 1 and parter matrix, the number of integrand evaluations of Algorithm 2 could be smaller than that of Algorithm 3 because Algorithm 2 can reuse all previous results when improving accuracy.

These observations show that Algorithm 2 can be a practical choice for the computation of the matrix logarithm by numerical quadrature.

4 Conclusion

In this paper, we focused on the DE formula as a new choice for the numerical quadrature formula of $\log(A)$. In order to utilize the DE formula, we proposed a method for selecting an appropriate finite interval based on error analysis, and we proposed two algorithms for practical computation.

We carried out two numerical experiments. The first experiment showed that the DE formula converged faster than the GL quadrature for ill-conditioned matrices. The second experiment demonstrated that the proposed adaptive quadrature algorithm worked well with appropriate stopping criteria.

Our future work will focus on three problems. The first one is the analyses of the convergence rate for the DE formula and the GL formula, the second one is a method of selecting appropriate stopping criteria, and the third one is the verification of the practical performance of the presented algorithms, when applied to large sparse matrices from current research problems.

Appendix A Effect of the parameter errors in Algorithm 1

In this section, in order to check the effects of the parameter errors in Algorithm 1, we present some numerical examples.

The convergence histories of the DE formula are shown in Figure 2.

Each graph shows two histories: one is obtained by using the eigvals and svdvals functions for computing $\rho(\tilde{A})$, $\|\tilde{A}-I\|_{2}$, and $\|\tilde{A}^{-1}\|_{2}$ as in Section 3; the other one is obtained by using Arpack.jl with tol=0.01333 The parameter tol defines the relative tolerance for convergence of Ritz values. In this examples, $\rho(\tilde{A})$ is obtained using eigs(A, nev=1, which=:LM, tol=0.01), $\|\tilde{A}-I\|_{2}$ is obtained using eigs((A-I)’*(A-I), nev=1, which=:LM, tol=0.01), and $\|\tilde{A}^{-1}\|_{2}$ is obtained using eigs(A’*A, nev=1, which=:SM, tol=0.01), where A is the preconditioned matrix $\tilde{A}$. . The figure shows that the behaviors of the histories are almost equal, and the effects of the errors did not appear.

In conclusion, the parameters used at Step 2 of Algorithm 1 can be computed roughly.

Appendix B Adaptive quadrature algorithm based on the GL quadrature

In this section we show an adaptive quadrature algorithm based on the GL quadrature for (3), using a technique from [2].

Let $G(m):=\sum_{i=1}^{m}w_{i}F_{\mathrm{GL}}(u_{i})$ be the results of the $m$-point GL quadrature, where $F_{\mathrm{GL}}(u_{i})=[(1+u)(A-I)+2I]^{-1}$. Using [2, Eq. (4.6)], the following error estimate can be applied:

[TABLE]

Based on (56), we present the algorithm designed to compute $\log(A)$ by automatically adding abscissas until the error is smaller than a given tolerance as Algorithm 3.

The computational cost of Algorithm 3 for dense $A$ is $[2(\sum_{i=0}^{k+1}m_{i})+2]n^{3}+\mathcal{O}(n^{2})=[(2^{k+3}-2)m_{0}+2]n^{3}+\mathcal{O}(n^{2})$, and the computational cost for $\log(A)\bm{b}$ with sparse $A$ is $(2^{k+2}-1)m_{0}c_{\mathrm{abscissa}}+c_{\mathrm{mul}}+c_{\mathrm{param}}$. Under the assumption that the total computational cost largely depends on the coefficient of $c_{\mathrm{abscissa}}$, the computational cost of the Algorithm 2 will be smaller than that of Algorithm 3 when the convergence ratios of the DE formula and the GL quadrature are the same.