Least-Squares Pad\'e approximation of parametric and stochastic Helmholtz maps

The present work deals with the rational model order reduction method based on the single-point Least-Square (LS) Pad\'e approximation technique introduced in [3]. Algorithmical aspects concerning the construction of the rational LS-Pad\'e approximant are described. In particular, the computation of the Pad\'e denominator is reduced to the calculation of the eigenvector corresponding to the minimal eigenvalue of a Gramian matrix. The LS-Pad\'e technique is employed to approximate the frequency response map associated to various parametric time-harmonic wave problems, namely, a transmission/reflection problem, a scattering problem and a problem in high-frequency regime. In all cases we establish the meromorphy of the frequency response map. The Helmholtz equation with stochastic wavenumber is also considered. In particular, for Lipschitz functionals of the solution, and their corresponding probability measures, we establish weak convergence of the measure derived from the LS-Pad\'e approximant to the true one. 2D numerical tests are performed, which confirm the effectiveness of the approximation method.


Introduction
Many applications require the fast and accurate numerical evaluation of Helmholtz frequency response functions, i.e., functions that map the wavenumber to the solution (or some quantity of interest related to the solution) of the corresponding time-harmonic wave-problem, for a large number of frequencies. In mid-and high-frequency regimes, very fine meshes or high polynomial degrees should be considered, in order to obtain accurate Finite Element (FE) solutions of the time-harmonic wave-problem. Moreover, low order FE schemes are affected by the pollution effect [2], namely, an increasing discrepancy between the best approximation error and the FE error, as the wave number increases. In the "many-queries" context, i.e., when many solutions of the underlying Partial Differential Equation (PDE) are needed, the "brute force" approach entails the solution of a large number of high-dimensional linear systems, and it is then out of reach.
Model order reduction methods aim at significantly reducing the computational cost by approximating the quantity of interest starting from evaluations at only few wavenumbers. They rely on a two-step strategy: the OFFLINE stage consists in the computation of a finite dimensional basis -e.g., the basis of snapshots (see, e.g., [5,14,19,25,26,28,29,15,8,22]), or evaluations of the frequency response map and its derivatives at fixed centers (Padé method, see, e.g., [7,13,10,9,3]); the output of this phase, whose computational cost may be very high, is stored, to be used during the ONLINE phase, in which the approximation of the frequency response map corresponding to a given new value of the parameter is constructed. This stage does not involve the numerical solution of any PDE, and is expected to provide the output in real time.
In this work, we focus on the Padé-based model order reduction technique introduced in In this paper, we describe in detail the algorithmical aspects of the construction of the single-point LS-Padé approximant. In particular, the identification of the LS-Padé denominator is proved to be equivalent to the identification of the normalized eigenvector corresponding to the smallest non-negative eigenvalue of the Gramian matrix of the set T (z 0 ), T 1,z0 , . . . , T N,z0 , where T α,z0 denotes the Taylor coefficient of T of order α at z 0 .
Moreover, we explore the effectiveness of the single-point LS-Padé technique when applied to parametric frequency response problems which go beyond the setting considered in [3], namely, a transmission/reflection problem, and a scattering problem. In both cases, we first prove that the frequency response map associated with the considered problem is meromorphic. 2D numerical results are provided, which demonstrate the convergence of the LS-Padé approximation. Moreover, 2D numerical tests in highfrequency regime are performed for the parametric problem presented in [3].
The stochastic Helmholtz boundary value problem is also considered. We refer to [23] for Uncertainty Quantification for frequency responses in vibroacoustics, to [6,17,16] for model order reduction for random frequency responses in structural dynamics, and to [24,12] for the stochastic Helmholtz equation with uncertainty arising either in the forcing term or in the boundary data or in the shape of the scatterer.
Within the present framework, we propose a novel approach to the stochastic Helm-holtz boundary value problem based on the LS-Padé technique, where the wavenumber k 2 is modeled as a random variable taking values into K = [k 2 min , k 2 max ]. We approximate the random variable X := L(S(k 2 )) with X P := L(S [M/N ] (k 2 )). Here, L : V → R is a Lipschitz functional representing a quantity of interest, S is the meromorphic frequency response map associated with the (stochastic) Helmholtz equation endowed with either homogeneous Dirichlet or homogeneous Neumann boundary conditions, and S [M/N ] is the LS-Padé approximation of S. An upper bound on the approximation error for the characteristic function is derived.
All the considered boundary value problems fall into the following general setting. Let D be an open connected bounded Lipschitz domain in R d (d = 1, 2, 3), and consider the following Helmholtz boundary value problem where the wavenumber k 2 is either a parameter or a random variable, which takes values into an interval of interest K : . Moreover we assume the functions in V to be complex-valued.
The outline of the paper is the following. In Section 2, we recall the definition of the single-point LS-Padé approximant and the main convergence result of [3]. In Section 3, we describe the algorithm to compute the LS-Padé approximant. Section 4 deals with a parametric transmission/reflection problem, whereas Section 5 deals with a parametric scattering problem. In Section 6, the LS-Padé approximation is tested in high-frequency regime, and in Section 7 the LS-Padé methodology is applied to the stochastic setting. Finally, conclusions are drawn in Section 8.

Least-Squares Padé approximant of the parametric model problem
This section deals with the Least-Squares (LS) Padé approximation of the following parametric Helmholtz problem: Problem 1 (Parametric Model Problem) The Helmholtz equation (1) has parametric wavenumber k 2 ∈ K := [k 2 min , k 2 max ] ⊂ R + , ε r = 1, and is endowed with either Dirichlet or Neumann homogeneous boundary conditions on ∂D, i.e., Γ R = ∅ and either Γ D = ∂D and g D = 0, or Γ N = ∂D and g N = 0.
The following result was proved in [3].
Theorem 2.1 Let S be the frequency response map which associates to each z ∈ C, the solution u z ∈ V of the weak formulation of Problem 1: Then, S is well-defined, i.e., problem (2) admits a unique solution for any z ∈ C \ Λ, Λ being the set of (real, non negative) eigenvalues of the Laplace operator with the considered boundary conditions. Moreover, S is meromorphic in C, with a pole of order one in each λ ∈ Λ. We recall now the definition and the convergence theorem of the LS-Padé approximant of the frequency response map S.
Let K = [k 2 min , k 2 max ] ⊂ R + be the interval of interest, and z 0 ∈ C \ Λ with Re (z 0 ) > 0. To fix the ideas we take z 0 = k 2 min +k 2 max 2 +δi, with δ ∈ R\{0} arbitrary. The LS-Padé approximant of S, centered in z 0 , is given by the ratio of two polynomials of degree M and N respectively: The denominator Q [M/N ] (z) is a function of z only, and belongs to the space P N (C) of all polynomials of degree at most In the following, we denote with P M (C; V ) the space of polynomials of degree at most M in z ∈ C with coefficients in V .
The construction of the LS-Padé approximant proposed in [3] relies on the minimization of the functional j E,ρ : P M (C; V ) × P N (C) → R, parametric in E ∈ N and ρ ∈ R + , defined as where the brackets · α,z0 denote the α-th Taylor coefficient of the Taylor series centered in z 0 (i.e., for a map T : dz α (z 0 )), and · V, √ Re(z0) denotes the weighted H 1 (D)-norm (equivalent to the standard one) defined as We recall the formal definition of the LS-Padé approximant of the solution map S, and we refer to Section 3 for the proof of the existence of a (not in general unique) LS-Padé approximant.
The following convergence result has been proved in [3].
Theorem 2.4 Let N ∈ N be fixed, and let R ∈ R + be such that the disk B(z 0 , R) contains exactly N poles of S. Then, for any z ∈ B(z 0 , R) \ Λ and for any |z| < ρ < R, it holds lim Then for any 0 < ρ < R such that B(z 0 , ρ) ⊃ K, there exists M ∈ N such that, for any M ≥ M and for any z ∈ K \ K α , it holds where the constant C > 0 depends on ρ, R, N , z 0 , λ min = min{λ ∈ Λ}, f L 2 (D) , and g(z) = λ∈Λ∩B(z0,R) (z − λ).

Remark 2.5 In [3] the bound
was proved, with a constant C that depends on and g(z) = λ∈Λ∩B(z0,R) (z − λ). Since the frequency response map S presents only simple poles (given by the Dirichlet/Neumann Laplace eigenvalues), and the interval of interest K contains a finite number of poles of S, it follows that there exists C g > 0 such that |g(z)| ≥ C g min λ∈Λ∩B(z0,R) |z − λ| ∀z ∈ K, hence g K\Kα ≥ C g α and the bound (7) follows.
We can draw the following consequences:

Algorithmical aspects
In this section, we describe an algorithm for the computation of a LS-Padé approximant (defined according to Definition 2.3) of the Helmholtz frequency response map S introduced in the previous section. We underline that the presented algorithm can be likewise applied to any V -valued meromorphic map T : C → V . As a first instructive step, we recall the proof of the existence of such an approximant, which was developed in [3, Proposition 4.1]. Proof. We want to show that the minimization problem (6) admits at least one solution.
Since P has degree M , then P (z) α,z 0 = 0 for all α > M . Hence, we can rewrite jE,ρ as problem (6) can be formulated as a minimization problem in Q only: find Q ∈ P N (C) such thatj Since the functionaljE,ρ is continuous and the set P N (C) is compact (being homeomorphic to the unit sphere in C N +1 ),jE,ρ has a global minimum on P N (C), and the minimization problem (8) admits at least one solution.
In the following proposition we express an equivalent formulation of the constrained minimization problem (8). (8) is equivalent to the identification of the (normalized) eigenvector corresponding to the smallest non-negative eigenvalue of the Hermitian positive-semidefinite matrix G E,ρ ∈ C (N +1)×(N +1) with entries
The Hermitian matrix G E,ρ defined in (10) is obtained as weighted sum of submatrices of the Gram matrix G ∈ C (N +1)×(N +1) associated with the solution map S, namely, the matrix with entries , for i, j = 0, . . . , N .
See Figure 1 for a graphical representation. By following the steps performed in the proof of Proposition 3.1, and applying Proposition 3.2, we devise Algorithm 1 for the computation of the LS-Padé approximant.

Remark 3.3
The choice of ρ impacts the algorithm only by determining the weights in the computation of G E,ρ . Specifically, small (respectively large) values of ρ emphasize the contributions from the sub-matrices located in the top-left (respectively bottom-right) portion of G. A fast version of the algorithm, where G E,ρ reduces just to the leading term (i.e., for ρ → +∞), is currently under investigation (see [4]). : Gram matrix (top) associated with the frequency response map S. To lighten the notation, we omit both the argument (z 0 ) of the Taylor coefficients S α , and the weight Re (z 0 ) of the scalar product ·, · V . In blue the sub-matrix corresponding to N = 2 and α = 3, which provides a contribution to G E,ρ (bottom) with weight ρ 6 .
Observe that a transposition with respect to the secondary diagonal is carried out before computing the sum.

Application to a transmission/reflection problem
We consider the transmission/reflection problem treated in [18], i.e., the transmission/reflection of a plane wave e iκx·d with wavenumber κ and direction d = (cos(θ), sin(θ)), across a fluid-fluid interface. In particular, the considered domain D = (−1, 1) 2 is divided into two regions with different refractive indices n 1 , n 2 ; we assume n 1 < n 2 . The Helmholtz problem is the following For any angle 0 ≤ θ < π/2, the following function is a solution of equation (12): K2+κn1d2 and T = 1 + R. We couple the Helmholtz equation (12) with Dirichlet boundary conditions derived from the exact solution (13), i.e., u| ∂D = u ex | ∂D .
Depending on the value of θ (angle of the incident wave), the solution may exhibit two types of behavior: • if θ < θ crit := arccos n2 n1 , then Im (K 2 ) = 0, and u ex decays exponentially for x 2 > 0. Physically, this phenomenon is called total internal reflection; • if θ > θ crit , then d is close to the normal incidence, and the wave is refracted at the interface.
The two behaviors are depicted in Figure 2.

Frequency response map
We are interested in the following boundary value problem: Problem 2 (Transmission/Reflection Problem) The wavenumber κ 2 ranges in the interval of interest K = [κ 2 min , κ 2 max ], and the Helmholtz equation is endowed with Dirichlet boundary conditions on Γ D = ∂D: where g D := u ex | ∂D , and u ex is given by formula (13) with κ = 11 and either θ = 29 • or θ = 69 • .

LS-Padé approximant of the frequency response map
Since the frequency response map is meromorphic, it is appropriate to use the LS-Padé technology to catch the singularities of S, and provide sharp approximations of S(z), when z is close to the center z 0 . We apply Algorithm 1, and compute the coefficients of the denominator as the entries of the eigenvector corresponding to the minimal eigenvalue of the Gram matrix (10). The Taylor coefficient of order β ≥ 1, Problem (20) admits a unique solution for all z ∈ C \ Λ, since the PDE operator is the same as in (15) and the right-hand side is a bounded linear form. Let K = [3,12] be the interval of interest and θ = 29 • . In Figure 3, the H 1 (D)weighted norm of the P 2 finite element approximation of S, S h , is compared with the norm of its LS-Padé approximant S h,P centered in z 0 = 7.5 + 0.5i, for various degrees (M, N ). We have empirically observed (see Figure 4) that the LS-Padé approximation delivers a better accuracy than that predicted in (7): where {λ j } j are the elements of Λ ordered according to: |λ 1 − z 0 | < |λ 2 − z 0 | < . . .. We refer to [4] for a formal derivation of (21), where S [M/N ] is computed by a fast version of Algorithm 1.

Application to a scattering problem
In this section, we consider the scattering of an acoustic wave at a scatter occupying the domain B((0, 0), 0.5) ⊂ R 2 . The incident wave u i is the time-harmonic plane wave traveling along the direction d = (cos(θ), sin(θ)) with wavenumber k, i.e., u i = e ikd·x . The total field u, given by the sum of the incident wave u i with the scattered wave u s , satisfies the following boundary value problem in the infinite domain R 2 \B ((0, 0) The finite element approximation of problem (22)  whose outer boundary will be denoted as Γ R . Approximating the Sommerfeld radiation condition at infinity in problem (22) by a first order absorbing boundary condition, we write the following parametric problem: Problem 3 (Scattering Problem) The wavenumber k 2 ranges in the interval of interest K := [k 2 min , k 2 max ] ⊂ R + , n is the outgoing normal vector field to Γ R , and g R := ∂u i ∂n − iku i is the impedance trace of the incoming wave u i . We consider the Helmholtz boundary value problem

Regularity of the frequency response map
We extend problem (23) to complex wavenumbers. Given a complex wavenumber z ∈ C, we introduce the incident plane wave u i = e izd·x and its impedance trace g z := ∂u i ∂n − izu i , and we define the frequency response map S : z → S(z) If z ∈ R, problem (24) admits a unique solution (see, e.g., [11]), which implies that the frequency response map is well-defined on R. The following Theorem extends this result to the complex half plane {z ∈ C : Im (z) ≥ 0}. Since the wavenumber in (24) is square of the parameter z, we will endow the Hilbert space V with the weighted H 1 (D)-norm, with weight w = Re (z 0 ) (and not w = Re (z 0 ), as was done before).
We recall here the following theorem, see [27,Theorem 1], which will be used in the proof of Proposition 5.3. Proof. We proceed as in [20,Proposition 2]. We add and subtract the term D uzvdx to the left-hand side of (24), and we get which can be written equivalently as where T (z), Gz : V → V are defined, respectively, as where Ctr is the continuity constant of the trace operator γ : H 1/2+ε (D) → L 2 (∂D) (see, e.g., [1, Theorem 5.36]). Applying Theorem 5.2, we conclude that (I − T (z)) −1 is meromorphic in all open bounded and connected subsets of C and, since Gz is linear in z (hence holomorphic in C), the same conclusion applies to the frequency response function S(z) = (I − T (z)) −1 Gz. Moreover, since Theorem 5.1 states that S is well defined in C + , we deduce that all poles of S must have negative imaginary part.

LS-Padé approximant of the frequency response map
We construct the LS-Padé approximant of the frequency response map S following Algorithm 1. Having fixed z 0 ∈ C + , N, M , and E ≥ M + N , the coefficients of the denominator are computed by identifying the eigenvector corresponding to the smallest eigenvalue of the matrix (10), where the β-th Taylor coefficient of S, S β,z0 , solves the following recursive problem: Since the PDE operator in (29) is the same as in (24), and the linear form at the right-hand side is bounded, by applying Theorem 5.1, we conclude that problem (29) is well-posed for any z ∈ C + . Let z = 3, and z 0 = 3 + 0.5i. In Figure 6 (left) we plot the relative LS-Padé approximation error versus the degree of the LS-Padé numerator, for different values of denominator degree. In Figure 6 (right), the relative error obtained by approximating the frequency response map with the Taylor polynomial (black dashed line), and with the LS-Padé approximant are compared. Also the diagonal LS-Padé approximant is considered (dashed purple line), where the LS-Padé numerator and denominator have the same degree. In Figure 6 (right), the errors are plotted versus the number of derivatives S β,z0 , β = 0, . . . , E computed (i.e., the number of PDEs solved offline). Since B(z 0 , |z − z 0 |), the disk with center z 0 = 3 + 0.5i and radius r 1 = |z − z 0 | = 0.5, is contained in the half plane where the frequency response map is holomorphic (see Proposition 5.3), the Taylor series centered in z 0 converges, and the Taylor approximation error is comparable to the LS-Padé approximation error. Figure 7 presents analogous plots as in Figure 6, for the point z = 2. In this case, B(z 0 , |z − z 0 |) ∩ {Im (z) < 0} = ∅, and the Taylor series diverges.
In Figure 8, we plot the numerical solution S h ∈ V of problem (2) with z = 51, computed via P 3 continuous finite elements. Observe that the relative finite element error is of the order of 10 −5 . In Figure 9, the LS-Padé approximant S h,P centered in z 0 = 47 + 0.5i evaluated in z = 51 is represented for two different values of the denominator degree. Due to the fact that more derivatives are employed in the right plot, more accurate results are obtained with higher denominator polynomial degrees.
In Figure 10, we plot the LS-Padé approximation error w.r.t. the exact solution, in z = 51, for different values of the degree of the denominator, and we compare it with the numerical rate (21). When N = 2, 4, the LS-Padé technique works as expected (or even better), whereas for N = 6 the error is no longer decreasing. We believe that this behavior is caused by the ill-conditioning of Step 7 in Algorithm 1), i.e., the computation of the (normalized) eigenvector of the Gramian matrix G E,ρ defined in (10).
We partition uniformly the interval of interest K in 100 subintervals. At each point z of the grid we have computed the numerical solution S h (z) of the Helmholtz problem (2), and the LS-Padé approximant S h,P (see Figure 11), as well as the relative error Figure 12). In Figure 13,

LS-Padé approximant of the stochastic model problem
This section deals with the stochastic counterpart of Problem 1: Problem 4 (Stochastic Model Problem) The wavenumber k 2 of the Helmholtz equation is modeled as a random variable with bounded density function F k 2 . In this section, either Dirichlet or Neumann or mixed Dirichlet/Neumann homogeneous boundary conditions on ∂D are considered.
We introduce a Lipschitz functional L : V → R representing a quantity of interest of the frequency response map S, and we define the following two random variables: and Im (z 0 ) = 0; this guarantees that z 0 / ∈ Λ, Λ being the set of eigenvalues of the Laplacian, with the considered boundary conditions. Let φ X , φ X P : R → C denote the characteristic functions of X and X P , respectively, i.e., φ X (t) := E e itX , φ X P (t) := E e itX P . We are interested in studying the LS-Padé approximation error on the characteristic function, i.e., we aim at proving an a priori bound for (32) Theorem 7.1 Let L : V → R be a Lipschitz functional with Lipschitz constant L, and let X, X P be the random variables defined in (30) and (31). Given α > 0, then it holds with the same definitions of R, ρ, and K α , and the same characterization of C > 0 as in Theorem 2.4, and |·| denoting the Lebesgue measure.
Proof. Using the definition of the characteristic function and the linearity of the expected value we find We bound the two integrals separately. For the integral over Kα, we have The conclusion follows from inequalities (34) and (35). In particular, there exists C > 0 such that for any t ∈ R err t ≤ C |t| This corollary establishes, in particular, uniform exponential convergence of φ X P to φ X on any compact subset of R. Remark 7.3 Theorem 7.1 and Corollary 7.2 can be generalized to derive an a priori upper bound on err ξ = |E [ξ(X)] − E [ξ(X P )]|, for any continuous and bounded functional ξ : R → R. Thus, the weak convergence of X P to X follows, as M → +∞.
Let us consider the case of ∂D = Γ D . Let K = [7,14] be the interval of interest (which contains three eigenvalues of the Dirichlet-Laplace operator: 8, 10, 13), and let the wavenumber be modeled as a random variable uniformly distributed on K, i.e., k 2 ∼ U(K). Given the functional L = · V, √ Re(z0) , where z 0 = 10+0.5i, we consider the random variables X = S h (k 2 ) V, √ Re(z0) and X P = S h,P (k 2 ) V, √ Re(z0) . We define as S h the P 3 finite element approximation of S; then S h,P is the LS-Padé approximant of S h , centered in z 0 and with polynomial degrees (M, N ). In Figure 14, we display the random variables X and X P evaluated at 100 sample points uniformly distributed in K. When the degree of the LS-Padé denominator is N = 3, all the poles are correctly identified by the LS-Padé approximant, provided that M is larger than 4. In Figure 15 we plot the characteristic function of the random variable X P , φ X P (t), where the degrees of the Padé denominator and denominator are N = 3, and M = 2, 4, 6, respectively. The expected value has been computed by the Monte Carlo method, using 10 5 samples.

Conclusions
The present paper concerns a model order reduction method based on the single-point LS-Padé approximation technique introduced in [3]. We have described an algorithm to compute the LS-Padé approximant of the Helmholtz frequency response map, and we have explored the applicability and potentiality of the method via 2D numerical experiments in various contexts. Moreover, the time-harmonic wave equation with random wavenumber has been analyzed. We are currently investigating the extension of the proposed methodology and of its convergence analysis to the case of multi-point LS-Padé expansions, where evaluations of the frequency response map S and of its derivatives at multiple frequencies are used. We believe that this technique will outperform the single-point one, when a large number of singularities of S need to be identified.