{'fraction_non_alphanumeric': 0.02965693600345389, 'fraction_numerical': 0.010660555943010859, 'mean_word_length': 4.825498162120333, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Extended Krylov Subspace Model-Order Reduction Technique to Simulate Wave Propagation in Unbounded Domains

In this paper we present a novel extended Krylov subspace reduced-order modeling technique to efficiently simulate time- and frequency-domain wavefields in open complex structures. To simulate the extension to infinity, we use an optimal complex-scaling method which is equivalent to an optimized perfectly matched layer in which the frequency is fixed. Wavefields propagating in strongly inhomogeneous open domains can now be modeled as a non-entire function of the complex-scaled wave operator. Since this function contains a square root singularity, we apply an extended Krylov subspace technique to construct fast converging reduced-order models. Specifically, we use a modified version of the extended Krylov subspace algorithm as proposed by Jagels and Reichel [Linear Algebra Appl., Vol. 434, pp. 1716 - 1732, 2011], since this algorithm allows us to balance the computational costs associated with computing powers of the wave operator and its inverse. Numerical experiments from electromagnetics and acoustics illustrate the performance of the method.


Introduction
Frequency-and time-domain wave problems in unbounded domains emerge in a wide range of applications: from seismic processing in geophysics to metamaterial development, from bioelectromagnetics to radar and antenna design. Furthermore, developing fast and robust forward modeling methods is not only a significant topic by itself, it is also of great importance in inversion and optimization algorithms. In all these different disciplines, conventional forward solution techniques are usually the method of choice. For time-domain problems, it is typically the Finite-Difference Time-Domain method (FDTD method, see, for example, [1]) that is used, while discretized frequency-domain problems are usually solved using linear system solvers. When solutions for multiple frequencies are required, conventional workflow consists of solving frequency-domain problems for each frequency separately. In many applications, such methods quite often leave much room for improvement. Indeed, discretization grids for 3D problems may consist of up to a billion nodes and, even with state-of-the-art preconditioners, obtaining the solution of discretized frequency-domain problems results in rather computationally intensive tasks. Projection-based model reduction is a well-established tool and allows us to obtain time-and frequency-domain solutions by projecting the large scale dynamic system on a small Krylov or rational Krylov subspace. In an FDTD setting, for example, optimal time steps can be determined, though rather implicitly, via polynomial Krylov subspace projection methods. These Krylov methods have been shown to provide significant speed up for lossy diffusion-dominated problems, e.g., for the diffusive Maxwell equations (see [18]). For lossless problems, however, e.g. wave problems in lossless media on bounded domains, such an optimality hardly provides any acceleration compared with time-stepping at the Courant limit [4].
The principal difference between wave problems in bounded and unbounded domains is that the latter one is essentially lossy. Indeed, loosely speaking, infinity acts like an absorber and for any bounded subdomain there is no wave permanently traveling inside it without energy loss. Therefore, an attempt to construct a Krylov subspace projection method to speed up the solution of such problems is not hopeless. In fact, two of the authors showed that for wavefield problems on unbounded domains, standard polynomial Krylov methods may indeed outperform FDTD on large time intervals provided we respect the delicate spectral properties of the corresponding continuous problem (see [2]). To be more precise, in any practical implementation the domain of computation is necessarily bounded and boundary conditions have to be imposed for domain truncation. If selfadjoint boundary conditions (Dirichlet boundary conditions, for example) are used then the solution to a selfadjoint problem with a continuous spectrum is approximated by a solution to a selfadjoint problem with a discrete spectrum. In the latter case, the solution consists of an infinite summation of periodic eigenfunctions that do not decay. Enlarging the domain of computation does not help, since the eigenfunctions remain periodic and do not show any attenuation. For general sources then, polynomial Krylov methods used on discretized bounded domains with selfadjoint boundary conditions cannot outperform FDTD running at the Courant limit. However, for non-selfadjoint boundary conditions that simulate the extension to infinity, polynomial Krylov methods may outperform FDTD on large time intervals. This can be explained by considering the scattering poles expansion of a wavefield u satisfying the wave equation on an unbounded domain [5]. In particular, Lax and Phillips showed that for any fixed point in R n , n odd, the solution u(t) can be expanded for large t in a sum of time-exponential modes. These modes decay exponentially in time and correspond to the so-called scattering resonances. The resonances are the poles of the resolvent for the wave equation and are located on the second Riemann sheet of the square root function. In practice, only a finite number of the corresponding modes essentially contribute to the solution for large t. Consequently, polynomial Krylov methods may outperform FDTD for large times if we can capture these resonant poles via a Krylov subspace method.
The scattering poles can be found by truncating the domain of interest via complex coordinate scaling also known as Aguilar-Balslev-Combes-Simon theory [6]. This coordinate stretching technique was introduced in the 1970s and is since then used extensively in atomic and molecular physics to calculate energies, widths, and cross-sections of open quantum systems. The technique coincides with Berenger's perfectly matched layer (PML, [7]) using a fixed frequency for coordinate stretching. The problem with complex scal-ing is, however, that resonances and anti-resonances of the system are found simultaneously resulting in instability of the traditionally computed solution via a sin or cos function of the operator. Two of the authors resolved this issue in [2] by introducing a so-called stability-corrected exponential function that depends on the square root of the wave operator. Any type of Krylov subspace method can now be applied and we split them into two groups depending on the application. The first group targets to applications where the offline cost of constructing the subspace is more important compared to the online cost of subsequent projection. For example, as a stand-alone forward solver, the latter cost is typically negligible. Polynomial Krylov subspace methods can be quite competitive here in terms of computation time, especially when the frequency range doesn't include small frequencies. However, when applied within an inversion/optimization framework, not only the forward problem solution is required but also the derivatives with respect to the parameters that are being optimized or inverted. Depending on the number of these parameters, the computation of the derivatives can make the projection cost dominating. For these applications, multiple solutions of the problem are required, and the same reduced order model can be reused for this purpose. In this case, it is beneficial to construct a subspace of smaller size for a given accuracy even at a larger offline cost. We put the methods of this kind into the second group. Taking this into account, we note that for non-entire operator functions (such as the stability-corrected operator function considered here), polynomial Krylov subspace methods are known to be not the best choice in terms of convergence speed. Indeed, rational Krylov subspace (RKS) methods are known to be much better suited for non-entire functions (see, e. g., [11]). Moreover, it has recently been shown that RKSs are beneficial for the solution of regular enough wave problems on bounded domains [9].
In this paper, we follow the approach based on the stability-corrected exponential function and construct optimal rational field approximants via a so-called Extended Krylov Subspace Method (EKSM). This Krylov method generates field approximations using both positive and negative matrix powers and produces rational field approximants with restricted poles at the origin and at infinity. The original EKSM of [11] allowed for short recursions and a structured projection matrix only for a fixed number of positive or negative matrix powers. More flexible EKSM modifications, allowing for sequences of nested extended subspaces with increasing positive and negative matrix powers, were suggested in [12,13]. Here we use a modified version of the algorithm proposed in [13], since it allows for an arbitrary ratio between positive and negative matrix powers. We illustrate the performance of our proposed solution technique in Section 5, where we present a number of numerical experiments related to the solution of the wave equation and solutions of Helmholtz's equation for a wide range of frequencies. In our experiments for time-domain problems, we show that our solution technique significantly outperforms polynomial Krylov subspace projection in terms of subspace dimension, while for Helmholtz's equation we show that with the PKS method we obtain solutions for a complete frequency range at the cost of an unpreconditioned linear solver applied to the problem with the smallest frequency. Our results clearly show that EKS significantly outperforms a polynomial Krylov approach when solutions for small frequencies are required.

Problem formulation and polynomial Krylov subspace projection
Let us consider the hyperbolic initial value problem for t ≥ 0. Its frequency-domain counterpart is In these equations, A is a selfadjoint nonnegative partial differential equation operator on an unbounded domain having an absolutely continuous spectrum. The solution of Eq. (1) can be written in terms of operator functions as where η(t) is the Heaviside unit step function. In [2,3] we discretized the above partial differential equation on a uniform second-order grid and used Zolotarev's optimal rational approximation of the square root to implement an optimal complex scaling method for domain truncation. (The optimal scaling approach of [2] deals with propagative modes only, while in [3] a more advanced implementation is presented in which evanescent and propagative modes are taken into account simultaneously.) Writing the resulting system asÃ where the tilde indicates that complex-stretching has been applied and the subscript N signifies that we are dealing with a discretized system, it is tempting to take the solution as an approximation to the exact solution as given by Eq. (3). The approximation of Eq. (5) is unstable, however, since the eigenvalues ofÃ N approximate both the resonances and the anti-resonances of the system. Therefore sin(Ã 1/2 N t) increases exponentially as t increases and u N is unstable. Similarly, important frequency-domain field properties are lost if we take as an approximation to u(s). In particular, the symmetry property u(s) = u(s) is broken and the limit lim s→0 u N (s) is not purely real as it should be, of course. Fortunately, we can correct for these failures by exploiting the symmetry of the exact solution with respect to (−Ã N ) 1/2 and (−Ã N ) 1/2 , where the overbar denotes complex conjugation. Specifically, setting and taking the principle value for the square root, a stabilized approximation to the continuous open problem of Eq. (1) is given by [2] This approximation is unconditionally stable because of the principle value convention. We note that if we replaceÃ N by A in Eq. (8), then the expression for the stability-corrected field becomes identical to the operator function expression for u in Eq. (3). In addition, by applying a Laplace transform to Eq. (8), we obtain the frequency-domain result where I N is the identity matrix of order N. This stability-corrected field satisfies the symmetry propertyũ N (s) =ū N (s) and becomes purely real as s → 0. Finally, we note that for frequency-domain wave fields operating at an angular frequency ω we take s = iω, of course.
Having the stability-corrected field approximations available, we can now construct polynomial reduced-order models in the usual way. In particular, the models are drawn from the Krylov subspace and are obtained by following a Petrov-Galerkin procedure. We remark that a basis for K m can be generated very efficiently via a Lanczos-type algorithm using a short three-term recurrence relation, since matrixÃ N is essentially complex-symmetric (see Section 4). Furthermore, the error of the reducedorder models based on the Krylov space of Eq. (10) is of the same order as the error of the PML discretization inÃ N on the frequency interval of the signal. Finally, we note that matrixÃ N is only required to form matrixvector products in this algorithm and the cost per iteration is close to the cost of a single FDTD iteration. Further details can be found in [2], where it is also shown that the reduced-order models generated in this way eventually do indeed outperform FDTD on large time intervals.

Extended Krylov Subspace Reduced-Order Models
Given that the stability-corrected field approximations of Eqs. (8) and (9) have a square root singularity, we observe that the reduced-order models taken from the Krylov subspace of Eq. (10) may not provide us with the fastest convergence. The reason is that models taken from this space are polynomial field approximations, but rational functions are usually much better at approximating functions near singularities (see [8]). This then naturally leads us to consider rational Krylov subspaces for stability-corrected field approximations. However, the main drawback of utilizing rational Krylov spaces is that large shifted systems of equations (Helmholtz systems) need to be solved to generate bases for these spaces. Without any effective preconditioners, this may simply be unfeasible.
To handle the square root singularity, we therefore resort to the extended Krylov subspaces An extended Krylov subspace can be seen as a special case of a rational Krylov subspace with one expansion point at zero and one at infinity. Obviously, the action ofÃ −1 N on a vector is required to generate a basis for such a space. This essentially amounts to solving a Poisson-type equation for which efficient solution techniques are available [10]. Computing matrixvector products with the inverse of the system matrix is generally still more expensive than computing matrix-vector products with the system matrix itself, however, and from a computational point of view we therefore prefer to deal with extended Krylov subspaces K m 1 ,m 2 with m 1 < m 2 . In addition, the basis vectors should be constructed via short term recurrence relations, since storage of all basis vectors is generally not practical for large-scale real-world problems.
In the original EKS method as introduced for Hermitian matrices in [11], the number of inverse powers m 1 is fixed and nested subspaces are recursively generated via short term recurrence relations. Furthermore, it is shown that for the matrix square root and m 2 = m 1 , the number of iterations required for convergence of the EKS reduced-order models is of . We stress that this estimate holds for Hermitian matrices A N and diagonal (m 1 = m 2 ) subspaces only. Generating a basis for such spaces obviously requires an equal number of iterations with the system matrix and its inverse.
In the EKS method proposed in [12], again a diagonal extended Krylov subspace is considered, but this time the sequence of nested subspaces is generated via short recurrence relations. The method converges at a rate as predicted by theory, but its computational costs are very large.
In [13,14], Jagels and Reichel developed a more generalized approach in which the sequence of subspaces is generated again via short term recurrences. Here, i is an integer that allows us to optimize the computational costs. In particular, since matrixvector products with the system matrix are usually more easily computed than matrix-vector products with its inverse, it may be advantageous to apply the extended Krylov approach of Jagels and Reichel, since the number of iterations with the system matrix is essentially i times larger than the number of iterations k that are carried out with the inverse of the system matrix. Different values of i do produce different extended Krylov subspace sequences, however, and the reduced-order models based on these spaces may have different approximating properties. Furthermore, the Krylov projection matrix is a pentadiagonal matrix of order k(i + 1) with a block-tridiagonal structure and consequently the EKS reduced-order models can be evaluated very efficiently.
In this paper, we use a variant of the extended Krylov scheme as proposed by Jagels and Reichel to construct EKS reduced-order models for wavefield problems in open domains. We exploit the complex symmetric structure of the system matrix to generate a basis of the extended Krylov space K k,ki+1 . The basis vectors are complex orthogonal by construction and can be computed via short term recurrence relations. This scheme is presented in Section 4 and provides us with an extended Krylov decomposition of the system matrixÃ N . The computed decomposition then allows us to construct the EKS reduced-order models for the stability-corrected field approximation of Eqs. (8) and (9).

Extended Krylov Subspace Algorithm
As mentioned above, we are interested in constructing EKS reduced-order models for wavefield quantities in unbounded domains that satisfy the wave equation with variable coefficients. To be specific, let us consider wavefields propagating in two-or three-dimensional unbounded domains. In the interior, we discretize the wave equation on an equidistant second-order grid and use optimal Zolotarev complex-scaling to simulate the extension to infinity. We do not provide any details about this discretization process here, since it is completely described in [2,3]. We do mention, however, that our approach is not limited to second-order schemes in the interior. Other higher-order discretization schemes can be used as well. In fact, since the optimal scaling-method has spectral accuracy, it would be preferable to use a discretization scheme that has spectral accuracy in the interior as well. Highorder spectral methods can be applied for example (see [15]) or optimal grids for interior domains as described in [16]. Since our focus is on constructing EKS reduced-order models, we use a straightforward second-order accurate scheme for simplicity only.
After the discretization procedure, we arrive at the semi-discretized system of Eq. (4). Due to complex-scaling, the system matrixÃ N belongs to C N ×N and is not selfadjoint (Hermitian). It does, however, satisfy the symmetry relationÃ where M is a diagonal matrix with products of the spatial stepsizes and medium parameters on its diagonal (for details, see [2]). Note that matrix M is not positive definite and does not even belong to R N ×N , since the optimal stepsizes used for complex-scaling are complex. Introducing now the bilinear form we have Ã N x, y = x,Ã N y for any x, y ∈ C N and two (nonzero) vectors x and y belonging to C N are said to be complex orthogonal if x, y = 0. Finally, we stress that Eq. (13) does not define an inner product, since M is not positive definite. Given the symmetry of matrixÃ N , we can now generate a complex orthogonal basis for K k,ki+1 by modifying the orthogonalization algorithm for K k,ki+1 and real symmetric matrices of Jagels and Reichel [13,14]. In particular, if we replace the inner products for vectors from R N by the bilinear form of Eq. (13), we arrive at an algorithm that generates a complex orthogonal basis of the Krylov space K k,ki+1 (Ã N , b N ), where each basis vector has a Euclidean norm equal to one. We note that the algorithm may break down, since the bilinear form does not define a definite inner product. With · denoting the Euclidean norm, the algorithm is as follows: end After a successful completion of this algorithm, we store the vectors of the generated basis in the N-by-k(i + 1) matrix and all iterations can be summarized into the equatioñ where and h ij is the (i, j) entry of matrix H k(i+1) . Furthermore, matrix H k(i+1) is a pentadiagonal matrix (see [13]) that satisfies where D k(i+1) is diagonal matrix with entries δ 0 , δ 1 , ..., δ i , δ −1 , ..., δ −k+1 , ..., δ ik . The entries of matrix H k(i+1) can be easily obtained in explicit form from the ones given in [13]. Equation (15) is the extended Krylov decomposition of the system ma-trixÃ N and can be used to construct the EKS reduced-order models both in the time-and frequency-domain. Specifically, setting d = k(i + 1) and introducing where again the principle value for the square root is taken, the time-domain EKS reduced-order model is given bỹ where e 1 is the first column of the identity matrix of order d. Its frequencydomain counterpart is given bỹ

Numerical Experiments
To illustrate the performance of the EKS method in the time-and frequencydomain, two sets of experiments are presented. In the first set an acoustic wavefield problem in the frequency-domain is considered, while in the second set we consider an electromagnetic photonic waveguide problem in the time-domain.
For the acoustic problem, we have considered a 3D SEG/Salt model and benchmarked our solver against an independently developed unpreconditioned BiCGStab algorithm that is a competitive conventional iterative Helmholtz solver (see [17] for details). Here we took the case m 1 = 1 corresponding to PKS. Figure 1 shows the solutions at different receiver locations for a fixed frequency of 7.5Hz. Both the BiCGStab and PKS model reduction solutions coincide with each other very well. When comparing the number of iterations, we note that the PKS approach requires one matrixvector multiplication per iteration, while BiCGStab needs two. Since this part constitutes the most computationally intensive part, it makes sense to compare the performance of the two methods in terms of matrix-vector multiplications rather than in terms of iterations. To this end, we have plotted the convergence rates of BiCGStab and model order reduction against the number of matrix-vector multiplications in Figure 2 (left) again for a frequency of 7.5Hz. We observe that for this single frequency problem both rates are rather close. We note that both approaches converge faster for higher frequencies and slower for lower ones. This is shown explicitly for the model reduction approach in Figure 2 (right), where the convergence of Krylov approximations of a fixed order is shown for different frequencies within a frequency range running from 2.5Hz to 7.5Hz. At this point, we stress that the principal difference between the PKS method and BiCGStab is that the former approach obtains solutions for the whole given frequency range at the convergence cost of BiCGStab for the lowest frequency (from that range).
Subsequently, we consider what adding negative powers gives us in terms of performance. In Figure 3 (left), we have plotted convergence curves (with respect to increasing k) for approximants obtained using EKS K k,ki+1 for different values of fixed i (i = ∞ corresponds to PKS) and for a frequency range that runs from from 2.5Hz to 7.5Hz. From this figure it is clear that adding negative powers (i < ∞) visibly improves EKS convergence compared with the PKS method. In terms of computation time, it took PKS 157500 seconds (4500 iterations) to reach 10 −4 accuracy, while EKS with i = 3, i = 5, and i = 7 required 154000 seconds (2300 iterations), 189000 seconds (3000 iterations), and 156000 seconds (300 iterations), respectively, to reach the same level of accuracy. From these experiments, it is clear that PKS converges faster for higher frequencies, while EKS improves convergence for smaller frequencies and, consequently, makes the convergence more uniform for the whole frequency range. This statement is confirmed in Figure 3 (right), where we have plotted the convergence curves for the EKS method with the same subspace parameters as before, but this time for a frequency range that runs from 0.1Hz to 7.5Hz. We observe that PKS significantly slows down in this case, while EKS performs almost as efficient as for the narrower frequency range. For example, to reach a 10 −4 accuracy, it takes about the same time for EKS to converge, while PKS hardly reached 10 −2 accuracy in 525000 seconds of computational time (15000 iterations).    shown in Figure 4. This configuration was also considered in [2], where the PKS method was used to compute the field distributions. We revisit this problem here so that we can compare the performances of EKS and PKS. As an independent check, we also compare our results with FDTD computations. Since FDTD is based on the first-order Maxwell system, we use a first-order formulation of the EKS method in our time-domain experiments as well.
The configuration of Figure 4 consists of a set of dielectric rods with vacuum as a background medium. The rods have a relative permittivity of 11.56 and the distance a between the rods is approximately 0.58 µm. The radius of each rod is 0.18a. By removing rods inside this crystal we create a bend as is illustrated in Figure 4. A source is placed at the corner of the bend and we record the electric field strength halfway one of the bends. The source pulse is a modulated Gaussian with its spectrum essentially contained in the frequency interval 9.8 · 10 14 ≤ ω ≤ 1.44 · 10 15 rad/s, which coincides with the ideal frequency bandgap of this structure.
In Figure 5, we present the electric field response at the receiver location for a diagonal extended Krylov subspace (m 1 = m 2 ) after 100 iterations with the system matrix and 100 iterations with its inverse (black line). The field response obtained after 8197 iterations of the FDTD method run at the Courant limit is also shown in Figure 5. Clearly, the EKS reduced-order model has not converged yet on the time interval of interest. We therefore double the number of iterations and carry out 200 iterations with the system matrix and 200 with its inverse. The resulting EKS reduced-order model for the electric field strength is plotted in Figure 6 along with the corresponding FDTD result. Even though there is some improvement, the reduced-order model still has not converged on our time interval of observation. Finally, Figure 4: Crystal consisting of dielectric rods. The cross and the triangle indicate the location of the source and receiver, respectively. The rods have a relative permittivity of 11.56. The distance a between the rods is 0.58 µm and the radius of each rod is 0.18a.
carrying out 400 iterations with the system matrix and 400 with its inverse we obtain the reduced-order model as shown in Figure 7 (black line). The model essentially overlaps with the FDTD result in this case and the number of iterations required by the EKS method is clearly much smaller than the number of FDTD iterations (800 for EKS, 8197 for FDTD). However, the amount of work that is involved at each iteration is not the same of course, since EKS requires matrix inversion, while FDTD does not. We therefore compare the computation times of both methods as well and find that FDTD requires 1764 seconds to solve this problem, while the EKS method requires only 888 seconds. In terms of offline costs, we have achieved a speed up of approximately a factor of two compared with FDTD. If we compare this data with the results presented in [2], we observe that EKS is actually slower than PKS in this case, since according to [2], PKS requires only 606 seconds to solve this problem. However, the order of the PKS reduced-order model is 3000 (see [2]), which is much larger than the order of the EKS reduced-order model (which is 800). We conclude that for this problem both EKS and PKS are faster than FDTD. When reduction of the offline costs is important then PKS should be used, since it has the shortest computation time, but if the online costs are dominant then the EKS model should be used, since it provides us with the smallest reduced-order model.

Conclusions
Based on model reduction techniques, we have developed a powerful tool for solving large-scale multi-frequency wave problems. The PKS method allows us to obtain wave field solutions on a whole frequency range at the cost of an unpreconditioned linear solver applied to the problem with the smallest frequency. Moreover, we have shown that EKSs significantly outperform the polynomial approach when solutions for small frequencies are required and we have demonstrated that the EKS approach allows for substantial dimensionality reduction in the time-domain. The time-domain models produced by EKS have a much smaller order than the models produced by PKS. On the other hand, subspace creation via EKS takes more time than subspace creation via PKS. Based on these findings, we conclude that PKS should be used if reduction of the offline computational costs is desired, while EKS is the method of choice when the online costs need to be minimized.
Finally we note that our approach is not limited to second-order schemes in the interior. In fact, since the optimal discrete PML has spectral accuracy (see [2,3]), it would be preferable to use a discretization scheme that has spectral accuracy in the interior as well. High-order spectral methods can be applied, for example, (see [15]) or optimal grids for interior domains as described in [16].

Acknowledgements
We'd like to thank Dr. Guangdong Pan for providing us with results for his BiCGStab finite-difference solver.