Solving inverse scattering problems via reduced-order model embedding procedures

We present a reduced-order model (ROM) methodology for inverse scattering problems in which the reduced-order models are data-driven, i.e. they are constructed directly from data gathered by sensors. Moreover, the entries of the ROM contain localised information about the coefficients of the wave equation. We solve the inverse problem by embedding the ROM in physical space. Such an approach is also followed in the theory of ``optimal grids,'' where the ROMs are interpreted as two-point finite-difference discretisations of an underlying set of equations of a first-order continuous system on this special grid. Here, we extend this line of work to wave equations and introduce a new embedding technique, which we call Krein embedding, since it is inspired by Krein's seminal work on vibrations of a string. In this embedding approach, an adaptive grid and a set of medium parameters can be directly extracted from a ROM and we show that several limitations of optimal grid embeddings can be avoided. Furthermore, we show how Krein embedding is connected to classical optimal grid embedding and that convergence results for optimal grids can be extended to this novel embedding approach. Finally, we also briefly discuss Krein embedding for open domains, that is, semi-infinite domains that extend to infinity in one direction.


Introduction
In this paper we discuss so-called reduced-order model (ROM) embedding procedures to solve inverse scattering problems.In such a procedure, data-driven ROMs are constructed from spectral impedance data collected at one end of a bounded interval of interest.Subsequently, a ROM is interpreted as a two-point finite-difference discretisation of an underlying set of first-order continuous wave equations.We refer to this step as embedding of the ROM in physical space and it is this embedding procedure that allows us to determine the medium parameters on the spatial interval of interest.
Two embedding procedures are discussed in this paper, the first, called optimal grid embedding, is based on the theory of (optimal) truncated spectral measure grids [1].We consider the recovery of a velocity profile from the first poles and residues of the impedance function, the so-called truncated spectral measure.An implementation of this procedure in terms of travel time coordinates is presented in [2].Here, we take this optimal grid procedure as a starting point and present optimal grid embedding in terms of spatial coordinates instead of travel time coordinates.
However, a drawback of optimal grid embedding is that it requires training for a known medium, data, and boundary conditions.More specifically, to retrieve the spatially varying medium parameters on a certain bounded interval of interest, we first have to determine an optimal grid for a homogeneous reference medium.Having this trained reference grid available, the position dependent medium parameters can be determined on the interval of interest.Moreover, the length of this interval must be included in training as well, but this information is not always available.In radar imaging, for example, the distance to the surface may not be known exactly and for problems on semi-infinite domains, the step sizes quickly become very large [4].
The use of a reference grid is avoided in the second embedding procedure presented in this paper, which is a new procedure that we call Krein embedding, since it is inspired by Krein's seminal work on vibrations of a string [5,6].Furthermore, the length of the interval of interest is not required either.We also show how Krein embedding is related to optimal grid embedding and that convergence results presented in [1] for optimal truncated spectral measure grids essentially carry over to the Krein embedding approach.
Since a trained reference grid and the length of a reconstruction interval are not required for Krein embedding, such an embedding approach may be applied to scattering problems on semi-infinite domains as well, that is, domains that extend to infinity in one direction.Furthermore, scattering problems on semi-infinite domains characterised by continuous spectral measures can equivalently be described by scattering poles that correspond to passive dissipative systems [9].Our approach is then to construct ROMs that can be interpreted as a finite-difference discretisation of a dissipative first-order system.We call this Krein-Nudelman embedding, since the case of a dissipative boundary condition was discussed by Krein and Nudelman in [6].
However, for problems on semi-infinite domains a uniqueness problem arises, since the method embeds an impedance function provided in pole-residue form, which cannot distinguish between a lossy bounded domain and a lossless open domain.In the latter case, the spectrum is not a point spectrum but is represented as such.Nevertheless, if we apply a Krein-Nudelman embedding approach to such a problem, we find that the medium profile is actually recovered up until the last reflector, where Krein-Nudelman embedding places an absorbing effective medium to match the (complex) point spectrum of the impedance function.A numerical example will be presented that illustrates this phenomenon and the uniqueness problem is discussed further in the Appendix.The remainder of this paper is organised as follows.In Section 2 we discuss the construction of the ROMs from spectral data and the optimal grid and Krein embedding procedures that may be used to retrieve the medium parameters on a bounded interval.Subsequently, these two procedures are discussed in detail in Sections 3 and 4, while in Section 5 a number of numerical examples are presented that illustrate the performance of both embedding procedures.Finally, Krein-Nudelman embedding on a semi-infinite domain is briefly discussed in Section 6 and the conclusions can be found in Section 7.

Embedding of reduced-order models
We are interested in reconstructing the wave speed c in the wave equation from boundary measurements.We formulate the problem in the temporal Laplace domain with complex Laplace frequency s ∈ C and reflecting boundary conditions.Specifically, on the bounded interval [0, L] of interest, the governing Laplace-domain equation is given by This problem is equivalent to the equation for a vibrating string studied by Kac and Krein [5] if c −2 (x) is replaced by the mass density of a string.The above equation also follows from the telegrapher equations for a transmission line that is short-circuited at the far-end of the line and a unit current is fed into the near-end of the line.In this case, c(x) represents the wave speed along the transmission line.For a given wave speed profile, equation (1) essentially represents a regular Sturm-Liouville problem.We are also interested in the corresponding singular case in which equation ( 1) is considered on the semi-infinite interval [0, ∞) with the boundary condition for u at x = L replaced by the condition that |u(x)| → 0 as x → ∞ for Re(s) > 0.
In both the regular and singular case, the associated inverse problem is to recover c(x) from measurements of the impedance function f (s) = u(0, s) for s on some curve in the complex plane.The spectral inverse problem that we consider in this paper consists of the reconstruction of c(x) from poles and residues of the impedance function f (s).In particular, we assume to have access to the first n complex-conjugate pairs of poles λ j and residues y j of the function f (s).In an application, this spectral information is not readily available from measurements of f (s).However, it can be retrieved from the measured impedance function f (s) using the vectorfit algorithm [3].In the numerical experiment section, we illustrate this procedure.
Assuming that we have n complex-conjugate pole-residue pairs available, we can construct the reduced-order model where the overbar denotes complex conjugation.For a regular Sturm-Liouville problem, the residues y j are real and positive and the λ j are purely imaginary.However, for the singular Sturm-Liouville case discussed later this may not be the case and we therefore include conjugation in the above spectral expansion of f ROM (s).Finally, we mention that if we introduce the diagonal matrix and the residue vector the reduced-order model can also be written as where I is the 2n × 2n identity matrix.

Building a reduced-order model from spectral data
The key idea behind optimal grid and Krein embedding is to interpret the reduced-order model f ROM (s) as the impedance function of a two-point finite-difference discretisation of an underlying set of first-order differential equations.Specifically, introducing a dual variable û and the staggered grid shown in Figure 1, the first-order finite-difference system that corresponds to equation ( 1) is given by ûj+1 − ûj γj+1 + su j = 0, ∀j = 0, . . ., n − 1, with u j = u(x j ) and ûj = û(x j ) and where γ j and γj are the edge weights, i.e. products of the step sizes x j − x j−1 and the medium parameters on the grid.The boundary conditions of the system are given by û(0) = û0 = 1 and u(L) = u n = 0. Introducing the 2n-by-1 vector of unknowns the finite-difference system of equation ( 6) can be written compactly as where e 1 is the first canonical basis vector of length 2n and matrix T is a tridiagonal matrix of order 2n given by To accommodate a standard complex-symmetric Lanczos implementation further on, we prefer to work with the transpose-symmetric matrix instead of the tridiagonal matrix T. Matrix T TS is related to matrix T by the similarity transform T TS = STS −1 , where the similarity matrix S is given by With f ROM (s) = u 0 (s) = e T 1 u we now find that by solving the system of equations ( 8) for u that To extract the finite difference parameters γ j , γj from the spectral parameters (λ j , y j ) we need to transform the diagonal transfer function representation of equation ( 5) to the tridiagonal one of equation (12).To that end, let Y denote the eigenvector matrix of T TS such that it satisfies T TS Y = YΣ with Y T Y = I and Σ is a diagonal matrix of order 2n with the eigenvalues of T TS on the diagonal.The impedance function can now be written as where y 1 = Y T e 1 is the vector containing the first components of all eigenvectors.
Setting y 1 = √ γ1 y and with Σ = Λ, we observe that the above reduced-order model coincides with the reduced-order model of equation (5).Since the residue vector y and the poles Λ are known, the problem turns into an inverse eigenvalue problem in which we attempt to reconstruct the tridiagonal matrix T TS from the first components of its eigenvectors and its known eigenvalues.As is well known, this problem can be solved via the Lanczos algorithm.Specifically, with orthogonalization in the transpose bilinear form [7,8], vector y/ y T y as a starting vector, and the matrix of eigenvalues Λ as iteration matrix, the Lanczos algorithm produces the desired tridiagonal matrix.The algorithm is given in Algorithm 1 and as soon as the tridiagonal matrix is obtained, the coefficients γ j and γj can be extracted using the recursive scheme shown in Algorithm 2. We note that the weight γ1 can be determined directly from the residues, since Finally, we also note that the algorithm explicitly computes the diagonal elements α j of matrix T TS , while we know that these elements vanish for lossless media.However, this very same algorithm can be used for problems on open domains or for problems involving lossy media for which the diagonal elements of matrix T TS may no longer vanish.Therefore, we prefer to work with this Lanczos algorithm since it accommodates all cases of interest.
Algorithm 1 The Lanczos algorithm for complex symmetric matrices to obtain the tridiagonal matrix T TS from the poles Λ and residues y.In this algorithm, Y j indicates the jth column of the eigenvector matrix Y, α j is the jth element on the diagonal (T TS ) j,j and β j the jth element on the super diagonal (T TS ) j−1,j of matrix T TS .

Krein and optimal grid embedding
Up to this point we constructed a ROM from the poles and residues of a boundary impedance function.The ROM can be interpreted as the impedance function of a two-point finite-difference discretization of a first-order system with primary and dual grid coefficients γ j and γj , respectively.How to interpret these coefficients depends on the dual variable û that is introduced to obtain a first-order system of ODEs from the original second-order ODE.In particular, introducing the dual variable û = −s −1 ∂ x u, the second-order system of equation ( 1) can be written as and this form will lead to what we call the Krein embedding interpretation.On the other hand, introducing the quantities w = c − 1 2 u and ŵ = c 1 2 û as well as the slowness coordinates (sometimes called travel time coordinates) we obtain the first-order system with w(T (L), s) = 0 and ŵ(0, s) = c 1 2 (0), and this system will leads to a standard optimal grid embedding interpretation.Note that in this latter case the impedance function is f (s) = u(0, s) = c 1 2 (0)w(0, s) and that standard optimal grid embedding is formulated in terms of travel time coordinates.
In the next two sections we will give detailed descriptions of optimal grid and Krein embedding procedures for wave propagation problems.Convergence of these procedures is also briefly discussed.

Optimal grid embedding
Optimal grid embedding of ROMs was developed for the diffusion equation in [1].In this section, we extend this embedding approach to wave propagation problems and rely on results obtained for the diffusion equation.
The main difficulty in ROM embedding is that each of the finite difference weights γ j and γj consist of a product of the unknown local medium parameter and an unknown grid step.Fortunately, for diffusion problems it has been shown that there exists a computable grid {γ 0 j , γ0 j } that is independent of the medium parameters in the limit n → ∞ (Lemma 3.2 of [1]).The existence of this grid allows us to obtain the local material parameters by taking ratios of the grid steps γ 0 j , γ0 j and ROM parameters γ j , γj thereby reconstructing the medium.Furthermore, it can be shown the medium estimates converge pointwise in L 1 (0, L) (Theorem 6.1 of [1]).
As a first step towards optimal grid embedding for wave propagation, we extend the results from [1] by considering equation ( 17) in spatial coordinates instead of travel time coordinates.Subsequently, we take the kinematic effects into account and discuss medium parameter retrieval based on equation (17).
To avoid confusion, we call the medium parameter ζ when considering (17) in spatial coordinates.In other words, we start by considering 0)v(0, s).A staggered finite-difference grid with primary grid steps δ j and dual grid steps δj has the grid points and a two-point finite-difference discretisation of the first-order system (18) gives where the notation v j = v(x j , s) and vj = v(x j , s) is used.Furthermore, ζ j and ζj denote averaged medium parameters at the grid locations x i and xi , respectively.After symmetrization with diagonal matrices the finite-difference pencil that corresponds to (18) can be written in terms of the transpose-symmetric tridiagonal matrix In ROM embedding we try to interpret the ROM constructed from the measurement data as a discretisation of the underlying equation.If we compare the above tridiagonal discretisation stencil to the tridiagonal ROM matrix from equation (10) we find that the discretisation has twice as many unknowns as the ROM has parameters, i.e. we can't disentangle the grid steps from the local medium parameters.
In [1] it was shown that a tridiagonal ROM that matches the lowest 2n poles and residues of the transfer function corresponds to a discretisation on a special grid, also known as the optimal grid or spectrally matched grid.This grid is independent of the medium parameter ζ in the asymptotic limit n → ∞ and can be computed from the ROM of a reference simulation with ζ 0 (x) = 1.Let this reference grid be characterised by the primary and dual weights γ 0 j and γ0 j , respectively, then pointwise estimates of ζ can be directly extracted from the ROM.To be more specific, let ζ ROM (x) and ζROM (x) be interpolants with interpolation properties with x optimal 0 = 0, then it can be shown that ζ ROM and ζROM converge pointwise in L 1 (0, L) to the true medium profile ζ as n → ∞, see [1] for details.
Having discussed equation (18) using the results of [1], let us now include kinematic wave effects and consider equation ( 17) to obtain an optimal grid reconstruction scheme for wave propagation problems.Initially following the same procedure as above, we can show that the ratios γ j /γ 0 j and γ0 j /γ j converge to the wave speed c[x(T )] parametrised in slowness coordinates and the optimal grids in equations ( 22) and (23) are primary and dual optimal slowness grids T opt j and T opt j , respectively.To obtain the wave speed in physical coordinates, the inverse slowness transform x(T ) : T → x given by needs to be extracted as well.This can be realized in two steps: First the optimal grid is adjusted to the average slowness of the medium and, second, the grid is locally adjusted to slowness coordinates.
The background grid should be computed for the domain [0, T (L)], however this requires knowledge of the average slowness as this defines T (L) = Lc −1 .There are many ways to extract the average slowness from the ROM and we choose to extracted it from the first ROM coefficient γ1 as for which obviously the wave speed c(0) at the sensor location is required.The effectiveness of the above equation is due to the fact that γ0 1 depends linearly on the domain size and in the limit n → ∞ the ratio γ0 1 γ1 c −1 converges to the wave speed at x = 0.Alternatively, T (L) can be extracted from the limit lim s→0 1 s f (s).With this in place the wave speed c(x) can be extracted from the ROM parameters as which are to be interpreted as pointwise estimates of the wave speed at the optimal grid points that are adjusted to the local slowness coordinates and These estimates converge pointwise in L 1 (0, T (L)) to the true wave speed which can be obtained as a corollary to Theorem 6.1 in [1].Essentially, the optimal grid embedding recovers c[x(T )] at points on the optimal grid in slowness coordinates.This recovered c[x(T )] then provides the inverse slowness transform transform T → x, to embed c(x) into physical space.

Krein embedding
From regular Sturm-Liouville theory it is well known that equation ( 1) is satisfied for infinitely many eigenpairs (ϕ i (x), λ i ), with real eigenfunctions ϕ i (x) and imaginary eigenvalues s = λ i .The true impedance function is thus a meromorphic function with an infinite number of poles, corresponding to the eigenvalues λ i and λi .Now Krein and optimal grid embedding approaches both utilise truncated spectral impedance data, and in [5] it was shown by Kac and Krein that there is a one-to-one correspondence between an n-term truncated spectral impedance function and a medium with a mass function M (x) with n points of increase.Krein embedding allows us to link this mass function to a continuous wave speed profile c(x).
The following proposition provides an explicit connection between the continuous wave equation and the discrete finite-difference problem obtained from truncated spectral impedance data.
Proposition 1 [Kac and Krein] Let w n satisfy the wave equation with the mass function M n in the weak sense, where The continuous function w n then interpolates the finite-difference approximation from equation (6), that is, we have w n (x i ) = u i for i = 0, 1, ..., n.

Proof:
The proof is straightforwardly obtained by substituting a (piecewise) linear interpolation of u i between the x i 's into the second-order finite-difference equation.More precisely, let ũ(x) be the linear interpolation of u i on the grid {x i }.For ũ(x) we have Substituting this into the second-order form of the finite-difference relation from equation ( 6) yields and after integration, we obtain the weak form of equation (32) x i +δ In Krein embedding we interpret the ROM as the impedance function at x = 0 of a finite-difference discretization of the first order system (15).In particular, a two-point Solving inverse scattering problems via reduced-order model embedding procedures 12 finite-difference discretization on a staggered grid with primary stepsizes δ j and dual stepsizes δj reads ûj+1 − ûj where we used the shorthand notation u j = u(x j , s) again.The primary ROM edgeweights γ j are therefore interpreted as a step size δ j and the dual edge-weights γj as c −2 j δj .To be more precise, we introduce the nondecreasing mass function and the coordinate transform T → x given by The mass function can now be written as and we observe that for j = 1, 2, ..., n can be interpreted as ROM quadrature rules of ( 41) and (42), respectively.Finally, if we let T Krein j denote the slowness coordinate corresponding to x Krein j then, using the results of the optimal grid case discussed in the previous section, we have and Convergence of the Krein embedding follows directly from equations ( 44) and (45).Since the optimal grid parameters converge, convergence of χ n (T ) → x(T ) and M n (x) → M (x) follows again as a corollary of Theorem 6.1 in [1].

Numerical examples
Optimal grid embedding leads to pointwise estimates for the wave speed c(x), whereas Krein embedding leads to an estimation of the mass function M (x) = x 0 c −2 (ξ)dξ.Therefore, it is natural to display the inversion result of the optimal grid embedding and Krein embedding in terms of these two quantities.In Figure 2 the dashed line signifies a smoothly varying velocity profile on the interval [0, 2] that we attempt to reconstruct using knowledge of poles and residues of the impedance function at x = 0.This spectral data is obtained by applying the vectorfit algorithm [3] to the impedance function at x = 0 to obtain a rational function representation in pole-residue form as indicated in equation ( 2).
We start with n = 12 pole-residue pairs of the impedance function f (s) to construct a ROM and realise optimal grid and Krein embeddings.In this example, the background grid in case of optimal grid embedding is extracted from a reference simulation with a constant c 0 = 1.With only n = 12 spectral points, the wave speed c(x) in case of an optimal grid embedding and the mass function M (x) in case of Krein embedding are accurately reconstructed as illustrated in Figures 2 and 3, respectively.We note that due to the low number of pole-residue pairs that are used, the estimates are slightly misplaced as the map T → x is estimated with only a few quadrature points.
Therefore, let us increase the number of pole-residue pairs to n = 50.In this case we show the performance of optimal grid and Krein embedding in terms of velocity profile and mass function reconstruction.Figure 4 shows the optimal grid embedding reconstruction for n = 50, and clearly indicates convergence of the optimal grid embedding approach.Figure 5 shows the velocity profile reconstruction with the Krein embedding approach, which indicates that Krein embedding converges as   well.Finally, mass function reconstructions for optimal grid and Krein embedding are shown in Figures 6 and 7, respectively.These reconstruction results indicate that both embedding approaches converge and benefit if more spectral data (pole-residue pairs) of the impedance function is available.
Finally, we show how the two embedding approaches reconstruct a smoothed step (piecewise constant) velocity profile.The dashed line Figure 8 the velocity profile along with the optimal grid embedding reconstruction for n = 50.Clearly, optimal grid embedding captures the exact wave speed away from the step and shows some Gibbs ringing around the step.The Krein embedding reconstruction for n = 50 of the corresponding mass function is shown in Figure 9.Here Gibbs ringing is not observed, since an integrated velocity profile is reconstructed.

Embedding of ROMs on semi-infinite domains
Up till now, we have considered scattering problems on bounded domains.For semiinfinite domains, that is, domains that extend to infinity in one direction, the optimal and Krein grids may be applied, but these approaches produce excessively large step sizes [4] as illustrated in Figure 10.This problem can be avoided, however, by constructing ROMs that can be interpreted as two-point finite-difference discretisations of dissipative first-order continuous wave equations.We refer to this approach as Krein-Nudelman embedding [6].
To be specific, consider the wave equation in the Laplace domain on a semi-infinite     x 0 c(ξ) −2 dξ from impedance measurements f (s) = u(0, s).The impedance function measured at x = 0 for an open domain does not have a point spectrum.To implement Krein-Nudelman embedding, we fit the impedance function with rational functions with n complex-conjugate pole-residue pairs.Since the impedance function is passive, the poles and residues have a positive real part.In the presented results, we use vector fit [3] to obtain the poles and residues from the impedance function.
With complex poles and residues the diagonal entries (called α in Algorithm 1) of the tridiagonal matrix obtained from the Lanczos algorithm are no longer zero.Since similar matrices have the same trace we can equate the sum of the eigenvalues in Λ to the sum of the diagonal elements of T TS , called α i 's in the Lanczos algorithm and see that the diagonal entries cannot all vanish.These diagonal elements are interpreted as a loss term in Krein-Nudelman embedding and absorb energy.To facilitate an embedding that allows for absorption, we introduce the complex symmetry tridiagonal matrix T TS;open which we obtain from running the Lanczos algorithm with the complex poles and residues obtained from fitting the impedance function of a semi-infinite domain.To interpret this tridiagonal matrix using Krein-Nudelman embedding consider a finite-difference discretization of equation ( 46) on the grid displayed in Figure 11 ûj+1 − ûj and compare it to finite difference discretization introduced by the tridiagonal matrix in equation (48) ûj+1 − ûj γj+1 + (α 2j+1 + s)u j = 0, ∀j = 0, . . ., n − 1, Note that the Sommerfeld condition is applied on the last dual grid point, whereas the Dirichlet condition in the bounded domain was applied on the last primary grid node.The diagonal of T TS;open appears as a loss terms in the equation, and can thus be used to define integrated primary and dual loss coefficients A Sommerfeld radiation condition would be equivalent to all r i and ri values vanishing except for rn .
To illustrate this embedding we consider the inverse problem of recovering the velocity profile in Figure 12 from impedance measurements at x = 0. We use n = 121 pole-residue pairs to fit the impedance function and compute T TS;open using the complex symmetric Lanczos algorithm.
The Krein-Nudelman embedding of the obtained mass function is shown alongside the true mass function in Figure 13.The smooth velocity variations are recovered well, whereas the reflector is not recovered.In Figure 14 the loss coefficients are embedded into space for visualisation purposes.The primary losses are displayed on the Krein grid and the dual losses in between Krein grid nodes.Almost all losses r/r vanish up until the last reflector in the medium, where the last thee ri contribute to an effective loss term.The recovered r and r terms do not exactly correspond to a Sommerfeld radiation condition, but rather an effective absorbing medium being placed at the end of the Krein embedding.Nonetheless, the velocity profile up to the last reflector is well recovered at which point the loss coefficients increase.Note that the last grid point of the Krein embedding is outside the interval of wave speed variations, since the Sommerfeld condition is prescribed on a dual grid node in between the Krein-grid nodes.
Here we reach a fundamental limitation of the Krein embedding.From the fitted pole residue form, the Krein embedding cannot distinguish between an infinite lossless medium, which does not have a point spectrum, and a finite-absorbing medium.See the Appendix for a further discussion.
However, up until the Krein embedding recovers an "effective" absorbing medium the embedding can be used to solve inverse problems.This is in some sense physical as one typically does not know from which point the infinite domain has constant medium parameters.Compared to an optimal grid approach, the Krein-Nudelman approach only produces a grid in the interval where the variations of the wave speed are supported, whereas the optimal grid approach places grid points beyond the variations of the medium, which is not useful for inversion.
characterised by σ = 0 everywhere, while if σ > 0 in some domain in space, we refer to the above system as a lossy first-order system.
Eliminating v in the lossless case, we obtain the Hermitian system Aũ = λũ + sg with A = LL * and λ = s 2 .Here, the transfer functions of the Hermitian and first-order system are related by f (λ) = λ −1/2 h(λ 1/2 ) and h satisfies the three criteria for a passive SISO system, that is, we have (i) h(s) has no poles in C + , (ii) h(s) = h(s) for every s ∈ C + , where the overbar denotes complex conjugation, and (iii) Re[h(iω + 0)] ≥ 0 for ω ∈ R.
For a lossy first-order system, the transfer function h(s) still satisfies the above three criteria, but we no longer arrive at a Hermitian form when v is eliminated from the first-order system.However, we can establish an isomorphism between passive transfer functions and Stieltjes functions.More precisely, we have Proposition 2. Using f (λ) = λ −1/2 h(λ 1/2 ), every Stieltjes function f (λ) can be transformed to a passive transfer function h(λ 1/2 ) and vice versa.
Proof: Consider the transformation f → h.First note that from the definition (52) of a Stieltjes function it follows that f (λ) is an analytic function of λ on C \ R − .Setting s = λ 1/2 , it follows from h(s) = sf (s 2 ) and the analyticity of f that h(s) is analytic in C + and therefore h(s) has no poles in C + .Furthermore, since ρ is real (positive), we have f (s 2 ) = f (s 2 ) and it follows that h(s) = h(s) with s ∈ C + .Finally, we have where we have used h(s) = sf (s 2 ) and Reasoning in reverse order gives the "vice versa" result.□ Now to each Stieltjes function there corresponds a unique Stieltjes-Krein string [5] and a consequence of the above proposition is that we cannot distinguish between lossy media and lossless media on a semi-infinite domain, since in both cases the transfer function is passive.

Figure 1 :
Figure1: Grid used to interpret the ROM as a finite-difference discretisation of the underlying differential operator.The crosses represent primary grid nodes and the circles dual grid nodes.On the grey nodes the boundary conditions are applied.

Figure 2 :
Figure 2: Velocity profile c(x) (solid line) and optimal grid reconstruction of this profile based on n = 12 pole-residue pairs of the impedance function at x = 0. Crosses: reconstructed velocity values at primary nodes.Circles: reconstructed velocity values at dual nodes.

Figure 3 :
Figure 3: Mass function M (x) (solid grey line) and Krein embedding reconstruction of this function (solid black line) based on n = 12 pole-residue pairs of the impedance function at x = 0.The Krein grid is visualized with small dots on the x-axis.

Figure 4 :
Figure 4: Optimal grid embedding reconstruction of the velocity profile (solid line) based on n = 50 pole-residue pairs of the impedance function at x = 0. Crosses: reconstructed velocity values at primary nodes, circles: reconstructed velocity values at dual nodes.

Figure 5 :
Figure 5: Krein embedding reconstruction (circles) of the velocity profile (solid line) based on n = 50 pole-residue pairs of the impedance function at x = 0.

Figure 6 :
Figure 6: Optimal grid embedding reconstruction of the mass function (solid line) based on n = 50 pole-residue pairs of the impedance function at x = 0. Crosses: reconstructed mass function values at primary nodes, circles: reconstructed mass function values at dual nodes.

Figure 7 :
Figure 7: Krein embedding reconstruction (black line) of the mass function (grey line) based on n = 50 pole-residue pairs of the impedance function at x = 0.The nodes of the Krein grid are shown as solid dots on the x-axis.

Figure 8 :
Figure 8: Optimal grid embedding reconstruction of a step velocity profile (solid line) based on n = 50 pole-residue pairs of the impedance function at x = 0. Crosses: reconstructed velocity values at primary nodes, circles: reconstructed velocity values at dual nodes.

Figure 9 :
Figure 9: Krein embedding reconstruction (black line) of the mass function (grey line) based on n = 50 pole-residue pairs of the impedance function at x = 0.The nodes of the Krein grid are shown as solid dots on the x-axis.

Figure 10 :
Figure10: Illustration of primary (crosses) and dual (circles) grid nodes of an optimal/Krein grid for a semi-infinite-domain.The grid steps become progressively large as we move away from the measurement point, which coincides with the left-most primary grid node.

Figure 11 :
Figure 11: Grid used to interpret the ROM as a finite-difference discretisation of the underlying differential operator on a semi-infinite domain.The crosses represent primary grid nodes and the circles dual grid nodes.Note that the open boundary condition is applied at the last dual node x = L.

Figure 12 :
Figure 12: Wave speed model for the semi-infinite domain example.A reflector is placed at the end of the medium.

Figure 13 :
Figure 13: Mass function M (x) (solid black line) and Krein embedding grid reconstruction of this function (grey solid line) based on n = 121 pole-residue pairs approximating the impedance function at x = 0.

Figure 14 :
Figure 14: Embedding: Primary r i (×) losses embedded on the Krein grid χ Krein and dual losses ri (•) embedded centered between Krein grid nodes.An effective absorbing inclusion is placed at the end of the domain by the Krein embedding.