rational Krylov, parameterized linear systems, companion linearization, shifted linear systems, inexact Krylov, hierarchical matrices, boundary element method, Helmholtz equation

1 Introduction↩︎

The motivation of this work arises from the field of integral equations. The Boundary Element Method (BEM) is a popular numerical technique for solving scattering problems in the frequency domain, based on an integral equation reformulation of the governing partial differential equation, such as the Helmholtz equation. Among many advantages of this approach, the resulting BEM system matrix is generally dense and nonlinearly dependent on the frequency. The dense matrix can still be approximated in a data-sparse way, using a range of algebraic and analytic techniques [1]–[8] and the resulting systems can be solved efficiently with iterative or direct solvers [5], [9]–[13].

Still, the issue of nonlinear frequency dependency remains. In a ‘frequency sweep’, the solution to a scattering problem is sought over a range of frequencies. While the above techniques speed up solution at a single frequency, repeating such procedures for all frequencies of interest may still be intractable. Earlier research on frequency sweeping has involved both the Finite Element Method (FEM) and BEM, see [14] for a recent survey. Frequency sweeping requires tackling both the efficient construction and representation of the frequency-dependent system, e.g.[15], [16], as well as efficiently resolving the solutions over the frequency range, e.g. [17]–[19].

The wavenumber-dependent system for Helmholtz BEM leads to a parametrized linear system, given in general by \[\label{eq:parametrized95linear95system} \mathsf{P}(\mu)\mathsf{x}(\mu) = \mathsf{b}(\mu)\tag{1}\] with \(\mu\in S\subseteq\mathbb{C}\) and \(\mathsf{P}: S \to \mathbb{C}^{n\times n}\), where \(\mathsf{x}(\mu)\) is sought for many values of \(\mu\). The study of parametrized linear systems is rich, see [20]–[23] and references therein.

Using companion linearization techniques, we formulate the BEM system as a set of shifted linear systems, linear in the parameter \(\mu\), with a matrix pencil \((\boldsymbol{\mathsf{A}},\boldsymbol{\mathsf{B}})\) as \[\label{eq:shifted95linear95system} (\boldsymbol{\mathsf{A}}-\mu \boldsymbol{\mathsf{B}})\boldsymbol{\mathsf{y}}= \boldsymbol{\mathsf{d}}.\tag{2}\] This involves the approximation of \(\mathsf{P}(\mu)\) using a basis, for which we use rational functions. The linearized system can be solved for many different values of \(\mu\) simultaneously using Krylov subspace methods [24], but it has (much) larger dimensions than the original one. We take inspiration from the compact rational Krylov (CORK) framework, originally formulated for nonlinear eigenvalue problems [25], to exploit the rich algebraic structure of the system and speed up the application of rational Krylov.

1.1 Related work↩︎

Relevant related work for parametrized systems was described in [23], [26], which combined a linearization using Taylor series expansions with a GMRES solver, and in [27], using Chebyshev interpolation of \(\mathsf{P}(\mu)\) in combination with BiCG. However, rational Krylov was not used in these references. In comparison, the CORK framework enables general companion linearizations, such as those associated with rational functions, and owing to rational Krylov it also allows for improved convergence of the solutions through the use of multiple shifts.

Rational Krylov for shifted systems was investigated in [28], leading to a scheme similar to the one we propose using left preconditioning, but for the special case \(\boldsymbol{\mathsf{B}}=\mathsf{I}\) and without initial shift. Rational Krylov for shifted systems with right preconditioning is equivalent to the so-called ‘flexible Krylov’ techniques of [29], [30]. We augment both settings with the CORK framework for faster computations.

1.2 Results↩︎

The contributions of this article are the following. We:

• demonstrate how the CORK framework can be used to solve parametrized linear systems with a general pencil \((\boldsymbol{\mathsf{A}},\boldsymbol{\mathsf{B}})\) for a range of the parameter \(\mu\);

• introduce a right-preconditioned variant of CORK;

• allow parameter dependence of the right-hand side of 1 ;

• incorporate inexact Krylov theory [31], [32] into CORK;

• apply the framework, combined with the compact wavenumber-dependent BEM representation from [16], [33], [34], to efficiently perform frequency sweeps.

The frequency sweeping method is valid for a continuous frequency range, based on solves at a small number of frequency values. They correspond to shifts of the rational Krylov method, and can be greedily chosen very effectively based on a cheaply available estimate of the residual across the range. This is a highly attractive feature of the scheme. Afterwards, each additional solve requires only the solution of a very small linear system, whose size relates to that of the constructed Krylov subspace. We show that the method is very performant for small to moderately large frequencies.

1.3 Outline↩︎

We summarize rational Krylov GMRES for shifted linear systems in §2, including a variant involving right preconditioning and a greedy shift selection procedure. We introduce compact rational Krylov (CORK) for parameterized linear systems with varying right hand sides in §3. We integrate inexact Krylov theory in §4. The methodogy is applied to BEM frequency sweeping problems in §5 and further illustrate the approach numerically in §6.

2 Rational Krylov GMRES for shifted linear systems↩︎

The rational Krylov method was developed in [35], [36] and generalizes shift-and-invert Arnoldi for the solution of (generalized) eigenvalue problems. This is achieved by allowing multiple different shifts throughout the iterative process.

Rational Krylov for matrix pencils works as follows. Given a starting vector \(\boldsymbol{\mathsf{v}}_1\), a matrix pencil \((\boldsymbol{\mathsf{A}},\boldsymbol{\mathsf{B}})\) and a set of shifts \(\boldsymbol{\sigma}_m=(\sigma_i)_{i=1}^{m}\), the algorithm constructs the rational Krylov subspace³ \[\label{eq:ratkrylov-subspace} \mathcal{K}_{m+1} = \mathrm{span}{ \left\{\boldsymbol{\mathsf{v}}_1,(\boldsymbol{\mathsf{A}}-\sigma_1\boldsymbol{\mathsf{B}})^{\smin 1}\boldsymbol{\mathsf{B}}\boldsymbol{\mathsf{v}}_1,\dots,\prod_{i=1}^{m}\left((\boldsymbol{\mathsf{A}}-\sigma_i\boldsymbol{\mathsf{B}})^{\smin 1}\boldsymbol{\mathsf{B}}\right)\boldsymbol{\mathsf{v}}_1 \right\} }\tag{3}\] and builds an orthogonal basis \(\boldsymbol{\mathsf{V}}_{m+1}:=[\boldsymbol{\mathsf{v}}_1\cdots\boldsymbol{\mathsf{v}}_{m+1}]\) satisfying the recurrence relation \[\label{eq:rat-krylov-reccurence} \boldsymbol{\mathsf{A}}\boldsymbol{\mathsf{V}}_{m+1} \underline{\mathsf{H}}_m = \boldsymbol{\mathsf{B}}\boldsymbol{\mathsf{V}}_{m+1} \underline{\mathsf{K}}_m\tag{4}\] with \(\underline{\mathsf{H}}_m,\underline{\mathsf{K}}_m\in\mathbb{C}^{(m+1)\times m}\) both upper Hessenberg. Letting \(\mathsf{\Sigma}_m:=\mathrm{diag}{(\boldsymbol{\sigma}_m)}\), \(\underline{\mathsf{K}}_m\) is related to the shifts as \(\underline{\mathsf{K}}_m=\underline{\mathsf{H}}_m\mathsf{\Sigma}_m+\underline{\mathsf{T}}_m\), in which upper triangular matrix \(\underline{\mathsf{T}}_m\) contains the selected continuation combinations, see [36] and [37].

Rational Krylov methods are typically applied to eigenvalue problems. Approximations to eigenpairs of \((\boldsymbol{\mathsf{A}},\boldsymbol{\mathsf{B}})\) can be found from the eigenpairs of the pencil \((\mathsf{K}_m,\mathsf{H}_m)\). The approximate eigenvectors live in the subspace spanned by \(\boldsymbol{\mathsf{V}}_{m+1}\underline{\mathsf{H}}_m\).

Here, instead, we aim to solve shifted linear systems 2 with pencil \((\boldsymbol{\mathsf{A}},\boldsymbol{\mathsf{B}})\). When trying to find eigenvalues, rational Krylov chooses the shifts \(\boldsymbol{\sigma}_m\) in the region of interest to steer convergence. Similarly, we can choose the shifts close to values of \(\mu\) for which we want to solve. In this way, we hope the rational Krylov subspace contains good approximations to the solutions for all these values of \(\mu\).

We discuss the use of rational Krylov as an algorithm to solve shifted systems. This leads first to left-preconditioned rational Krylov GMRES (LRK-GMRES).

2.1 Left-preconditioned RK-GMRES↩︎

A shift-and-invert strategy can be combined with GMRES to solve shifted linear systems. One chooses a shift \(\sigma\) in the range of interest for \(\mu\), with associated preconditioner \((\boldsymbol{\mathsf{A}}-\sigma\boldsymbol{\mathsf{B}})^{\smin 1}\). The proposition below, due to [38], formalizes how to exploit the shift-invariance of Krylov subspaces to solve the shifted system for any value of \(\mu\). The solution is found from a small least-squares system, whose accuracy relates to the preconditioned residual.

In rational Krylov, instead of applying \((\boldsymbol{\mathsf{A}}-\sigma \boldsymbol{\mathsf{B}})^{\smin 1}\boldsymbol{\mathsf{B}}\) in each iteration for a fixed \(\sigma\), the shift varies. This results in recurrence relation 4 instead of 5 . The following proposition can be seen to generalize the rational Krylov subspace method for shifted linear systems from [28]; the latter is obtained by taking \(\boldsymbol{\mathsf{B}}=\mathsf{I}\) and \(\sigma_0=\infty\).

An initial shift \(\sigma_0\) remains in the scheme and appears in the preconditioned residual. Hence, its choice plays a role in the convergence of the scheme.

2.2 Right-preconditioned RK-GMRES↩︎

The question arises whether the choice of an initial shift can be avoided. A shift-invariant Krylov subspace can also be obtained using right-preconditioning instead. As it turns out, this setting is exactly the same as that of ‘flexible GMRES for shifted systems’ [30], see also the flexible GMRES algorithm of [29]. For this reason, we shall be brief in the description here.⁴

Right preconditioning for the solution of \((\boldsymbol{\mathsf{A}}-\mu \boldsymbol{\mathsf{B}})\boldsymbol{\mathsf{y}}= \boldsymbol{\mathsf{d}}\) means \[(\boldsymbol{\mathsf{A}}-\mu \boldsymbol{\mathsf{B}})(\boldsymbol{\mathsf{A}}-\sigma \boldsymbol{\mathsf{B}})^{\smin 1}\bar{\boldsymbol{\mathsf{y}}} = \boldsymbol{\mathsf{d}}, \quad with \quad \boldsymbol{\mathsf{y}}= (\boldsymbol{\mathsf{A}}-\sigma \boldsymbol{\mathsf{B}})^{\smin 1}\bar{\boldsymbol{\mathsf{y}}}.\] This is equivalent with \(\left(\mathsf{I}+ (\sigma-\mu)\boldsymbol{\mathsf{B}}(\boldsymbol{\mathsf{A}}-\sigma \boldsymbol{\mathsf{B}})^{\smin 1}\right)\bar{\boldsymbol{\mathsf{y}}} = \boldsymbol{\mathsf{d}}\) and suggests applying Arnoldi on \(\boldsymbol{\mathsf{B}}(\boldsymbol{\mathsf{A}}-\sigma \boldsymbol{\mathsf{B}})^{\smin 1}\) instead of \((\boldsymbol{\mathsf{A}}-\sigma \boldsymbol{\mathsf{B}})^{\smin 1}\boldsymbol{\mathsf{B}}\). With multiple shifts, this leads to \[\label{eq:ratkrylov-subspace-right} \mathcal{K}_{m+1}^{\textrm{right}} = \mathrm{span}{ \left\{\boldsymbol{\mathsf{v}}_1,\boldsymbol{\mathsf{B}}(\boldsymbol{\mathsf{A}}-\sigma_1\boldsymbol{\mathsf{B}})^{\smin 1}\boldsymbol{\mathsf{v}}_1,\dots,\prod_{i=1}^{m}\left(\boldsymbol{\mathsf{B}}(\boldsymbol{\mathsf{A}}-\sigma_i\boldsymbol{\mathsf{B}})^{\smin 1}\right)\boldsymbol{\mathsf{v}}_1 \right\} }.\tag{12}\]

This time around, we can use the starting vector \(\boldsymbol{\mathsf{v}}_1:=\boldsymbol{\mathsf{d}}/\|\boldsymbol{\mathsf{d}}\|\). In this context a matrix \(\boldsymbol{\mathsf{Z}}_m\) arises due to the presence of the flexible right preconditioner. It satisfies \[\boldsymbol{\mathsf{B}}\boldsymbol{\mathsf{Z}}_m = \boldsymbol{\mathsf{V}}_{m+1}\underline{\mathsf{H}}_m, \quad and \quad \boldsymbol{\mathsf{A}}\boldsymbol{\mathsf{Z}}_m - \boldsymbol{\mathsf{B}}\boldsymbol{\mathsf{Z}}_m\mathsf{\Sigma}_m = \boldsymbol{\mathsf{V}}_m\mathsf{T}_m.\] Multiplying the first equation by \(\mathsf{\Sigma}_m-\mu \mathsf{I}_m\) shows that \[\label{eq:rrkgmres-shifted-recurrence} (\boldsymbol{\mathsf{A}}-\mu \boldsymbol{\mathsf{B}})\boldsymbol{\mathsf{Z}}_m = \boldsymbol{\mathsf{V}}_{m+1}\left(\underline{\mathsf{T}}_m + \underline{\mathsf{H}}_m(\mathsf{\Sigma}_m-\mu \mathsf{I}_m)\right) = \boldsymbol{\mathsf{V}}_{m+1}(\underline{\mathsf{K}}_m - \mu \underline{\mathsf{H}}_m).\tag{13}\] Thus, the approximate solution lies in the column space of \(\boldsymbol{\mathsf{Z}}_m\). It is obtained as \(\boldsymbol{\mathsf{y}}_m(\mu)=\boldsymbol{\mathsf{Z}}_m\mathsf{y}^\mathrm{red}_m(\mu)\), where \(\mathsf{y}^\mathrm{red}_m\) solves 8 as before.

One advantage of right preconditioning is that the residual of the projected least-squares system, which is minimized, is equal to the residual of the original system, \[\label{eq:rrk95residual} \boldsymbol{\mathsf{r}}_m = \boldsymbol{\mathsf{d}}- (\boldsymbol{\mathsf{A}}-\mu \boldsymbol{\mathsf{B}})\boldsymbol{\mathsf{y}}_m = \beta\mathsf{e}_1-(\underline{\mathsf{K}}_m-\mu\underline{\mathsf{H}}_m)\mathsf{y}_m^\mathrm{red},\tag{14}\] \(\beta=\|\boldsymbol{\mathsf{d}}\|\), whereas with left preconditioning one obtains a preconditioned residual. On the other hand, the matrix \(\boldsymbol{\mathsf{Z}}_m\) needs to be stored in addition to \(\boldsymbol{\mathsf{V}}_{m+1}\).

2.3 Greedy selection of shifts in rational Krylov↩︎

The distribution of the shifts in rational Krylov relative to the values of \(\mu\) for which one wants to solve the system determines the effectiveness of the approach. It can generally be expected that the introduction of a shift \(\sigma_k\) at iteration \(k\) is most beneficial for nearby values of \(\mu\). Logically, one may choose the next shift where the error is currently large.

One of the compelling features of rational Krylov is that a cheap estimate of the residual norm is available for all values of \(\mu\), precisely because one can solve a small least-squares system. That system involves an upper Hessenberg matrix \(\underline{\mathsf{K}}_k-\mu\underline{\mathsf{H}}_k\), which can be reduced to upper triangular form in \(\mathcal{O}(k^2)\) operations. A practical strategy is to choose a finite subset \(\hat{S} \subset S\) of values for \(\mu\) in its range of interest. A straightforward greedy approach consists in choosing the next shift as the maximizer of the residual norm over \(\hat{S}\). For right and left-preconditioned RK-GMRES this means \[\sigma_{k+1} = \mathop{\mathrm{arg\,max}}_{\mu\in \hat{S}}{\|\boldsymbol{\mathsf{r}}_k(\mu)\|}, \quad or \quad \sigma_{k+1} = \mathop{\mathrm{arg\,max}}_{\mu\in \hat{S}}{\|\boldsymbol{\mathsf{r}}_k^{\sigma}(\mu)\|}.\] All that remains is to select an initial shift \(\sigma_1\) and, for LRK-GMRES, also \(\sigma_0\).

Such a greedy scheme was recently investigated in [28], in which shifts are referred to as ‘poles’. It was compared with the so-called ADM shift selection strategy of [39], which is based on a representation of the residual norm \(\|\boldsymbol{\mathsf{r}}_k(\mu)\|\) as a rational function in \(\mu\). It was found that the simple greedy strategy performs significantly better.

The way we apply the greedy scheme differs in two ways from that of [28]. First, in the setting of parametrized linear systems we are primarily interested in the residual of the parametrized system, not of its linearization. The second difference is that we use inexact iterative solvers at each shift, which means that the computed residuals no longer match the true ones. These are the topics of the next two sections.

3 Compact rational Krylov for parametrized linear systems↩︎

We move from rational Krylov as a means to solve families of shifted linear systems to its application for the solution of parametrized linear systems 1 . The exposition is based on the original paper [25], in which compact rational Krylov (CORK) was proposed for nonlinear eigenvalue problems. As mentioned before, other methods have been described using linearization techniques for the solution of parametrized linear systems, see [23], [26], [27]. In comparison to those, CORK accommodates general companion linearizations rather than specific choices. It also permits the use of multiple shifts to steer convergence, along with a cheap greedy selection strategy.

3.1 Linearization and left-preconditioned CORK-GMRES↩︎

The starting point is a matrix-valued expansion for \(\mathsf{P}(\mu)\) in a basis, say of the form \(\mathsf{P}(\mu) = \sum_{i=1}^s \mathsf{C}_i g_i(\mu)\). This is always possible with \(s=n^2\) basis functions \(g_i\), and often with fewer, but these functions are application-specific and may be unstructured.

Linearization requires a basis with more structure. Say \(d\) functions \(f_i: \mathbb{C}\to \mathbb{C}\), with \(\mathsf{f}(\mu)={[f_1(\mu) \,\cdots\, f_d(\mu)]}^\mathsf{T}\), satisfy the linear relationship \[\label{eq:linear-relation} (\mathsf{M}-\mu\mathsf{N})\mathsf{f}(\mu) = \mathsf{0},\tag{15}\] with \(\mathsf{M},\mathsf{N}\in\mathbb{C}^{(d-1)\times d}\) and \(\mathrm{rank}{(\mathsf{M}-\mu\mathsf{N})} = d-1\). Examples include polynomials and rational functions. If the spaces of the bases are the same, we obtain \[\label{eq:matrix-function-add} \mathsf{P}(\mu) = \sum_{i=1}^d (\mathsf{A}_i-\mu \mathsf{B}_i)f_i(\mu), \qquad \mathsf{A}_i,\mathsf{B}_i\in\mathbb{C}^{n\times n}.\tag{16}\] We refer to [25] for several examples. If the spaces differ, \(\mathsf{P}(\mu)\) can also be approximated using polynomials or rational functions, e.g.using AAA [40], [41].

Such linearization of \(\mathsf{P}\) leads to a large matrix pencil \(\boldsymbol{\mathsf{A}}-\mu\boldsymbol{\mathsf{B}} \in \mathbb{C}^{dn\times dn}\), with \[\label{eq:linearization} \boldsymbol{\mathsf{A}}= \left[\begin{array}{cccc} \mathsf{A}_1 & \mathsf{A}_2 & \cdots & \mathsf{A}_d \\ \noalign{\vskip 2pt} \hline \noalign{\vskip 2pt} \multicolumn{4}{c}{\mathsf{M}\otimes\mathsf{I}_n } \end{array}\right],\quad\quad \boldsymbol{\mathsf{B}}= \left[\begin{array}{cccc} \mathsf{B}_1 & \mathsf{B}_2 & \cdots & \mathsf{B}_d \\ \noalign{\vskip 2pt} \hline \noalign{\vskip 2pt} \multicolumn{4}{c}{\mathsf{N}\otimes\mathsf{I}_n } \end{array}\right].\tag{17}\] The following theorem shows the connection between the shifted linear system \(\boldsymbol{\mathsf{A}}-\mu\boldsymbol{\mathsf{B}}\) and the parametrized linear system \(\mathsf{P}(\mu)\mathsf{x}=\mathsf{b}\) (with constant right hand side). A more general, but also more abstract, statement was given in [20].

1 reveals that by using the pencil \((\boldsymbol{\mathsf{A}},\boldsymbol{\mathsf{B}})\), together with the appropriate right-hand side, we can employ methods for shifted linear systems to solve the parametrized linear systems 1 . We shall consider variable right hand sides in §3.3.

For nonlinear eigenvalue problems [25] CORK finds eigenpairs by applying the rational Krylov method [35], [36] on the companion pencil \((\boldsymbol{\mathsf{A}},\boldsymbol{\mathsf{B}})\). The CORK framework exploits the Kronecker structure present in the large pencil to efficiently employ Krylov methods. See [25], [34] for more details. Because CORK uses the rational Krylov subspace 3 , it corresponds to LRK-GMRES from §2.1. We will further refer to it as left-preconditioned CORK-GMRES (LCORK-GMRES).

3.2 Right-preconditioned CORK-GMRES (RCORK-GMRES)↩︎

In this setting we mimick the strategy of CORK and aim to exploit the Kronecker structure in the companion pencil for right preconditioning as well. For the sake of conciseness, the structure of this subsection parallels that of section 4 in [25]. We also refer to [34], which includes an explicit overview of steps that are only implicit in literature as well as some material omitted here for brevity.

We start with restating Theorem 2.3 from [25].

In addition to the pencil having block structure, so do the matrices \(\boldsymbol{\mathsf{V}}_{m+1}\) and, in the right-preconditioned case, \(\boldsymbol{\mathsf{Z}}_{m}\). With blocks \(\mathsf{V}_m^{[i]}\in\mathbb{C}^{n\times m}\) and \(\mathsf{v}_{m+1}^{[i]}\in\mathbb{C}^n\), for \(i=1,\dots,d\), we write \(\boldsymbol{\mathsf{V}}_{m+1}\in\mathbb{C}^{dn\times(m+1)}\) as \[\boldsymbol{\mathsf{V}}_{m+1} = \left[\begin{array}{cc} \boldsymbol{\mathsf{V}}_m & \boldsymbol{\mathsf{v}}_{m+1} \end{array}\right] = \left[\begin{array}{cc} \mathsf{V}_m^{[1]} & \mathsf{v}_{m+1}^{[1]} \\ \vdots & \vdots \\ \mathsf{V}_m^{[d]} & \mathsf{v}_{m+1}^{[d]} \end{array}\right].\] We consider a similar subdivision for \(\boldsymbol{\mathsf{Z}}_m\in\mathbb{C}^{dn\times m}\) into \(\mathsf{Z}_m^{[i]}\in\mathbb{C}^{n\times m}\) (\(i=1,\dots,d\)).

Using 1, \(\boldsymbol{\mathsf{V}}_{m+1}\) and \(\boldsymbol{\mathsf{Z}}_m\) can be expressed as \[\boldsymbol{\mathsf{V}}_{m+1} = (\mathsf{I}_d \otimes\mathsf{Q}_{m+1})\boldsymbol{\mathsf{U}}_{m+1},\quad \boldsymbol{\mathsf{Z}}_m = (\mathsf{I}_d \otimes\mathsf{Q}_{m+1})\boldsymbol{\mathsf{S}}_m,\] \[\text{with}\qquad \boldsymbol{\mathsf{U}}_m = \left[\begin{array}{c} \mathsf{U}_m^{[1]} \\ \vdots \\ \mathsf{U}_m^{[d]} \end{array}\right] \in \mathbb{C}^{dr_m\times m},\quad \boldsymbol{\mathsf{S}}_m = \left[\begin{array}{c} \mathsf{S}_m^{[1]} \\ \vdots \\ \mathsf{S}_m^{[d]} \end{array}\right] \in \mathbb{C}^{dr_{m+1}\times m}.\] This leads to the introduction of the CORK quintuple; compare with the CORK quadruple in [25].

Through additional theorems we show that this representation is compact in the sense that its storage scales as \(\mathcal{O}(nm+dm^2)\) instead of \(\mathcal{O}(dnm)\). The compact block structure is inherited from that of the solution of the shifted system.

When the shifted system \((\boldsymbol{\mathsf{A}}-\mu\boldsymbol{\mathsf{B}})\boldsymbol{\mathsf{y}}=\boldsymbol{\mathsf{d}}\) is used to solve the parametrized linear system \(\mathsf{P}(\mu)\mathsf{x}=\mathsf{b}\), we have \(\boldsymbol{\mathsf{d}}=\mathsf{e}_1\otimes\mathsf{b}\) and \(\tilde{d}=1\). In general, with \(c=2\), RCORK-GMRES requires roughly twice the amount of storage to store its vectors compared to LCORK-GMRES. This is similar to how regular RRK-GMRES requires twice the storage compared to LRK-GMRES by having to store \(\boldsymbol{\mathsf{Z}}_m\).

In [34], the main author further details how the block ULP decomposition combined with the compact representation of \(\boldsymbol{\mathsf{V}}_{m+1}\) and \(\boldsymbol{\mathsf{Z}}_m\) leads to efficient implementations of both LCORK-GMRES and RCORK-GMRES. The computationally intensive parts are the ones that scale with the problem size \(n\). They are, in iteration \(k\):

1. Solving a system involving \(\mathsf{P}(\sigma_k)\) with shift \(\sigma_k\).

2. A matrix-vector product with \((\bar{\boldsymbol{\mathsf{A}}}_2-\sigma_k\bar{\boldsymbol{\mathsf{B}}}_2)\in\mathbb{C}^{(d-1)n\times n}\) and one with \(\bar{\boldsymbol{\mathsf{B}}}\in\mathbb{C}^{dn\times n}\). This can be recast to a single matrix-vector product with \(\bar{\boldsymbol{\mathsf{C}}}:=[\mathsf{C}_1 \,\cdots\, \mathsf{C}_s]\in\mathbb{C}^{sn\times n}\) in the case of LCORK-GMRES. RCORK-GMRES requires two separate products with \(\bar{\boldsymbol{\mathsf{C}}}\) (unless \(\bar{\boldsymbol{\mathsf{B}}}=\mathsf{b}^\mathsf{T}\otimes\mathsf{I}_n\)).

3. Expanding \(\mathsf{Q}_k\) to \(\mathsf{Q}_{k+1}\) by means of orthogonalization. It is expanded by one vector for LCORK-GMRES and by two for RCORK-GMRES (unless \(\bar{\boldsymbol{\mathsf{B}}}=\mathsf{b}^\mathsf{T}\otimes\mathsf{I}_n\)).

The cost of the orthogonalization depends on the choice of orthogonalization process, but generally results in \(\mathcal{O}(nm^2)\) operations over \(m\) iterations. In contrast, the cost of the other operations mostly depends on the form of \(\mathsf{P}(\mu)\). Note that when \(\bar{\boldsymbol{\mathsf{B}}}=\mathsf{b}^\mathsf{T}\otimes\mathsf{I}_n\), RCORK-GMRES practically has the same cost as LCORK-GMRES. Earlier, the same was also true for their memory footprint.

3.3 Parameter dependency of the right-hand side↩︎

At first, a parameter-dependent right hand side \(\mathsf{b}:=\mathsf{b}(\mu)\) seems like a major obstacle. Krylov methods that rely on the shift-invariance property to solve shifted linear systems are not designed to handle this situation, since there is no single starting vector that is uniformly valid in \(\mu\). To the best of our knowledge, this case is not treated in the literature on linearization.

Fortunately, there is an original solution to this problem, which is also relatively simple. We can make a companion-like extension of the system to a slightly larger one, with size \((n+1)\times(n+1)\), with a constant right-hand side. To that end, we assume that \(\mathsf{b}\) can be represented with the same scalar functions as \(\mathsf{P}\), i.e., \[\label{eq:cond-param-dep-rhs} \mathsf{b}(\mu) = \sum_{i=1}^d (\mathsf{b}_{\mathsf{A},i} - \mu \mathsf{b}_{\mathsf{B},i}) f_i(\mu).\tag{18}\] We aim to construct an augmented system \(\underline{\mathsf{P}}(\mu)\underline{\mathsf{x}}(\mu)=\underline{\mathsf{b}}(\mu)\) of which the right-hand side is only dependent on \(\mu\) through a scalar. We do so using a scalar function \(\gamma(\mu)=\sum_{i=1}^d (\alpha_i-\mu\beta_i)f_i(\mu)\) and defining \[\underline{\mathsf{P}}(\mu):=\left[\begin{array}{cc} \mathsf{P}(\mu) & -\mathsf{b}(\mu) \\ \mathsf{0}_{1\times n} & \gamma(\mu) \end{array}\right], \quad \underline{\mathsf{b}}(\mu) := \gamma(\mu)\mathsf{e}_{n+1}=\left[\begin{array}{c} 0 \\ \gamma(\mu) \end{array}\right].\] where \(\underline{\mathsf{x}}(\mu) = [\mathsf{x}(\mu)^\mathsf{T}\; 1]^\mathsf{T}\). The augmented system results in a linearization with the original matrices \(\{\mathsf{A}_i\}_{i=1}^d\) and \(\{\mathsf{B}_i\}_{i=1}^d\) also augmented to slightly larger ones. Importantly, we can perform CORK-GMRES with starting vector \(\mathsf{e}_{n+1}\). The solution to the original system is obtained from the first \(n\) elements of the augmented solution, scaled by \(\gamma(\mu)\).

Any scalar function \(\gamma(\mu)\) of the right form will do, but generally one can choose the coefficients \(\alpha_i\) and \(\beta_i\) such that \(\gamma(\mu)\) is a constant. The choice could impact the stability of CORK. Strictly speaking, we have also changed the meaning of the residual of the system, which is a criterion for the greedy selection of shifts. Still, we found the approach to be very effective in the application of this paper.

4 Inexact rational Krylov↩︎

Thus far we have assumed that the systems \(\mathsf{P}(\sigma)\mathsf{x}=\mathsf{b}\) can be solved exactly. This is often not the case, or it might be impractical, for large-size problems. Instead, iterative methods are used to obtain an approximate solution with a suitably small residual norm. This leads to inexact Krylov methods, for which there is a rich body of literature [42]–[52].

The choice of the tolerance with which to solve the systems affects accuracy and computational cost. It was shown in [31] that one can loosen the tolerances of inexact matrix-vector products in GMRES, CG and related methods, throughout the iterations. More specifically, the required tolerances can be inversely proportional to the computed residual norm of the previous iteration. This observation also applies to shift-and-invert Arnoldi for shifted linear systems [31].

We adopt the strategy for RK-GMRES, similar to [31], [32], and we also apply it to CORK-GMRES. We derive conditions on the residual norms in both settings.

4.1 Inexact RK-GMRES↩︎

Inexact solves in LRK-GMRES transform the exact recurrence relations 4 into a form that includes the inner residuals, \[\label{eq:inexact-lrk-gmres-relation} \boldsymbol{\mathsf{A}}\boldsymbol{\mathsf{V}}_{m+1} \underline{\mathsf{H}}_m = \boldsymbol{\mathsf{B}}\boldsymbol{\mathsf{V}}_{m+1} \underline{\mathsf{K}}_m - \boldsymbol{\mathsf{P}}_m.\tag{19}\] Using \(\boldsymbol{\mathsf{W}}_m=\boldsymbol{\mathsf{V}}_{m+1}\underline{\mathsf{H}}_m\), with columns \(\boldsymbol{\mathsf{w}}_k\), \(\boldsymbol{\mathsf{P}}_m=[\boldsymbol{\mathsf{p}}_1 \, \cdots \, \boldsymbol{\mathsf{p}}_m]\) contains the residuals \[\boldsymbol{\mathsf{p}}_k := \boldsymbol{\mathsf{B}}\boldsymbol{\mathsf{v}}_{\mathsf{t},k} - (\boldsymbol{\mathsf{A}}-\sigma_k \boldsymbol{\mathsf{B}})\boldsymbol{\mathsf{w}}_k, \quad \boldsymbol{\mathsf{v}}_{\mathsf{t},k}:=\boldsymbol{\mathsf{V}}_k\mathsf{t}_k, \quad k=1,\dots,m .\] The final true residual of LRK-GMRES is composed of a contribution from \(\boldsymbol{\mathsf{P}}_m\) and the residual \(\tilde{\boldsymbol{\mathsf{r}}}^\sigma_m\) of the least-squares problem, which can be cheaply computed, \[\tilde{\boldsymbol{\mathsf{r}}}^\sigma_m:=\boldsymbol{\mathsf{V}}_{m+1}\left(\hat{\beta}\mathsf{e}_1-(\underline{\mathsf{K}}_m-\mu\underline{\mathsf{H}}_m)\mathsf{y}^\mathrm{red}_m\right).\] An appropriate choice for the inner tolerances can be established.

As was to be expected, there is a gap between the computed residual \(\tilde{\boldsymbol{\mathsf{r}}}^\sigma_m\) and the true final residual \(\boldsymbol{\mathsf{r}}_m\). [prop:inexact-lrk-gmres] does not take into account solving the first system with \(\sigma_0\) inexactly. It is clear that an inner residual \(\boldsymbol{\mathsf{p}}_0\) introduced when solving the system with \((\boldsymbol{\mathsf{A}}-\sigma_0\boldsymbol{\mathsf{B}})\), propagates into a contribution to \(\|\boldsymbol{\mathsf{r}}_m\|\) of at most \(\|\boldsymbol{\mathsf{p}}_0\|\).

The case of RRK-GMRES is similar, though somewhat simpler as the residual is not preconditioned. With inexact solves the recurrence relation 13 becomes \[(\boldsymbol{\mathsf{A}}-\mu \boldsymbol{\mathsf{B}})\boldsymbol{\mathsf{Z}}_m = \boldsymbol{\mathsf{V}}_{m+1}(\underline{\mathsf{K}}_m-\mu \underline{\mathsf{H}}_m) - \boldsymbol{\mathsf{P}}_m, \qquad \boldsymbol{\mathsf{p}}_k:=\boldsymbol{\mathsf{v}}_{\mathsf{t},k}-(\boldsymbol{\mathsf{A}}-\sigma_k \boldsymbol{\mathsf{B}})\boldsymbol{\mathsf{z}}_k.\] For a solution \(\boldsymbol{\mathsf{y}}_m(\mu)=\boldsymbol{\mathsf{Z}}_m\mathsf{y}^\mathrm{red}_m(\mu)\), the final true residual \(\boldsymbol{\mathsf{r}}_m:=\boldsymbol{\mathsf{r}}_m(\mu)\) becomes \[\boldsymbol{\mathsf{r}}_m = \boldsymbol{\mathsf{d}}- (\boldsymbol{\mathsf{A}}-\mu \boldsymbol{\mathsf{B}})\boldsymbol{\mathsf{Z}}_m\mathsf{y}^\mathrm{red}_m = \boldsymbol{\mathsf{V}}_{m+1}(\beta \mathsf{e}_1 - (\underline{\mathsf{K}}_m-\mu\underline{\mathsf{H}}_m)\mathsf{y}^\mathrm{red}_m) + \boldsymbol{\mathsf{P}}_m\mathsf{y}^\mathrm{red}_m = \tilde{\boldsymbol{\mathsf{r}}}_m + \boldsymbol{\mathsf{P}}_m\mathsf{y}^\mathrm{red}_m\] where \(\tilde{\boldsymbol{\mathsf{r}}}_m:=\tilde{\boldsymbol{\mathsf{r}}}_m(\mu)\) denotes the computed residual of the small least-squares problem.

The following is similar to the result derived in [23].

4.2 Inexact CORK-GMRES↩︎

When \((\boldsymbol{\mathsf{A}},\boldsymbol{\mathsf{B}})\) is the companion linearization of a matrix-valued function \(\mathsf{P}\), the residual \(\boldsymbol{\mathsf{p}}=\boldsymbol{\mathsf{y}}-(\boldsymbol{\mathsf{A}}-\mu \boldsymbol{\mathsf{B}})\boldsymbol{\mathsf{z}}\), \(\|\boldsymbol{\mathsf{p}}\|\le\tau\) originates from an inexact solve with \(\mathsf{P}(\mu)\). Under the assumption that all other operations in efficiently applying \((\boldsymbol{\mathsf{A}}-\mu \boldsymbol{\mathsf{B}})^{\smin 1}\) are performed exactly, the following holds.

In the CORK setting, emphasis must be placed on the fact that \(\boldsymbol{\mathsf{r}}_m\) is not actually the residual of our original problem. Instead, it is residual \(\mathsf{r}_m^\mathsf{P}:=\mathsf{b}-\mathsf{P}(\mu)\mathsf{x}_m\) that is of importance. Due to the expected form of \(\boldsymbol{\mathsf{y}}_m\), see 1, one generally takes \(\mathsf{x}_m:=\mathsf{y}_m^{[1]}/f_1(\mu)\). If the approximate solution is exactly representable as \(\boldsymbol{\mathsf{y}}_m = \mathsf{f}(\mu)\otimes \mathsf{x}_m\), then \(\|\mathsf{r}^\mathsf{P}_m\|=\|\boldsymbol{\mathsf{r}}_m\|\) (\(\mathsf{r}_m^{[1]} = \mathsf{r}_m^\mathsf{P},\; \mathsf{r}_m^{[2:d]}=\mathsf{0}\)), thus we expect \(\|\mathsf{r}^\mathsf{P}_m\|\approx\|\boldsymbol{\mathsf{r}}_m\|\).

When the right-hand side is parameter-dependent as discussed in 3.3, solving the augmented system can clearly be achieved through solving a system with \(\mathsf{P}(\mu)\) with some additional cheap operations. However, when considering the inexact solves, the previous conditions relate to the augmented system, not the original.

The final true residual on the augmented system can be written as follows, assuming that CORK-GMRES is applied with starting vector \(\mathsf{e}_{n+1}\), \[\begin{align} \underline{\mathsf{r}}_m^\mathsf{P}(\mu) := \dfrac{\underline{\mathsf{b}}(\mu)-\underline{\mathsf{P}}(\mu)\underline{\mathsf{x}}_m(\mu)}{\gamma(\mu)} &= \left[\begin{array}{c} \mathsf{0} \\ 1 \end{array}\right] - \left[\begin{array}{cc} \mathsf{P}(\mu) & -\mathsf{b}(\mu) \\ \mathsf{0}_{1\times n} & \gamma(\mu) \end{array}\right] \left[\begin{array}{c} \mathsf{x}_m/\gamma(\mu) \\ 1/\gamma(\mu) \end{array}\right] \\ &= \left[\begin{array}{c} (\mathsf{b}(\mu) - \mathsf{P}(\mu)\mathsf{x})/\gamma(\mu) \\ 0 \end{array}\right] = \left[\begin{array}{c} \mathsf{r}_m^\mathsf{P}(\mu)/\gamma(\mu) \\ 0 \end{array}\right]. \end{align}\] This shows that the true residual norm on the original system is scaled by \(1/\gamma(\mu)\) and thus the bounds on the inner residual norms should scale by \(1/|\gamma(\mu)|\). As the final element of the approximate solution is not exactly equal to \(1/\gamma(\mu)\) in general, we only expect an approximate equality of \(\|\underline{\mathsf{r}}_m^\mathsf{P}(\mu)\|\approx\|\mathsf{r}_m^\mathsf{P}(\mu)\|/|\gamma(\mu)|\).

Besides relating the final true residual \(\underline{\mathsf{r}}_m^\mathsf{P}\) of the augmented system to that of the original system, we must do the same for the conditions on the inner residuals of the inexact solves. Consider the residual \(\underline{\mathsf{q}}= \underline{\mathsf{w}}^{[1]} - \underline{\mathsf{P}}(\mu)\underline{\mathsf{z}}^{[1]}\) on the augmented system⁵, \[\underline{\mathsf{q}}= \left[\begin{array}{c} \mathsf{w} \\ w_{n+1} \end{array}\right] - \left[\begin{array}{cc} \mathsf{P}(\mu) & -\mathsf{b}(\mu) \\ \mathsf{0}_{1\times n} & \gamma(\mu) \end{array}\right] \left[\begin{array}{c} \mathsf{z}\\ z_{n+1} \end{array}\right] = \left[\begin{array}{c} \mathsf{w}+z_{n+1}\mathsf{b}(\mu)-\mathsf{P}(\mu)\mathsf{z}\\ w_{n+1}-\gamma(\mu) z_{n+1} \end{array}\right].\] It is clear that we can choose \(z_{n+1}=w_{n+1}/\gamma(\mu)\) such that \(\|\underline{\mathsf{q}}\|=\|\mathsf{q}\|\) where \(\mathsf{q}\) is the residual on the inexact system solve of \(\mathsf{P}(\mu)\mathsf{z}\approx \mathsf{w}+\frac{w_{n+1}}{\gamma(\mu)}\mathsf{b}(\mu)\). Note that the residual norms are equivalent but the right-hand sides are not, such that the required relative tolerance is different from the non-augmented case.

We conclude by summarizing the practical bounds. In LCORK-GMRES applied to \(\mathsf{P}\) with right-hand side \(\mathsf{b}\), the condition \[\|\mathsf{p}_k\| \le \min_{\mu\in S}\dfrac{|{\mathsf{e}}^\mathsf{T}_1\mathcal{P}\mathsf{f}(\sigma_k)|}{|\mu-\sigma_0|} \dfrac{\text{sval}_\text{min}(\underline{\mathsf{K}}_m-\mu\underline{\mathsf{H}}_m)}{m} \dfrac{\|\mathsf{b}\|}{\|\tilde{\boldsymbol{\mathsf{r}}}^{\sigma}_{k-1}(\mu)\|} \epsilon,\quad k=1,\dots,m,\] where \(\mathsf{p}_k\) is the inner residual of the inexact solve with \(\mathsf{P}(\sigma_k)\) at iteration \(k\), guarantees⁶ \(\|\mathsf{r}_m^\mathsf{P}(\mu)\|\le \|(\boldsymbol{\mathsf{A}}-\sigma_0 \boldsymbol{\mathsf{B}})\tilde{\boldsymbol{\mathsf{r}}}^{\sigma}_m(\mu)\| + \epsilon \|\mathsf{b}\|\) for all \(\mu\in S\). Similarly for RCORK-GMRES, \[\|\mathsf{p}_k\| \le \min_{\mu\in S}|{\mathsf{e}}^\mathsf{T}_1\mathcal{P}\mathsf{f}(\sigma_k)| \dfrac{\text{sval}_\text{min}(\underline{\mathsf{K}}_m-\mu\underline{\mathsf{H}}_m)}{m} \dfrac{\|\mathsf{b}\|}{\|\tilde{\boldsymbol{\mathsf{r}}}_{k-1}(\mu)\|} \epsilon,\quad k=1,\dots,m,\] guarantees \(\|\mathsf{r}_m^\mathsf{P}(\mu)\|\le\|\tilde{\boldsymbol{\mathsf{r}}}_m(\mu)\|+\epsilon \|\mathsf{b}\|\) for all \(\mu\in S\). For parameter-dependent right-hand sides, when the methods are applied to \(\underline{\mathsf{P}}\) with right-hand side \(\mathsf{e}_{n+1}\), the bounds should be divided by \(|\gamma(\mu)|\) and \(\|\mathsf{b}\|=\|\mathsf{b}(\mu)\|\) is also \(\mu\)-dependent.

4.3 Numerical experiment: parametrized Helmholtz↩︎

We illustrate the effectiveness of LCORK- and RCORK-GMRES, and the incorporation of inexact Krylov theory, by comparing to the inexact infinite GMRES method from [23]. We repeat the experiment from Figure 6.7(b) of [23], which originates from a finite element discretization of a 2D parametrized Helmholtz equation (see §6.2 for details⁷).

We perform CORK-GMRES for \(\mu\in[0,2.5]\) with shifts \(\boldsymbol{\sigma}=(0.25,1.25,2.25)\). We use both the regular versions with LU factorizations and the inexact versions with BiCGStab. In the latter, bounds in §4.2 are used for the inner tolerances, with \(\epsilon=10^{-9}\) and the unknown smallest singular values replaced by one. In LCORK-GMRES, \(\sigma_0\) is set to \(1.25\) and this initial inexact solve is performed using a relative tolerance of \(10^{-10}\). The sine function in the formulation is approximated with AAA [53], using a tolerance of \(10^{-13}\), which yields a rational approximant with \(8\) support points. The resulting linearization differs significantly from the Taylor series linearization in [23].

The results are shown in 1. We consider the maximum relative residual norm \(\|\mathsf{r}_k^\mathsf{P}(\mu)\|/\|\mathsf{b}\|=\|\mathsf{b}-\mathsf{P}(\mu)\mathsf{x}_k(\mu)\|/\|\mathsf{b}\|\) over \(\mu\in[0,2.5]\). Comparing to (inexact) infinite GMRES, the number of required iterations to achieve similar residual norms is significantly decreased, owing to the absence of initial stagnation. In the exact case, a relative residual norm of \(10^{-10}\) is achieved in \(20\) iterations rather than \(40\). The inner relative tolerances are shown to increase throughout the iterations. The iteration at which the inexact version stagnates between \(10^{-9}\) and \(10^{-10}\) aligns well with the inner tolerance reaching a value around 1.

5 Application to Helmholtz BEM↩︎

As main application of this manuscript, we consider the solution of scattering problems defined by the 3D Helmholtz equation, through the use of the Boundary Element Method (BEM).

5.1 Problem definition↩︎

The application under consideration is that of acoustic (exterior) scattering problems in 3D, solved efficiently using their integral equation formulation. The Helmholtz equation \(\Delta u + \kappa^2 u = 0\) is considered in an open domain \(\Omega\subseteq\mathbb{R}^3\) where total field \(u\) is composed of a given incident field \(u^I\) and an unknown scattered field \(u^S\), i.e.\(u:=u^I+u^S\). Boundary conditions are imposed on \(u\) and \(u^S\) must satisfy the Sommerfeld radiation condition \(\lim_{\mathbf{r}\to\infty}{\mathbf{r}\cdot\left(\partial u^S/\partial\mathbf{r} - i\kappa u^S\right)}=0\). The elliptic boundary value problem can be reformulated as a boundary integral equation (BIE) [54] with an integral operator of the form \[(\mathcal{A}_\kappa\phi)(\mathbf{x}) = \int_{\Gamma}g_\kappa(\mathbf{x},\mathbf{y})\phi(\mathbf{y})\,\mathrm{d}s_\mathbf{y}, \;\; \forall \mathbf{x}\in\Gamma:=\partial\Omega\] where the kernel function \(g_\kappa\) involves the 3D Helmholtz Green’s kernel \(G_\kappa(\mathbf{x},\mathbf{y}) = {\mathrm{exp}({i\kappa\|\mathbf{x}-\mathbf{y}\|})}/{4\pi\|\mathbf{x}-\mathbf{y}\|}\) and/or its derivative(s).

To solve the BIE, it is discretized. Galerkin BEM employs two function bases \(\{\psi_i\}_{i=1}^n\) and \(\{\phi_j\}_{j=1}^n\) and results in a dense matrix \(\mathsf{A}:=\mathsf{A}_\kappa\) with elements \[a_{ij} = \left\langle \mathcal{A}_\kappa\phi_j,\psi_i \right\rangle = \int_\Gamma{\int_\Gamma{g_\kappa(\mathbf{x},\mathbf{y})\psi_i(\mathbf{x})\phi_j(\mathbf{y})\mathrm{d}s_\mathbf{x}\mathrm{d}s_\mathbf{y}}}, \quad i,j=1,\dots, n.\]

5.2 System matrix representation↩︎

While dense, the matrix \(\mathsf{A}\) arising from Galerkin discretization exhibits data sparsity in the fact that it contains many large subblocks of approximate low rank. The \(\mathcal{H}\)-matrix format, introduced in [2], exploits this by representing these blocks explicitly by an approximate low-rank factorization. The \(\mathcal{H}\)-matrix approximation \(\mathsf{A}^{\mathcal{H}}\) of \(\mathsf{A}\) is characterized as \[{\mathsf{A}^{\mathcal{H}}}_{|b} := {\mathsf{A}}_{|b} \;\;\forall b\in\mathcal{L}^-, \quad\quad {\mathsf{A}^{\mathcal{H}}}_{|b} := \mathsf{X}_b\mathsf{Y}_b^\mathsf{H}\approx {\mathsf{A}}_{|b} \;\;\forall b\in\mathcal{L}^+,\] where \(b:=t\times s,\; t,s\subset \{1,\dots,n\}\), \(\mathsf{X}_b\in\mathbb{C}^{|t|\times k_b},\mathsf{Y}_b\in\mathbb{C}^{|s|\times k_b},\; k_b\le k=\mathcal{O}(1)\), and the set \(\mathcal{L}:=\mathcal{L}^+\cup \mathcal{L}^-\) tiles the whole matrix. This representation allows an approximate representation of \(\mathsf{A}\) in \(\mathcal{O}(n\log n)\) complexity instead of \(\mathcal{O}(n^2)\) [55]. The same complexity is achieved for the matrix-vector product.

In [16], a compact representation is presented for matrix-valued functions such as \(\mathsf{A}(\kappa):=\mathsf{A}_\kappa\), based on a ‘\(\mathcal{H}\)-tensor’ format and a ‘frequency extraction’ technique to analytically remove most of the wavenumber dependence. This technique is not compatible with linearization, hence we employ the representation without. The \(\mathcal{H}\)-matrix format of the BEM matrix for a given frequency is extended into a \(\mathcal{H}\)-tensor format, in which the continuous wavenumber dependence is approximated by low-degree rational functions. We keep the discussion on \(\mathcal{H}\)-tensors here brief, full technical details are provided in [34].

The \(\mathcal{H}\)-tensor approximation \(\mathsf{A}^\mathcal{H}(\kappa)\) of a matrix-valued function \(\mathsf{A}(\kappa)\) is defined using the same subdivision into subblocks as with \(\mathcal{H}\)-matrices, such that \[\label{eq:htensor-original} {\mathsf{A}^{\mathcal{H}}}_{|b}(\kappa) := \mathsf{U}_b\left(\sum_{\nu=1}^{R_b} \mathsf{G}_{b,\nu} \cdot \mathsf{h}_{b,\nu}^\mathsf{T}\mathsf{g}(\kappa) \right)\mathsf{V}_b^\mathsf{H}, \quad \forall b\in\mathcal{L}^+\tag{20}\] with \(\mathsf{U}_b\in\mathbb{C}^{|t|\times \ell_b},\mathsf{V}_b\in\mathbb{C}^{|s|\times r_b}\) orthogonal, \(\{\mathsf{G}_{b,\nu} \}_{\nu=1}^{R_b}\) possibly non-square matrices, and \(\mathsf{h}_{b,\nu}=[h_{b,\nu}(z_1) \,\cdots\, h_{b,\nu}(z_d)]^\mathsf{T}\) and \(\mathsf{g}(\mu)={[g_1(\mu) \,\cdots\, g_d(\mu)]}^\mathsf{T}\), where \[\label{eq:bary-basis-funs} g_j: \mathbb{C}\to \mathbb{C}: \mu \mapsto \left. \dfrac{w_j}{\mu-z_j} \middle/ \sum_{i=1}^d \dfrac{w_i}{\mu-z_i} \right. \quad \forall j\in\{1,\dots,d\}\tag{21}\] are the barycentric rational basis functions with support points \(\{z_j\}_{j=1}^d\) and weights \(\{w_j\}_{j=1}^d\). For \(b\in\mathcal{L}^-\), we have \({\mathsf{A}^{\mathcal{H}}}_{|b}(\kappa):=\sum_{\nu=1}^{R_b} \mathsf{A}_{b,\nu} \mathsf{h}_{b,\nu}^\mathsf{T}\mathsf{g}(\kappa)\).

Similarly to how blocks of \(\mathcal{H}\)-matrices can be constructed algebraicly in \(\mathcal{O}(n\log n)\) time through the use of Adaptive Cross Approximation (ACA)[56], these \(\mathcal{H}\)-tensor blocks can be efficiently constructed using tensorial ACA variants. After performing QR-AAA [57], an accelerated version of set-valued AAA [40] for simultaneous rational approximation of multiple scalar functions, we obtain the representation as given above, where all blocks share the \(d\) basis functions \(\mathsf{g}\).

5.3 Incorporating the excitation↩︎

The right-hand side of a scattering BIE stems from the incident field \(u^I\). However, generally, \(u^I\) depends on the wavenumber \(\kappa\), i.e.\(u^I(\mathbf{x}):=u^I(\mathbf{x},\kappa)\), just like the integral operator. When Galerkin BEM is employed, this results in a wavenumber-dependent right-hand side vector \(\mathsf{b}(\kappa)\).

To employ the augmented system technique (3.3), the right-hand side \(\mathsf{b}(\kappa)\) must satisfy condition 18 , or at least, an approximation to \(\mathsf{b}\) of this form must be constructed. This can be ensured by incorporating it into the input to QR-AAA. The end result is an approximation involving the same basis functions 21 , such that \[\label{eq:rhs-affine} \mathsf{b}(\kappa) \approx \tilde{\mathsf{b}}(\kappa) := \mathsf{B}^\mathsf{T}\mathsf{g}(\kappa),\quad \mathsf{B}\in \mathbb{C}^{d\times n}, \quad b_{ij}=b_j(z_i).\tag{22}\]

The addition of the right-hand side into the QR-AAA approximation may increase the degree of the final rational approximant, which impacts the size of the linearization and thus the cost of CORK. We will not examine the degree increase in this manuscript as we are specifically interested in the CORK-based frequency sweep, which requires including the right-hand side anyway.

5.4 Linearization and CORK application↩︎

With wavenumber-dependent representations of the system matrix 20 and right-hand side vector 22 , a linearization can be constructed and CORK-GMRES performed. The parametrized system matrix of the augmented system is defined as \[\underline{\mathsf{P}}(\kappa) = \left[\begin{array}{cc} \mathsf{A}^{\mathcal{H}}(\kappa) & -\tilde{\mathsf{b}}(\kappa) \\ \mathsf{0}_{1\times n} & \gamma(\kappa) \end{array}\right] = \sum_{i=1}^d \underline{\mathsf{C}}_i g_i(\kappa), \quad \underline{\mathsf{C}}_i = \left[\begin{array}{cc} \mathsf{C}_i & -\mathsf{b}_i \\ \mathsf{0}_{1\times n} & \gamma_i \end{array}\right]\] and for all \(i=1,\dots,d\), \(\mathsf{b}_i = \mathsf{b}(z_i)\), \[{\mathsf{C}_i}_{|b} = \mathsf{U}_b \left(\sum_{\nu=1}^{R_b}h_{b,\nu}(z_i) \mathsf{G}_{b,\nu}\right){\mathsf{V}}^\mathsf{H}_b \;\; \forall b\in\mathcal{L}^+, \quad {\mathsf{C}_i}_{|b} = \sum_{\nu=1}^{R_b} h_{b,\nu}(z_i) \mathsf{A}_{b,\nu} \;\; \forall b\in\mathcal{L}^-.\] To avoid scaling issues, the scalar functions \(\mathsf{f}\) in \(\ref{eq:matrix-function-add}\) are defined as \(f_1 := 1\) and \(f_i := g_{i-1}\), \(i=2,\dots,d\). One can map between \(\mathsf{g}\) and \(\mathsf{f}\) because \(\sum_{i=1}^d g_i(\kappa) \equiv 1\). The coefficient matrices in \(\ref{eq:matrix-function-add}\) are easily derived. We choose \(\gamma(\kappa) \equiv 1\) by defining \(\gamma_i:=1\) for all \(i=1,\dots,d\). The particular matrix pencil \((\mathsf{M},\mathsf{N}\)) such that \((\mathsf{M}-\kappa\mathsf{N})\mathsf{f}(\kappa)\), is found in [40], which is the final ingredient of the linearization \((\boldsymbol{\mathsf{A}},\boldsymbol{\mathsf{B}})\).

CORK relies on the capability to solve a system with \(\mathsf{A}^{\mathcal{H}}\) and perform a matrix-vector product with \(\bar{\boldsymbol{\mathsf{C}}}\). The system solve is performed inexactly through an iterative solver. For each shift we obtain a cheap matrix-vector product by ‘slicing’ the tensor: we construct a \(\mathcal{H}\)-matrix from the \(\mathcal{H}\)-tensor at the shift (see [16], [33], [34]). This is inexpensive compared to constructing the \(\mathcal{H}\)-matrix from scratch.

With the \(\mathcal{H}\)-matrix in hand, flexible GMRES (FGMRES) [58] is employed on the \(\mathcal{H}\)-tensor, where the flexible preconditioner is BiCGStab(2) [59], [60] applied to the \(\mathcal{H}\)-matrix. This is necessary as the recompression performed to obtain the \(\mathcal{H}\)-matrix results in a numerical difference between it and the \(\mathcal{H}\)-tensor.⁸ Using FGMRES allows us to use a fixed tolerance for both the \(\mathcal{H}\)-matrix compression as well as the inner BiCGStab(2) solution, throughout all CORK iterations. This lends itself nicely to the comparison with naively solving the systems at each wavenumber separately. In the numerical experiments, we show that FGMRES needs very few iterations to converge.

Performing the matrix-vector product \(\mathsf{y}=\bar{\boldsymbol{\mathsf{C}}}\boldsymbol{\mathsf{x}}\) for \(\mathsf{y}\in\mathbb{C}^n\) and \(\boldsymbol{\mathsf{x}}\in\mathbb{C}^{dn}\) is the second ingredient. We refer to this as coefficient matrix application (CMA). The operation can be defined in terms of matrix-vector products on the level of the blocks. The operation to perform for each block is \[{\mathsf{y}}_{|t}^{(b)} = \sum_{i=1}^d {\mathsf{C}_i}_{|b} {\mathsf{x}}_{|s}^{[i]}, \quad \forall b = t\times s\in \mathcal{L}\] after which all the computed subvectors are appropriately summed, similar to a normal \(\mathcal{H}\)-matrix-vector product. For admissible blocks, one can reorder the operation into \[{\mathsf{y}}_{|t}^{(b)} = \sum_{i=1}^d \mathsf{U}_b \left(\sum_{\nu=1}^{R_b}h_{b,\nu}(z_i) \mathsf{G}_{b,\nu}\right){\mathsf{V}}^\mathsf{H}_b {\mathsf{x}}_{|s}^{[i]} = \mathsf{U}_b \left[\begin{array}{ccc}\mathsf{G}_{b,1}&\cdots&\mathsf{G}_{b,R_b}\end{array}\right]\mathrm{vec}\left({\mathsf{V}}^\mathsf{H}_b{\mathsf{X}}_{|s}\mathsf{H}_b\right),\] \[{\mathsf{X}}_{|s} = \left[\begin{array}{ccc} {\mathsf{x}}_{|s}^{[1]} & \cdots & {\mathsf{x}}_{|s}^{[d]} \end{array}\right], \quad \mathsf{H}_b :=\left[\begin{array}{ccc} h_{b,1}(z_1) & \cdots & h_{b,R_b}(z_1) \\ \vdots & \ddots & \vdots \\ h_{b,1}(z_d) & \cdots & h_{b,R_b}(z_d) \end{array}\right] \quad\forall b\in\mathcal{L}^+.\] The most expensive part of this operation is the matrix-matrix product of \({\mathsf{X}}_{|s}\in\mathbb{C}^{|s|\times d}\) with \({\mathsf{V}}^\mathsf{H}_b\). In a \(\mathcal{H}\)-matrix-vector product with \(\mathsf{A}^\mathcal{H}\), one would have to perform a matrix-vector product with \({\mathsf{V}}^\mathsf{H}_b\) instead. The other expensive part is a matrix-vector product with \(\mathsf{U}_b\), also done in the regular \(\mathcal{H}\)-matrix-vector product. Thus, in terms of floating point operations, the total application of \(\bar{\boldsymbol{\mathsf{C}}}\) has a cost of approximately \((d+1)/2\) matrix-vector products with \(\mathsf{A}^\mathcal{H}\). However, the matrix-matrix products can be much more efficient on hardware, reducing the actual relative cost in practice.

6 Numerical experiments↩︎

We aim to evaluate the effectiveness of CORK-GMRES for frequency sweeping of exterior scattering problems via Helmholtz BEM. Using multiple examples, the convergence behavior over the frequency range is investigated for both the left- and right-preconditioned methods, and their cost is compared with naively solving systems at individual frequencies. Direct comparisons against other frequency sweeping approaches, as well as tackling high-frequency regimes, are out of the scope of this manuscript.

6.1 Experimental setup↩︎

In the experiments, we consider the exterior Dirichlet problem and use the indirect combined-field integral equation (CFIE) formulation [61]. The boundary integral equation involves the boundary integral operator \(\mathcal{A}(\kappa) + \frac{1}{2}\mathcal{I}\), with the identity operator \(\mathcal{I}\) and kernel function \[g_\kappa(\mathbf{x},\mathbf{y})=G_\kappa(\mathbf{x},\mathbf{y})-i\eta\dfrac{\partial G_\kappa(\mathbf{x},\mathbf{y})}{\partial n_\mathbf{y}}.\] We choose \(\eta=\kappa\) and construct the \(\mathcal{H}\)-tensor for \(\mathcal{A}\). The discretized identity results in a contribution of a sparse mass matrix \(\mathsf{M}\).⁹ The right-hand side \(\mathsf{b}\) contains the Dirichlet boundary values of the excitation. The excitation source is either a plane wave with direction \(\mathbf{d}\) (\(\|\mathbf{d}\|=1\)) or a point source located at \(\mathbf{p}\), s.t.for all \(i\in I\), \[b_{i}(\kappa) = -\int_\Gamma u^I(\mathbf{x};\kappa)\psi_i(\mathbf{x})\mathrm{d}s_\mathbf{x}\;\;\text{with}\;\; u^I(\mathbf{x};\kappa)=\mathrm{exp}(i\kappa\mathbf{x}\cdot\mathbf{d}) \;\;\text{or}\;\; u^I(\mathbf{x};\kappa)=G_\kappa(\mathbf{x},\mathbf{p}).\]

The discretization we employ uses linear discontinuous Galerkin BEM where a mesh of \(N_t\) triangles and \(N_v\) vertices approximates the boundary. The operators that are discretized have three basis functions for each triangle. Using a sparse matrix \(\mathsf{T}\), these are combined into appropriate continuous linear functions for each of the mesh vertices. This means that the \(\mathcal{H}\)-tensors and \(\mathcal{H}\)-matrices we construct are of size \(3N_t \times 3N_t\) while the actual system matrix is of size \(N_v\times N_v\). Applying it involves only additional cheap operations with \(\mathsf{T}\). See [33] for more details.

For all the experiments, we apply the inexact greedy versions of LCORK- and RCORK-GMRES. For the inexactness, we use the bounds from 4.2 with \(\epsilon=10^{-5}\) and unknown smallest singular value replaced by \(1\). The relative tolerance used for \(\sigma_0\) in LCORK-GMRES is set to \(\epsilon/10\). \(\sigma_1\) and \(\sigma_2\) are fixed at the endpoints of the considered frequency range, and \(\sigma_0\) is put at the midpoint.

The system solution and coefficient matrix application are performed as discussed in 5.4. In BiCGStab(2), we use \(\mathsf{M}^{\smin 1}\) as right preconditioner. Unless specified otherwise, the same construction tolerance \(\tau=10^{-5}\) is used for \(\mathcal{H}\)-tensor and direct \(\mathcal{H}\)-matrix construction, as well as the construction of \(\mathcal{H}\)-matrices through the \(\mathcal{H}\)-tensor. To compute the integrals, Sauter quadrature [62] of degrees \((7,11,13,15)\) are employed. For both AAA approximation and greedy shift selection, we employ 10 samples per dimensionless wavenumber (\(\kappa\cdot\mathrm{diam}(\Gamma)\)).

All experiments are performed on a compute node consisting of 2 Intel Xeon Platinum 8360Y (Ice Lake) CPUs with 36 cores each, one L3 cache per CPU and 256 GiB total RAM.¹⁰ Julia code is used for the CORK framework with open-source C package BEACHpack providing efficient implementation of the BEM functionality covered in 5.¹¹

6.2 Performance of greedy CORK-GMRES↩︎

To assess the effectiveness of CORK-GMRES using the greedy shift selection, a first experiment is done on a trefoil knot where the incident field originates from a plane wave going through the hole of the knot. LCORK-GMRES and RCORK-GMRES is applied to the frequency range \(\kappa\cdot\mathrm{diam}(\Gamma)\in[10,30]\) for \(30\) iterations.

2 shows the maximum residual norm over the frequency range, together with the residual norm perceived by the greedy selection (i.e.the residual norm at the wavenumber where the estimate is largest). Both methods show good convergence, attaining a relative residual norm below \(10^{-5}\) at iteration \(22\) (LCORK) and \(23\) (RCORK). In addition, the convergence perceived through the greedy shift selection corresponds well to the actual convergence. 2 also reports the cumulative number of BiCGStab(2) solves, showing that \(22\) (\(23\)) iterations require \(38\) (\(41\)) system solutions for LCORK (RCORK). Recall that the number of system solutions can be seen as a proxy for the cost of CORK-GMRES in terms of naive system solutions at individual wavenumbers, although this assumes the system solves encompass the majority of the runtime.

In 3, the residual norms are plotted as a function of the frequency, for iteration numbers \(10\), \(20\) and \(30\) and a stricter inexactness tolerance of \(10^{-10}\). The greedy selection is effective in reducing the residual norm across the range.¹²

6.3 Convergence results on multiple meshes↩︎

The robustness of CORK-GMRES is tested by performing 30 iterations on six different shapes with \(\kappa\cdot\mathrm{diam}(\Gamma)\in[20,30]\). The number of triangles and vertices range from 98304 to 195014, and 49154 to 97509 respectively. ¿tbl:tab:meshes? reports on the required CORK-GMRES iterations (Out), the number of required solves, i.e.FGMRES iterations (Sol), and the total number of matrix-vector products (MV), to attain a maximal relative residual norm below \(\epsilon=10^{-1}\), \(10^{-3}\) and \(10^{-5}\). The time required to attain \(\epsilon=10^{-5}\) is reported as well.

The results of LCORK- and RCORK-GMRES are comparable, both in terms of BiCGStab(2) solves and matrix-vector products, as well as time spent. The required number of ‘naive solves’ to attain the desired accuracy over the whole wavenumber interval ranges from \(17\) to \(33\). Note that the inexact Krylov bounds in 4 are partly heuristic, and the smallest singular value is substituted by \(1\). We have observed that some cases require to tighten the inexactness tolerance. Scaling also influences where convergence stagnates, see Remark [rmk:scaling].

6.4 Computational and memory cost of the \(\mathcal{H}\)-tensor↩︎

The computational and memory cost in using CORK-GMRES to perform frequency sweeping, is largely determined by the \(\mathcal{H}\)-tensor. \(\mathcal{H}\)-tensors are constructed for the indented cube, hit by a point source, at three construction tolerances and over four different frequency ranges. ¿tbl:tab:abs-htensor-cost? compiles the memory footprint, degree and construction time of the \(\mathcal{H}\)-tensors, as well as the three important \(\mathcal{H}\)-tensor operations: slicing, matrix-vector multiplication (MV) and coefficient matrix application (CMA). More extensive experiments on the \(\mathcal{H}\)-tensor representation are performed in [34].

Both widening the wavenumber interval, and moving it to higher frequencies, increase the cost in all aspects, although it is mainly the former that has a significant impact here. Lowering the tolerance also has a large effect, especially on the construction time, because the quadrature degrees are increased as well.¹³ Notice that the computation time of CMA compared to MV, is lower than the degree would indicate.

In ¿tbl:tab:rel-htensor-cost?, the memory, construction and matrix-vector product of the \(\mathcal{H}\)-tensor are compared to those of the \(\mathcal{H}\)-matrices at the frequency boundaries. Memory footprint of the \(\mathcal{H}\)-tensor is around 2.7–3.5 times greater than a \(\mathcal{H}\)-matrix. The matrix-vector product shows a similar ratio, as expected. The increase in construction time is more significant, ranging between 13 and 24 times longer. This discrepancy is partly because tensorial ACA first overestimates the ranks of the \(\mathcal{H}\)-tensor, which is then alleviated by a recompression step. Additionally, whereas memory usage depends only linearly on the rank, the construction also has terms with quadratic and cubic rank dependence. Thus, the rank increase from \(\mathcal{H}\)-matrix to \(\mathcal{H}\)-tensor impacts the construction time more strongly.

Using Tables ¿tbl:tab:meshes?, ¿tbl:tab:abs-htensor-cost? and ¿tbl:tab:rel-htensor-cost?, one is able to more concretely compare the total cost of CORK-GMRES against a naive solve. Let us investigate the frequency sweeping of the indented cube using RCORK, over \(\kappa\cdot\mathrm{diam}(\Gamma)\in[20,30]\) with tolerances of \(10^{-5}\). Considering \(19\) outer iterations, \(33\) solves, and \(897\) matrix-vector products, the total cost, neglecting orthogonalization, is: \(19\) CMAs and slices, \(33\) MVs with the \(\mathcal{H}\)-tensor, and \(897\) MVs with \(\mathcal{H}\)-matrices. This results in a total time of \(137.3\) seconds, which corresponds well with the total time measured in ¿tbl:tab:meshes?. A naive solve, assuming \(27\) matrix-vector products are required (In/Sol), would take around \(2.94\) seconds. Again, actual measured naive solves take \(2.92\)–\(3.50\) seconds. The estimated number of desired frequencies such that CORK-GMRES is superior, is only around \(40\) to \(53\).

6.5 Frequency sweeping acoustic scattering by a submarine↩︎

As a final experiment, we perform a frequency sweep on the larger mesh of a submarine (\(N_t=\num{400886}\), \(N_v=\num{200445}\)).¹⁴ The source is a plane wave with frequency between \(100\) and \(250\) Hz. This corresponds to \(\kappa\cdot\mathrm{diam}(\Gamma)\in[25.97,64.93]\) or about 4 to 10 wavelengths over the submarine’s length. During the greedy shift selection, the true residual norm at the next shift is recorded. Once this norm is below the tolerance three times in a row, we consider CORK-GMRES to have converged, which was the case here.

4 displays the time spent in the different parts of CORK-GMRES. The orthogonalization cost is negligible and thus not shown. The main cost is in the system solution, around 80% of the total time. The cost of LCORK- and RCORK-GMRES is similar. To improve greedy shift selection and detection of convergence, one may compute the true residual norm at multiple wavenumbers. Say we test \(\ell\) wavenumbers, the cost of greedy shift selection would be roughly \(\ell\) times higher.

Finally, we compare LCORK-GMRES to naively solving at a set of wavenumber samples in the range. 5 plots the cost in time as a function of the number of samples, including construction costs. For LCORK-GMRES, the time also considers the computation of the approximate solution \(\mathsf{x}_m(\kappa)\) for each sample. There are two types of naive solves: constructing a \(\mathcal{H}\)-matrix at each sample directly, or constructing one \(\mathcal{H}\)-tensor and slicing it at each sample.

First, we mention that on average, the construction of the \(\mathcal{H}\)-tensor took approximately \(52.2\) times as long as the \(\mathcal{H}\)-matrix construction at the middle frequency while the ratio in memory footprint is only around \(3.35\). The discrepancy is larger than in ¿tbl:tab:rel-htensor-cost?, which is explained by worse load balancing. It is observed that for the \(\mathcal{H}\)-tensor, construction of the largest blocks takes longer, relative to the small blocks, than for the \(\mathcal{H}\)-matrices. If a few blocks, less than the number of cores, take most of the time, load balancing degrades unless parallelization is done within those blocks. Further investigation, e.g.comparing to sequential construction, is necessary.

The ratio in construction time influences the results in 5. Using the \(\mathcal{H}\)-tensor, either with CORK-GMRES or without, pays off once over 50 to 55 samples. If CORK-GMRES is used, almost all the cost is ‘offline’, because computing the approximate solutions is very cheap. Thus, it remains feasible to evaluate the solution at an increasingly large number of frequencies, in contrast to using only the \(\mathcal{H}\)-tensor.

Finally, we note that, in this experiment, the additional cost of CORK-GMRES is low relative to the \(\mathcal{H}\)-tensor construction, partially due to suboptimal parallelization. However, other parameters also affect this ratio, e.g., lowering the quadrature degrees reduces construction cost. Therefore, the relative cost between \(\mathcal{H}\)-tensor construction and CORK-GMRES execution may vary across setups. Frequency sweeping with the combination stays efficient regardless, given that the cost of repeated system solution depends on the parameter setup in a similar fashion.

7 Conclusions and future work↩︎

We have developed a new left- and right-preconditioned Krylov method for solving parametrized linear systems through linearization. While this results in very large systems, using the CORK framework the structure within the companion pencil is exploited and efficiency is maintained. The methods are more general than existing linearization methods and Krylov methods for shifted systems. Moreover, convergence is improved through the use of multiple shifts, which can be chosen greedily and fully adaptively in a cheap manner. Incorporating inexact Krylov theory reduces computational cost, since the tolerances in later iterations can be loosened while maintaining accuracy of the overall solution. Finally, we adapt the methodology to allow for parameter dependence in the right-hand side.

In combination with the \(\mathcal{H}\)-tensor format, the framework enables efficient frequency sweeping of three-dimensional scattering problems through Helmholtz BEM. It performs well for many shapes and moderate frequencies, and already outperforms the naive approach computationally once a modest number of frequency samples are desired, albeit with increased memory usage.

Future work may include exploration of broader applications, model order reduction (e.g.including field evaluation in Helmholtz BEM), development of compatible data-sparse representations for high-frequency scattering problems (for example inspired by [15]), the computation of scattering poles (i.e.eigenvalues of integral operator \(\mathcal{A}(\kappa)\)), and further investigation of the effects of scaling in linearization.

References↩︎

---

1. Department of Computer Science, KU Leuven, Leuven, 3000 (, , ).↩︎

2. To be submitted to the journal’s Software, High-Performance Computing, and Computational Methods in Science and Engineering section 2026.↩︎

3. One may consider ‘rational Krylov subspace’ an incorrect term, as it is not a Krylov subspace. Nevertheless we will use it as shorthand for the subspace spanned by a ‘rational Krylov sequence’.↩︎

4. We use the convention of defining shifted systems by subtracting the second matrix instead of adding. This slightly changes the equations compared to [30].↩︎

5. Note that we use similar notation as in [prop:inexact-cork].↩︎

6. Remember that these “guarantees" hinge on the assumption that the full approximate solution \(\boldsymbol{\mathsf{y}}_m\) is representable as \(\mathsf{f}(\mu)\otimes\mathsf{x}_m\) for some \(\mathsf{x}_m\). These assumptions are practically never satisfied exactly, only approximately.↩︎

7. Implementation note: to obtain the large system, we have set ‘200’ to ‘400’ on line 103 of fem_matrices.py at github.com/siobhanie/FlexInfGMRES.↩︎

8. Solving a system with the \(\mathcal{H}\)-matrix directly is possible as well, as one can expand the inexact Krylov theory to also take into the account the error on the system matrix.↩︎

9. The mass matrix is appropriately incorporated in all the operations we perform with the system matrix. Its computational contribution is negligible.↩︎

10. The compute node is part of the ‘wICE’ cluster provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation – Flanders (FWO) and the Flemish Government.↩︎

11. The version of BEACHpack used to produce the results of the numerical experiments can be found at gitlab.kuleuven.be/numa/software/beachpack/-/tree/kb-exp-paper-corkgmres. The Julia code is available upon request from the authors.↩︎

12. The scale of the computed residual norm differs from that of the true one, yet it is the shape over the frequency range that matters. The discrepancy in scale is because the computed residual is that of the pencil, not the actual parametrized system.↩︎

13. For tolerances \(\tau=10^{-3}\) (\(10^{-7}\)), quadrature degrees are decreased (increased) by \(2\).↩︎

14. Here, we used a VSC ‘big memory’ node to fit the \(\mathcal{H}\)-tensor in memory. The hardware difference is the RAM increase to 2048 GiB and that two NUMA domains are present instead of one.↩︎