Inexact rational Krylov method for evolution equations

Abstract

Linear and nonlinear evolution equations have been formulated to address problems in various fields of science and technology. Recently, methods using an exponential integrator for solving evolution equations, where matrix functions must be computed repeatedly, have been investigated and refined. In this paper, we propose a new method for computing these matrix functions which is called an inexact rational Krylov method. This is a more efficient version of the rational Krylov method with appropriate shifts, which was proposed by Hashimoto and Nodera (ANZIAM J 58:C149–C161, 2016). The advantage of the inexact rational Krylov method is that it computes linear equations that appear in the rational Krylov method efficiently while guaranteeing the accuracy of the final solution.

Appendices

A Proof of Proposition 3.1

We use the Householder reflectors for transforming the matrix $\tilde{K}_m$ into an upper Hessenberg matrix. Let $u_j:=(\tilde{k}_{j+1:m,j}-\eta _j e_1)/\Vert \tilde{k}_{j+1:m,j}-\eta _j e_1\Vert$ and let $\tilde{Q}_{j+1}:=-2u_ju_j^*$. Then, $I_{m-j}+\tilde{Q}_{j+1}$ is a unitary matrix and satisfies $(I_{m-j}+\tilde{Q}_{j+1})k_{j+1:m,j}=\eta _je_1$. Thus, the matrix $Q_m$ which is defined as $Q_m:=(I_m+\hat{Q}_{m-1})\cdots (I_m+\hat{Q}_2)$ is a unitary matrix, and there exists an upper Hessenberg matrix $H_m$ such that $Q_m\tilde{K}_mQ_m^*=H_m$. Here, $\hat{Q}_{j+1}:={\text {diag}}\{O_j,\ \tilde{Q}_{j+1}\}$ and $O_j$ is the $j\times j$ zero matrix.

Our first goal is to evaluate the magnitude of each element of the matrix $\tilde{Q}_{j+1}$. By the assumption (3.11), the inequality $\Vert \tilde{k}_{j+1:m,j}-\eta _j e_1\Vert \ge \eta$ holds. In addition, since the vector $\tilde{k}_{j+1:m,j}$ satisfies the condition (3.10), we have

$$\begin{aligned} \vert \eta _j\vert =\Vert \tilde{k}_{j+1:m,j}\Vert \le \sqrt{\sum _{k=1}^{\infty }(\hat{\alpha }\hat{\lambda }^k)^2}= \hat{\alpha }\sqrt{1/(1-\hat{\lambda }^2)}\hat{\lambda }. \end{aligned}$$

The first element of the vector $\tilde{k}_{j+1:m,j}-\eta _je_1$ is equal to $\tilde{k}_{j+1,j}-\eta _j$, and since the identity $\vert \eta _j\vert =\Vert \tilde{k}_{j+1:m,j}\Vert$ holds, we obtain $\vert \tilde{k}_{j+1,j}-\eta _j\vert \le 2\vert \eta _j\vert$. All the elements except for the first element of the vector $\tilde{k}_{j+1:m,j}-\eta _je_1$ are the same as the element $\tilde{k}_{j+1:m,j}$. For these reasons, on the basis of the assumption (3.10), $|u_{i,j}|<\tilde{\alpha }\hat{\lambda }^i$ is satisfied for

$$\begin{aligned} \tilde{\alpha }:=\max (\hat{\alpha }/\eta ,2\hat{\alpha } \sqrt{1/(1-\hat{\lambda }^2)}\hat{\lambda }/\eta ), \end{aligned}$$

where $u_{i,j}$ is the ith element of the vector $u_j\in \mathbb {C}^{m-j}$. Thus, we have

$$\begin{aligned} |(\tilde{Q}_{j+1})_{k,l}|=2|(u_ju_j^*)_{k,l}|=2|u_{k,j}u_{l,j}|\le 2\tilde{\alpha }^2\hat{\lambda }^{k+l}. \end{aligned}$$

Next, we evaluate the magnitude of each element of the matrix $Q_m$. Let $\check{\alpha }:=2\tilde{\alpha }^2$. For $l\ge 2$ and $l<k_1<k_2<\cdots <k_r$, we have

$$\begin{aligned} |(\hat{Q}_{k_r}\cdots \hat{Q}_{k_1}\hat{Q}_l)_{i,j}|&\le \displaystyle \frac{\check{\alpha }^{r+1}\hat{\lambda }^{-2(l-r-1)}}{(1-\hat{\lambda }^2)^r} \hat{\lambda }^{i+j}=\check{\alpha }^{r+1}\alpha ''(l,r)\hat{\lambda }^{i+j}\quad (i,j\le k_r),\nonumber \\ |(\hat{Q}_{k_r}\cdots \hat{Q}_{k_1}\hat{Q}_l)_{i,j}|&=0\quad (i>k_r\ \text {or}\ j>k_r), s \end{aligned}$$

(A.1)

where $\alpha ''(l,r):=\hat{\lambda }^{-2(l-r-1)}/(1-\hat{\lambda }^2)^r$. The inequality (A.1) is proved by the induction on r. For $r=1$, we have

$$\begin{aligned}&|(\hat{Q}_{k_1}\hat{Q}_l)_{i,j}| \le \sum _{a=k_1}^m\check{\alpha }\hat{\lambda }^{i-k_1+1+a-k_1+1} \check{\alpha }\hat{\lambda }^{a-l+1+j-l+1}\\&\qquad =\check{\alpha }^2\hat{\lambda }^{i+j}\hat{\lambda }^{-2(k_1+l-2)} \sum _{a=k_1}^m\hat{\lambda }^{2a}\le \frac{\check{\alpha }^2 \hat{\lambda }^{-2(l-2)}}{1-\hat{\lambda }^2}\hat{\lambda }^{i+j} \quad (i,j\le k_1). \end{aligned}$$

For $r\ge 2$, if the inequality (A.1) is satisfied with $r-1$, then we obtain

$$\begin{aligned}&|(\hat{Q}_{k_r}\cdots \hat{Q}_{k_1}\hat{Q}_l)_{i,j}| \le \sum _{a=k_r}^{m}\check{\alpha }\hat{\lambda }^{i-k_r+1+a-k_r+1} \frac{\check{\alpha }^r\hat{\lambda }^{-2(l-r)}}{(1-\hat{\lambda }^2)^{r-1}} \hat{\lambda }^{a+j}\\&\qquad \le \hat{\lambda }^{i+j}\frac{\check{\alpha }^{r+1} \hat{\lambda }^{-2(k_r+l-r-1)}}{(1-\hat{\lambda }^2)^{r-1}}\sum _{a=k_r}^m \hat{\lambda }^{2a} =\frac{\check{\alpha }^{r+1}\hat{\lambda }^{-2(l-r-1)}}{(1-\hat{\lambda }^2)^r}\hat{\lambda }^{i+j} \quad (i,j\le k_r), \end{aligned}$$

and the inequality (A.1) is also satisfied with r. This is the proof of the inequality (A.1). In the inequality (A.1), if $\hat{\lambda }\le 1/\sqrt{2}$ is satisfied, then the inequality $\alpha ''(l,r+1)\le \alpha ''(l,r)$ holds for any l. This results in $\alpha ''(l,r)\le \alpha ''(l,1)$ for any r and l. The matrix $Q_m$ can be represented as

$$\begin{aligned} Q_m&=(I_m+\hat{Q}_{m-1})\cdots (I_m+\hat{Q}_2)\nonumber \\&=I_m+\sum _{k=3}^{m-1}\sum _{l=2}^{k-1}\sum _{(a_1,a_2,\ldots , a_{k-l-1})\in \{0,1\}^{k-l-1}}\hat{Q}_{k}\hat{Q}_{k-1}^{a_1} \hat{Q}_{k-2}^{a_2}\cdots \hat{Q}_{l+1}^{a_{k-l-1}} \hat{Q}_l+\sum _{k=2}^{m-1}\hat{Q}_k. \end{aligned}$$

(A.2)

As a result, for $2\le \min {\{i,j\}}\le m-1$, the Eq. (A.2) and inequality (A.1) derive

$$\begin{aligned} |(Q_m-I_m)_{i,j}|&\le \check{\alpha }^2\sum _{k=3}^{\min {\{i,j\}}}\sum _{l=2}^{k-1} \left( 1+(k-l-1)\check{\alpha }+\left( \begin{array}{c}k-l-1\\ 2\end{array} \right) \check{\alpha }^2+\cdots +\check{\alpha }^{k-l-1}\right) \\&\qquad \cdot \alpha ''(l,1)\hat{\lambda }^{i+j} +\check{\alpha }^2\sum _{k=2}^{\min {\{i,j\}}}\alpha ''(k,1)\hat{\lambda }^{i+j}\\&\le \check{\alpha }^2\sum _{k=3}^{\min {\{i,j\}}}\sum _{l=2}^{k-1} (1+\check{\alpha })^{k-l-1} \alpha ''(l,1)\hat{\lambda }^{i+j}+\check{\alpha }^2 \sum _{k=2}^{\min {\{i,j\}}}\alpha ''(k,1)\hat{\lambda }^{i+j}. \end{aligned}$$

Therefore, for $2\le \min {\{i,j\}}\le m-1$ and $i\le j$, under the assumptions of (3.11) and (3.12), we have

$$\begin{aligned}&|(Q_m-I_m)_{i,j}|\\&\quad \le \sum _{k=3}^i\frac{(1+\check{\alpha })^{k-1}\check{\alpha }^2 \hat{\lambda }^4}{1-\hat{\lambda }^2} \hat{\lambda }^{i+j}\frac{((1+\check{\alpha })\hat{\lambda }^2)^{-2}}{((1+\check{\alpha })\hat{\lambda }^2)^{-1}-1}\left( ((1+\check{\alpha }) \hat{\lambda }^2)^{-k+2}-1\right) \\&\quad + \frac{\check{\alpha }^2\hat{\lambda }^4}{1-\hat{\lambda }^2}\hat{\lambda }^{i+j} \frac{\hat{\lambda }^{-4}}{\hat{\lambda }^{-2}-1}((\hat{\lambda }^{-2})^{i-1}-1)\\&\quad \le \frac{(1+\check{\alpha })\check{\alpha }^2\hat{\lambda }^4((1 +\check{\alpha })\hat{\lambda }^2)^{-1}}{(1-\hat{\lambda }^2)(1-(1+\check{\alpha })\hat{\lambda }^2)} \hat{\lambda }^{i+j}\sum _{k=3}^i(\hat{\lambda }^{-2})^{k-2} +\frac{\check{\alpha }^2\hat{\lambda }^4}{(1-\hat{\lambda }^2)^2} \hat{\lambda }^{i+j}\hat{\lambda }^{-2i}\\&\quad \le \frac{\check{\alpha }^2\hat{\lambda }^2}{(1-\hat{\lambda }^2)(1-(1+\check{\alpha })\hat{\lambda }^2)} \hat{\lambda }^{i+j}\frac{\hat{\lambda }^{-2}}{\hat{\lambda }^{-2}-1}(\hat{\lambda }^{-2})^{i-2} +\frac{\check{\alpha }^2\hat{\lambda }^4}{(1-\hat{\lambda }^2)^2} \hat{\lambda }^{i+j}\hat{\lambda }^{-2i}\\&\quad =\left( \frac{\check{\alpha }^2\hat{\lambda }^2}{(1-\hat{\lambda }^2) (1-(1+\check{\alpha })\hat{\lambda }^2)}\frac{\hat{\lambda }^4}{1-\hat{\lambda }^2} +\frac{\check{\alpha }^2\hat{\lambda }^4}{(1-\hat{\lambda }^2)^2}\right) \hat{\lambda }^{j-i}\\&\quad =:\alpha '\hat{\lambda }^{j-i}, \end{aligned}$$

where the sum $\sum _{k=3}^i$ becomes 0 for $i=2$, and we use the assumption (3.12) to derive the second inequality. In a similar manner, it can be deduced that $|(Q_m-I_m)_{i,j}|\le \alpha '\hat{\lambda }^{i-j}$ for $i>j$. In the case of $\min {\{i,j\}}=m$, the identity $i=j=m$ holds and we obtain

$$\begin{aligned} |(Q_m-I_m)_{m,m}|\le \check{\alpha }^2\sum _{k=3}^{m-1} \sum _{l=2}^{k-1}(1+\check{\alpha })^{k-l+1}\alpha ''(l,1) \hat{\lambda }^{i+j}+\check{\alpha }^2\sum _{k=2}^{m-1} \alpha ''(k,1)\hat{\lambda }^{i+j}\le \alpha '. \end{aligned}$$

For $\min {\{i,j\}}=1$, by the definition of the matrix $Q_m$, we have

$$\begin{aligned} \left\{ \begin{aligned} |(Q_m)_{1,1}|&=1\\ |(Q_m)_{i,j}|&=0\quad (i\ne 1\ \text {or}\ j\ne 1). \end{aligned} \right. \end{aligned}$$

Since the matrix $I_m$ is a diagonal matrix, redefining the constant $\alpha '$ as the sum of 1 and the previous $\alpha '$ completes the proof.

B Proof of Proposition 3.2

First, we apply Lemma 3.1 and Proposition 3.1 to the matrix $\phi _k\left( D_m-H_m^{-1}T_m\right) H_m^{-1}$. Since the matrix $H_m$ is an upper Hessenberg matrix satisfying the condition (3.13), by setting $f(z)=z^{-1}$ and applying Lemma 3.1, there exist constants $\hat{\alpha }>0$ and $0<\hat{\lambda }<1$ which satisfy the following inequality:

$$\begin{aligned} |(H_m^{-1})_{i,j}|\le \hat{\alpha }\hat{\lambda }^{i-j}\quad (i\ge j). \end{aligned}$$

(B.1)

Since the matrix $D_m$ is a diagonal matrix and the matrix $T_m$ is defined as the Eq. (2.5), redefining the constant $\hat{\alpha }$ as the sum of $\Vert D_m\Vert =N-h$ and the previous $\hat{\alpha }$ leads to

$$\begin{aligned} |(D_m-H_m^{-1}T_m)_{i,j}|\le \hat{\alpha }\hat{\lambda }^{i-j}\quad (i\ge j), \end{aligned}$$

where the constants $\hat{\alpha }$ and $\hat{\lambda }$ do not depend on m. Let

$$\begin{aligned} \mathscr {G}^{{\text {exp}}}(\alpha ,\lambda ):=\{A:\text {a square matrix}\mid |(A)_{i,j}| \le \alpha \lambda ^{|i-j|}\quad (\forall i,j)\}. \end{aligned}$$

By Proposition 3.1, there exists a unitary matrix $Q_m$ and an upper Hessenberg matrix $\tilde{H}_m$ such that $D_m-H_m^{-1}T_m=Q_m^*\tilde{H}_mQ_m$ and $Q_m\in \mathscr {G}^{{\text {exp}}}(\alpha ',\hat{\lambda })$, where $\alpha '>0$ is a constant which does not depend on m. Thus, it is deduced that there exists a matrix $\hat{H}_m\in \mathscr {G}^{{\text {exp}}}(\hat{\alpha },\hat{\lambda })$ such that

$$\begin{aligned} \phi _k\left( D_m-H_m^{-1}T_m\right) H_m^{-1}e_1&=\phi _k(Q_m^*\tilde{H}_mQ_m)H_m^{-1}e_1 =Q_m^*\phi _k(\tilde{H}_m)Q_m\hat{H}_me_1. \end{aligned}$$

The second equality holds, since by the inequality (B.1), there exists a matrix $\hat{H}_m\!\in \mathscr {G}^{{\text {exp}}}(\hat{\alpha },\hat{\lambda })$ which satisfies $H_m^{-1}e_1=\hat{H}_me_1$.

Our task is now to derive an upper bound of the magnitude of each element of the matrix $Q_m^*\phi _k(\tilde{H}_m)Q_m\hat{H}_m$, which is composed of the products of the matrices that have the decay properties. According to Benzi and Boito [2, Theorem 9.2], there exist constants $\alpha ''>0$ and $\lambda ''$ which do not depend on m and satisfy $Q_m\hat{H}_m\in \mathscr {G}^{{\text {exp}}}(\alpha '',\lambda '')$. In addition, because the function $\phi _k$ is an entire function, by setting $f=\phi _k$ in Proposition 3.1, there exist constants $\check{\alpha }>0$ and $0<\check{\lambda }<1$ satisfying $\vert (\phi _k(\tilde{H}_m))_{i,j}\vert \le \check{\alpha }\check{\lambda }^{i-j}$ for $i\ge j$. Let $\varSigma :=\bigcup _{m=1}^n W(D_m-H_m^{-1}T_m)$. Then, the set $\varSigma$ is closed and bounded, and the matrix $\tilde{H}_m$ satisfies

$$\begin{aligned} W(\tilde{H}_m)=W\left( Q_m(D_m-H_m^{-1}T_m)Q_m^*\right) =W(D_m-H_m^{-1}T_m) \subseteq \varSigma . \end{aligned}$$

(B.2)

Therefore, $|e^z|\le C'$ for a constant $C'>0$ is satisfied for $z\in W(\tilde{H}_m)$, and we have

$$\begin{aligned} |\phi _k(z)|=\left| \int _0^1e^{sz}\frac{(1-s)^{k-1}}{(k-1)!}ds \right| \le |e^z|\left| \int _0^1\frac{(1-s)^{k-1}}{(k-1)!}ds \right| \le \frac{C'}{k!}\quad (z\in W(\tilde{H}_m)). \end{aligned}$$

(B.3)

Using the result by Crouzeix [7], the condition (B.2), and the inequality (B.3), there exists a constant $2\le C\le 1+\sqrt{2}$ such that

$$\begin{aligned} \Vert \phi _k(\tilde{H}_m)\Vert \le C\sup _{z\in W(\tilde{H}_m)}|\phi _k(z)|\le \frac{CC'}{k!}. \end{aligned}$$

Redefining the constant $\check{\alpha }$ as the sum of the constant $CC'/(k!)$ and the previous $\check{\alpha }$ results in

$$\begin{aligned}&\left| \left( \phi _k(\tilde{H}_m)\right) _{i,j}\right| \le \check{\alpha }\check{\lambda }^{i-j}\quad (i\ge j), \end{aligned}$$

(B.4)

$$\begin{aligned}&\left| \left( \phi _k(\tilde{H}_m)\right) _{i,j}\right| \le \Vert \phi _k(\tilde{H}_m)\Vert \le \check{\alpha }\quad (i<j). \end{aligned}$$

(B.5)

By the upper bounds (B.4) and (B.5), we obtain

$$\begin{aligned}&\left| \left( \phi _k(\tilde{H}_m)Q_m\hat{H}_m\right) _{i,1}\right| \le \sum _{k=1}^{i}\check{\alpha }\check{\lambda }^{i-k}\alpha ''\lambda ''^{k-1} +\sum _{k=i+1}^m\check{\alpha }\alpha '' \lambda ''^{k-1}\nonumber \\&\qquad \le i\check{\alpha }\alpha ''\bar{\lambda }^{i-1}+\check{\alpha }\alpha '' \frac{\lambda ''^i}{1-\lambda ''}\le i\check{\alpha }\alpha ''\left( 1+\frac{\lambda ''}{1-\lambda ''}\right) \bar{\lambda }^{i-1} =i\bar{\alpha }\bar{\lambda }^{i-1}, \end{aligned}$$

(B.6)

where $\bar{\alpha }:=\check{\alpha }\alpha ''/(1-\lambda '')$ and $\bar{\lambda }:=\max \{\check{\lambda },\lambda ''\}<1$. As a result, using the fact $Q_m\in \mathscr {G}^{{\text {exp}}}(\alpha ',\hat{\lambda })$ and the upper bound (B.6), we have

$$\begin{aligned}&\left| \left( \phi _k(D_m-H_m^{-1}T_m)H_m^{-1}\right) _{i,1}\right| =\left| \left( Q_m^*\phi _k(\tilde{H}_m)Q_m\hat{H}_m\right) _{i,1}\right| \\&\quad \le \sum _{k=1}^{i}\alpha '\hat{\lambda }^{i-k}k\bar{\alpha }\bar{\lambda }^{k-1} +\sum _{k=i+1}^{m}\alpha '\hat{\lambda }^{k-i}k\bar{\alpha }\bar{\lambda }^{k-1}\\&\quad \le \frac{1}{2}(i+1)i\alpha '\bar{\alpha }\lambda ^{i-1}+\alpha '\bar{\alpha } \frac{i+1}{(1-\lambda ^2)^2}\lambda ^{i+1}\\&\quad \le \frac{1}{2}(i+1)i\alpha '\bar{\alpha }\left( 1+ \frac{2\lambda ^2}{(1-\lambda ^2)^2}\right) \lambda ^{i-1} =\frac{1}{2}(i+1)i\alpha \lambda ^{i-1}, \end{aligned}$$

where $\alpha :=\alpha '\bar{\alpha }(1+2\lambda ^2/(1- \lambda ^2)^2)$ and $\lambda :=\max \{\hat{\lambda },\bar{\lambda }\}<1$. This completes the proof of Proposition 3.2.