{'fraction_non_alphanumeric': 0.04661240590016389, 'fraction_numerical': 0.017917164365676825, 'mean_word_length': 4.082592121982211, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Functional-Input Gaussian Processes with Applications to Inverse Scattering Problems

Surrogate modeling based on Gaussian processes (GPs) has received increasing attention in the analysis of complex problems in science and engineering. Despite extensive studies on GP modeling, the developments for functional inputs are scarce. Motivated by an inverse scattering problem in which functional inputs representing the support and material properties of the scatterer are involved in the partial differential equations, a new class of kernel functions for functional inputs is introduced for GPs. Based on the proposed GP models, the asymptotic convergence properties of the resulting mean squared prediction errors are derived and the finite sample performance is demonstrated by numerical examples. In the application to inverse scattering, a surrogate model is constructed with functional inputs, which is crucial to recover the reflective index of an inhomogeneous isotropic scattering region of interest for a given far-field pattern.


Introduction
Computer experiments, the studies of real systems using mathematical models such as partial differential equations, have received increasing attention in science and engineering for the analysis of complex problems. Typically, computer experiments require a great deal of time and computing. Therefore, based on a finite sample of computer experiments, it is crucial to build a surrogate for the actual mathematical models and use the surrogate for prediction, inference, and optimization. The Gaussian process (GP) model, also called kriging, is a widely used surrogate model due to its flexibility, interpolating property, and the capability of uncertainty quantification through the predictive distribution. More discussions on computer experiments and surrogate modeling using GP models can be found in Santner et al. (2018) and Gramacy (2020). This paper is motivated by an inverse scattering problem in computer experiments, where the computer experiments involve functional inputs and therefore the analysis and inference rely on a surrogate model that can take into account functional inputs. Figure 1 illustrates the idea of inverse scattering. Let the functional input g represent the material properties of an inhomogeneous isotropic scattering region of interest shown in the middle of Figure 1. For a given functional input, the far-field pattern, u s , is obtained by solving partial differential equations (Cakoni et al., 2016) which is computationally intensive. Given a new far-field pattern, the goal of inverse scattering is to recover the functional input using a surrogate model. Therefore, a crucial step to address this problem is to develop a surrogate model applicable to functional inputs. Problems with functional inputs are frequently found in engineering applications of non-destructive testing where measurements on the surface or exterior of an object is used to infer the interior structure. Similar problems also come up in electrical impedance tomography where one wishes to recover the functional input representing the electric conductivity, from the measured current to voltage mapping; see, e.g., Mueller and Siltanen (2020), for the electrical impedance tomography model and Cakoni et al. (2016) for the model studied here for inverse scattering. Despite extensive studies on surrogate modeling using GPs (Gramacy, 2020), the developments for functional inputs are scarce. To the best of our knowledge, most of the existing development on GPs are restrictive to specific applications. For example, Nguyen and Peraire (2015) propose a functional-input GP with bilinear covariance operators and apply it to linear partial differential equations. Morris (2012) develops a kriging model with a covariance function specifically for time-series data. Chen et al. (2021) propose a spectraldistance correlation function and apply it to 3D printing. However, a systematical study of surrogate modeling by GPs, including a general framework with theoretical justifications, is still lacking.
A conventional approach to model functional inputs is to incorporate the idea of basis expansions. There are extensive studies on basis expansion in the functional data analysis literature to approximate functional data (Fan and Zhang, 1999;James, 2002;Malfait and Ramsay, 2003;Hall and Horowitz, 2007;Wang et al., 2010;Chen et al., 2011;Wang et al., 2016). For surrogate modeling, there are some recent studies along this direction where functional inputs are approximated by basis expansion approaches (Tan, 2019;Li and Tan, 2022). After dimension reduction of functional inputs through basis truncation, GP surrogates can be applied. Although ideas along this line are intuitive and easy to implement, there are three drawbacks. First, basis expansion requires an explicit specification of basis functions. Second, basis expansion approximates the functional input and achieves dimension reduction by a finite truncation of the basis functions, which can introduce additional bias to the model. Third, scaling up the techniques developed by basis expansion to high dimensional functional inputs is challenging due to the curse of dimensionality.
To tackle these problems with functional inputs, a new GP surrogate is proposed by introducing a new class of kernel functions that are directly defined on a functional space.
It is shown that the proposed kernels are closely connected to the idea of basis expansion without the need of specifying individual basis and without the loss of efficiency due to finite truncation. The procedure is general and easy to scale up. Numerical comparisons with the idea of basis expansion for functional inputs are conducted in the simulation studies as well as in the application of inverse scattering problem. The empirical results of the proposed surrogate model appear to outperform those based on basis expansion in terms of prediction accuracy and uncertainty quantification.
Although the proposed surrogate models extend the conventional GPs to functional inputs, the theoretical results, including the convergence rates of the mean squared prediction errors (MSPE) and the connections to experimental design, are nontrivial extensions.
Directly defining the kernels on a functional space reduces the model bias as compared to basis expansion, but posts technical challenges to the theoretical derivations. Additional scattered data approximation techniques, such as the local polynomial reproduction , have to be rigorously applied to the study of convergence rates. The convergence rates are further explored by the notion of fill distances, which provides a concrete connection between the performance of the proposed model and the experimental design in a functional space.
The remainder of the paper is organized as follows. In Section 2, a functional-input GP model is introduced. A new class of kernel functions, including a linear and a nonlinear kernel, and their theoretical properties are discussed in Section 3. Numerical analysis is conducted in Section 4 to examine the prediction accuracy of the proposed models.
In Section 5, the proposed framework is applied to construct a surrogate model for an inverse scattering problem. Concluding remarks are given in Section 6. Detailed theoretical proofs and the R code for reproducing the numerical results are provided in Supplementary Materials.

Functional-Input Gaussian Process
Suppose that V is a functional space consisting of functions defined on a compact and convex region Ω ⊆ R d , and all functions g ∈ V are continuous on Ω, i.e., V ⊂ C(Ω). A functional-input GP, f : V → R, is denoted by where µ is an unknown mean and K(g, g ) is a semi-positive kernel function for g, g ∈ V .
A new class of kernel K(g, g ) for functional inputs is discussed in Section 3.
Given a properly defined kernel function, the estimation and prediction procedures are similar to the conventional GP. Assume that there are n realizations from the functional-input GP, where g 1 , . . . , g n are the inputs and f (g 1 ), . . . , f (g n ) are the outputs. We have f (g 1 ), . . . , f (g n ) following a multivariate normal distribution, N n (µ n , K n ), with mean µ n = µ1 n and covariance K n , where 1 n is a size-n all-ones vector and (K n ) j,k = K(g j , g k ). The unknown parameters, including µ and the hyperparameters associated with the kernel function, can be estimated by likelihood-based approaches or Bayesian approaches. We refer the details of the estimation methods to Santner et al. (2018) and Gramacy (2020).
Suppose g ∈ V is an untried new function. By the property of the conditional multivariate normal distribution, the corresponding output f (g) follows a normal distribution with the mean and variance, where y n = (y 1 , ..., y n ), y i = f (g i ), and k n (g) = (K(g, g 1 ), ..., K(g, g n )) T . The conditional mean of (2) is used to predict f (g), and the conditional variance of (3) can be used to quantify the prediction uncertainty.

A New Class of Kernel Functions
In this section, a new class of kernel functions for the functional-input GP is introduced.
Based on the proposed models, the asymptotic convergence properties of the resulting mean squared prediction errors are derived. Section 3.1 focuses on the discussions of a linear kernel and Section 3.2 extends the discussions to a nonlinear kernel. A practical guidance on the selection of optimal kernel is discussed in Section 3.3. For notational simplicity, the mean in (1) is assumed to be zero in this section but the results can be easily extended to non-zero cases.

Linear kernel for functional inputs
We first introduce a linear kernel for functional inputs g 1 and g 2 : where g 1 , g 2 ∈ V and Ψ is a positive definite function defined on Ω × Ω. It can be shown in the following proposition that this kernel function is semi-positive definite.
Proposition 3.1. The linear kernel K defined in (4) is semi-positive definite on V × V .
By Mercer's theorem (Rahman, 2007), we have where x, x ∈ Ω, and λ 1 ≥ λ 2 ≥ ... > 0 and {φ k } k∈N are the eigenvalues and the orthonormal basis in L 2 (Ω), respectively. Given the positive definite function Ψ, we can construct a GP via the Karhunen-Loève expansion: where Z j 's are independent standard normal random variables, and φ j , g L 2 (Ω) is the inner product of φ j and g, which is φ j , g L 2 (Ω) = Ω φ j (x)g(x)dx. It can be shown that the covariance function of the constructed GP in (6) is K(g 1 , g 2 ) defined in (4), i.e., for any g 1 , g 2 ∈ V .
The proposed surrogate model is equivalent to a basis expansion through the Karhunen-Loève expansion in (6), but it is worth noting that the proposed method only requires the specification of kernel function in (4) instead of an explicit specification of individual basis φ j . Furthermore, there is no dimension reduction or approximation applied to the functional input, and thus there is no additional bias introduced to the surrogate. These advantages as compared to basis expansion are commonly seen in kernel-based methods, such as support-vector machines (SVM), kernel principal components analysis (KPCA), and kernel ridge regression (KRR) (Hastie et al., 2009).
Proof of Proposition 3.2. For any a, b ∈ R and g 1 , g 2 ∈ V , it follows that which finishes the proof.
The proposed kernel function has an intuitive interpretation that connects to Bayesian modeling. In a Bayesian linear regression, the conditional mean function is assumed to be f (x) = x T w, where w is typically assumed to have a multivariate normal prior, i.e., w ∼ N (0, Σ w ). Hence, for any two points, x and x , the covariance of f (x) and f (x ) is Cov(f (x ), f (x )) = x T Σ w x , which can be interpreted as a weighted inner product of x and x . The proposed model (6) can be viewed as an analogy to the Bayesian linear model with functional inputs and the covariance (7) can be viewed as a weighted inner product of the two functions g 1 and g 2 .
To understand the performance of the proposed predictor of (2) with the kernel function defined in (4), we first characterize the mean squared prediction error by the following theorem. Denote the reproducing kernel Hilbert space (RKHS) associated with a kernel Ψ as N Ψ (Ω).
By Proposition 10.28 in , it is can be shown that T g ∈ N Ψ (Ω) and therefore the right-hand side of (8) exists. According to Theorem 3.3, the MSPE can be represented as the distance between T g and its projection on the linear space spanned by {T g 1 , ..., T g n }. This distance can be reduced if g j 's can be designed so that the spanned space can well approximate the space V . We highlight some designs of g j 's in the following two corollaries where the convergence rates of MSPE can be explicitly discussed.
In the following corollaries, the kernel function Ψ is assumed to be Matérn (Stein, 1999): where Θ is a scale parameter, which is a d × d positive diagonal matrix, · 2 denotes the Euclidean norm, σ 2 is a positive scalar, B ν is the modified Bessel function of the second kind, and ν represents a smoothness parameter. Without loss of generality, we assume that Θ is an identity matrix and σ 2 = 1 for the theoretical developments in this section. More detailed discussions of these parameters, including Θ, σ 2 and ν, are given in Section 4.
For g ∈ V ⊂ L 2 (Ω), there exists a constant C 1 > 0 such that Furthermore, if g ∈ N Ψ (Ω), then there exists a constant C 2 > 0 such that Corollary 3.4 shows that if we can design the input functions to be eigenfunctions of Ψ, the convergence rate of MSPE is polynomial. If the functional space is further assumed to be the RKHS associated with the kernel Ψ (i.e. g ∈ N Ψ (Ω)), which is smaller than L 2 (Ω), the convergence rate becomes faster as indicated by (12). This result indicates a significant difference between the proposed GP defined on a functional space and the conventional one defined on a Euclidean space. That is, the convergence results on the right hand side of (11) and (12) depend on the norm of the functional space that the input g lies in, which is different from the conventional GP where only the Euclidean norm is involved.
Instead of selecting the input functions to be eigenfunctions, an alternative is to design the input functions through a set of knots in Ω, i.e., X n ≡ {x 1 , . . . , x n }, where x j ∈ Ω for j = 1, . . . , n, and the convergence rate is derived in the following corollary.
(2) For a quasi-uniform design, there exists a positive constant C such that h Xn,Ω ≤ Cn −1/d Müller, 2009). Therefore, there exists a constant C 4 > 0 such that As compared with the results in Corollary 3.4, the convergence rate of quasi-uniform designs as shown in (13) is slower than the choice of eigenfunctions as shown in (12). Despite a slower rate of convergence, designing functional inputs by Ψ(x j , ·) with space-filling x j 's can be relatively easier in practice than finding eigenfunctions of Ψ.
The proposed linear kernel can be naturally modified to accommodate the potential nonlinearity in f by enlarging the feature space using a pre-sepecified nonlinear transformation M on g, i.e., M : V → V 1 , where V 1 is a function class. The resulting kernel function can be written as and the corresponding GP can be constructed by The convergence results of MSPE can be extended to (14). There are many possible ways to define M so that the feature space can be enlarged; however, the flexibility comes with a higher estimation complexity. In the next section, we propose an alternative to address the non-linearity through a kernel function, which is computationally more efficient.

Nonlinear kernel for functional inputs
In this section, we introduce a new type of kernel function for functional inputs that takes into account the non-linearity via a radial basis function. Let ψ(r) : R + → R be a radial basis function whose corresponding kernel in R d is strictly positive definite for any d ≥ 1.
Note that the radial basis function of (10), whose corresponding kernel is Matérn, satisfies where · L 2 (Ω) is the L 2 -norm of a function, defined by g L 2 (Ω) = ( g, g L 2 (Ω) ) 1/2 , and γ is a parameter that controls the decay of the kernel function with respect to the L 2 -norm.
Proposition 3.6. The function K defined in (15) is positive definite on V × V .
Assume that there exists a probability measure P on V such that V g(t) 2 dP (g) < ∞, for t ∈ Ω (Ritter, 2007). Based on the positive definite function K as in (15), we can construct a GP via the Karhunen-Loève expansion: where ϕ j 's are the orthonormal basis obtained by a generalized version of Mercer's theorem, (Steinwart and Scovel, 2012) with respect to the probability measure P .
The nonlinear kernel in (15) can be viewed as a basis expansion of the functional input based on the fact that are orthonormal basis functions in L 2 (Ω). By a finite truncation of basis expansion, the input g can be approximated by {φ j } M j=1 in R M , for a positive integer M , and therefore f (g) can be approximated by a GP with the correlation function, ψ(γ · 2 ), where · 2 is the Euclidean norm on R M . However, similar to the discussions in Section 3.1, this introduces additional bias due to the finite truncation and requires an explicit specification of the orthonormal functions {φ j } M j=1 and M in advance. Instead, the proposed method directly evaluates the correlation on a functional space through the nonlinear kernel without approximation, and the application only requires the selection of a proper kernel function.
It is also worth noting that the L 2 norm in (15) can be replaced by any Hilbert space norm, such as RKHS norm. Therefore, the nonlinear kernel (15) is flexible and can be generalized to any target space of interest in practice. Nevertheless, the L 2 norm can be approximated by numerical integration methods, such as Monte Carlo integration (Caflisch, 1998), which is computationally more efficient compared to, for example, the RKHS norm that requires the computation of inverting an N × N matrix, where N is the size of grid points.
Based on the proposed nonlinear kernel, the convergence rate of the MSPE is studied in the following theorem.
Theorem 3.7. Suppose that Φ is a Matérn kernel function with smoothness ν 1 , and ψ is the radial basis function of (10), whose corresponding kernel is Matérn with smoothness ν.
Let τ = min(ν, 1). For any n > N 0 with a constant N 0 , there exist n input functions such Based on (17), it appears that the convergence rate is slower than the conventional GP where the inputs are defined in the Euclidean space. Although potentially there are rooms to further sharpen this rate, a slower rate of convergence for functional inputs is not surprising because the input space is much larger than the Euclidean space. Note that since the RKHS generated by a Matérn kernel function with smoothness ν 1 is equivalent to the (fractional) Sobolev space H ν 1 +d/2 (Ω) , the assumption of g ∈ N Φ (Ω) in Theorem 3.7 is equivalent to g ∈ H ν 1 +d/2 (Ω).
If Φ is a squared exponential kernel, the corresponding RKHS is within the RKHS generated by a Matérn kernel function with any smoothness ν 2 . Thus, one can choose a large ν 2 > ν 1 and apply Theorem 3.7 to obtain the same convergence rate as in (17) by replacing ν 1 with ν 2 . Therefore, one can conclude that the convergence rate of RKHS generated by a squared exponential kernel is faster than that of the RKHS generated by a Matérn kernel function with a fixed ν 1 .

Selection of kernels
Based on Sections 3.1 and 3.2, the linear kernel of (4) results in a less flexible model leading to a lower prediction variance but higher bias, while the nonlinear kernel of (15) results in a more flexible model leading to a higher variance but lower bias (Hastie et al., 2009).
To find an optimal kernel function that balances the bias-variance trade-off, the idea of cross-validation is adopted that allows us to select the kernel by minimizing the estimated prediction error.
Although cross-validation methods are typically expensive to implement in many situations, the leave-one-out cross-validation (LOOCV) of GP models can be expressed in a closed form, which makes computation less demanding (Zhang and Wang, 2010;Rasmussen and Williams, 2006;Currin et al., 1988). Specifically, denoteỹ i as the prediction mean based on all data except ith observation and y i as the real output of ith observation. One can show that, for a kernel candidate K, which can be either the linear kernel (4) or the nonlinear kernel (15), the LOOCV error is where Λ n is a diagonal matrix with the element (Λ n ) j,j = (K −1 n ) j,j . Thus, between the linear and nonlinear kernels, the optimal one can be selected by minimizing the LOOCV error.

Numerical Study
In Section 4.1, the sample paths of the functional-input GP with different settings of the parameters are explored. In Section 4.2, numerical experiments are conducted to examine the emulation performance of the proposed method.

Sample path 4.1.1 Linear kernel
The focus of this section is to study how the unknown parameters in the proposed linear kernel (4) affect the generated sample paths. We focus on the Matérn kernel function, which has the form of (9). Three types of unknown parameters are studied, including the d positive diagonal elements of the diagonal matrix Θ, denoted by θ, the positive scalar σ 2 , and the smoothness parameter ν.
The sample paths are generated by the input functions g(x) = sin(αx), where x ∈ Ω = [0, 2π] and α ∈ [0, 1]. The value α indicates the frequency of the periodic function and the RKHS norm of g, i.e., g N Ψ (Ω) , increases monotonically with respect to α. As a result, this input function creates an analogy to the sample paths in conventional GP by studying the paths as a function of α with different parameter settings. In Figure 2, the sample paths are demonstrated with different settings of the three types of parameters. The first row illustrates the sample paths with three different settings of ν, given θ = 1 and σ 2 = 1.
It appears that the smoothness of the resulting sample paths is not significantly affected by the setting of ν, which typically controls the smoothness in conventional GPs. This is mainly because, unlike the conventional Matérn kernel function (Stein, 1999), the derivative of (4) with respect to the input g is not directly related to the parameter ν. The middle panels in Figure 2 demonstrate the sample paths with different settings of θ, given ν = 2.5 and σ 2 = 1. It shows that, as θ increases, the number of the local maxima and minima increases, which agrees with the observations in conventional GPs. Lastly, the bottom three panels show the sample paths with different settings of σ 2 , given ν = 2.5 and θ = 1. Similar to conventional GPs, σ 2 controls the amplitude of the resulting sample paths.

Nonlinear kernel
Based on the nonlinear kernel of (15) with the Matérn kernel function defined in (10), the sample paths are studied with respect to different settings of the γ, ν and σ 2 . As shown in Figure 3, the results appear to be consistent with the observations in conventional GPs where ν controls the smoothness of the function, θ controls the number of the local maxima and minima, and σ 2 controls the amplitude of the functions.
Consider eight functional inputs, which are shown in the first row of Table 1, and assume that these eight functional inputs are implemented for each of the synthetic examples and their outputs are given in Table 1.
Both of the linear kernel (4) and the nonlinear kernel (15) are performed for the proposed functional input GP and their LOOCV errors are reported in Table 2. According to Section 3.3, LOOCV is then used to identify the optimal kernel. By minimizing the LOOCV errors, the linear kernel is identified as the optimal choice for the linear synthetic example, f 1 (g), and the nonlinear kernel is identified as the optimal choice for the nonlinear synthetic examples, f 2 (g) and f 3 (g). For the three synthetic simulations, their MSEs are summarized in Table 2. It appears that the optimal kernels selected by LOOCV are consistent with the selections based on minimizing the MSEs, which shows that LOOCV is a reasonable indicator of the optimal kernel when the ground truth is unknown.
As a comparison, we consider a basis-expansion approach discussed in Section 3.2. That is, consider a functional principal component analysis (FPCA) with truncated components Measurements Kernel f 1 (g) = Ω Ω g f 2 (g) = Ω Ω g 3 f 3 (g) = Ω Ω sin ( Table 2: The leave-one-out cross-validation errors (LOOCVs) and the mean squared errors (MSEs) for three testing functions, in which the errors corresponding to the optimal kernel are boldfaced. (Wang et al., 2016): with the leading M eigenfunctions {ψ j (x)} M j=1 and the corresponding coefficients {z ij } given by: where n = 8 in this example. The number of components, M = 3, is chosen that explains 99.46% of variance. We refer more details regarding the expansion to Wang et al. (2016) and Mak et al. (2018). Given the training input-output pairs, {z i , y i } n i=1 where z i = (z i,1 , . . . , z i,M ), a conventional GP (with a Matérn kernel) is adopted to fit the training data.
The test input, {z i } 11 i=9 can be similarly obtained by (20), and their outputs are predicted by the fitted GP.
In addition to MSEs, we further consider two numerical measurements to evaluate the prediction performance and quantify the uncertainty: an average coverage rate of the 95% Measurements Method f 1 (g) = Ω Ω g f 2 (g) = Ω Ω g 2 f 3 (g) = Ω Ω sin(g)  Table 3: Prediction results of the FIGP and basis-expansion approach, including MSEs, average coverage rates of the 95% prediction intervals, and the proper scores from Gneiting and Raftery (2007), in which the values with better performances are boldfaced.
prediction intervals and an average proper scores. Coverage rate is the proportion of the times that the interval contains the true value, and the proper score is the scoring rule from equation 27 in Gneiting and Raftery (2007), which is an overall measure of the accuracy of the combined prediction mean and variance predictions. A larger proper score indicates a better prediction. The results are summarized in Table 3, which shows that the proposed method, FIGP, outperforms the basis-expansion approach in terms of both predictions and uncertainty quantification. The average scores of FIGP are close to the nominal coverage 95%, while Basis yields very poor intervals that contain less than 80% of the true values.
The scores of FIGP are all positive and significantly higher than the scores from Basis which are negative for nonlinear operators.

Applications to Inverse Scattering Problem
In this section, we revisit the inverse scattering problem shown in Figure 4. Let D ⊂ R 2 denote an inhomogeneous isotropic scattering region of interest, and the functional input g, the support of which is D, is related to the refractive index for the region D of the unknown scatterer. Given a set of finite element simulations as the training data, the goal of inverse scattering is to recover the functional input from a given far-field pattern. To achieve this goal, an important task is to construct a surrogate model for functional inputs.
In this study, 10 functional inputs, including 1 + x 1 , 1 − x 1 , 1 + x 1 x 2 , 1 − x 1 x 2 , 1 + x 2 , 1 − x 2 , 1 + x 2 1 , 1 − x 2 1 , 1 + x 2 2 , and 1 − x 2 2 , are conducted in the training set and the corresponding far-field patterns are shown in Figure 4. The pre-processing of using principal component analysis (PCA) is first applied to reduce the dimensionality of the output images. The first three principal components, denoted by u l ∈ R 1024 , l = 1, 2, 3, are shown in Figure 5, which explain more than 99.9% variability of the data. Therefore, given the functional input g i for i = 1, . . . , 10, the output of far-field images can be approximated by 3 After the dimension reduction, the three-dimensional outputs f 1 (g), f 2 (g) and f 3 (g) are assumed to be mutually independent and follow the functional-input GP as the surrogate model. For any untried functional input g ∈ V , by the results of (2) and (3), the far-field pattern can by predicted by the following normal distribution: . . , f l (g n )), k l (g) = (K l (g, g 1 ), . . . , K l (g, g n )) T , (K n,l ) i,j = K l (g i , g j ), and K l is the kernel function with the hyperparameters estimated based on f l .
The linear kernel of (4) and the nonlinear kernel of (15) are both performed. The optimal kernel is selected by comparing the LOOCV errors in predicting the far-field pattern.
The LOOCV error based on the linear kernel is 9.6 × 10 −8 , which is smaller than that of the nonlinear kernel, 3.1 × 10 −5 . Therefore, we apply the linear kernel and examine its prediction performance for the test function, g(x) = 1 − sin(x 2 ). Similarly to Section 4, the proposed method, FIGP, is compared with the basis-expansion approach, Basis. The images of the true far-field patterns and their predictions are illustrated in Figure 6. Compared to the ground truth, the predictions of FIGP can capture the underlying structures reasonably well with some discrepancies appearing on the lower right corner. On the other hand, the predictions of Basis appear to deviate more from the ground truth, especially the lower There are extensive studies in experimental design for conventional GP surrogates, but the study of optimal designs for GPs with functional inputs remains scarce. In this paper, it is shown that space-fillingness is a desirable property in controlling the convergence rate of MSPE. An interesting topic for future research is to explore the construction procedure for efficient space-filling designs with functional inputs. Apart from experimental designs, another important future work is to explore systematic approaches to efficiently identify the functional input given an observed far-field pattern, which is the ultimate goal of inverse scattering problems. Based on the proposed GP surrogate, we will explore a Bayesian inverse framework and a calibration procedure that integrates computer experiments and real observations.

Supplemental Materials Additional supporting materials can be found in Supplemental
Materials, including the theoretical proofs of Propositions 3.1 and 3.6, Theorems 3.3 and 3.7, and Corollaries 3.4 and 3.5, and the data, R code, and workflow for reproducing the results in Sections 4 and 5.

Supplementary Materials for "Functional-Input Gaussian Processes with Applications to Inverse
Scattering Problems" S1 Proof of Proposition 3.1 For any non-zero (α 1 , ..., α n ) and (g 1 , ..., g n ), it follows the quadratic form: and the quadratic form is strictly greater than zero if g 1 , ..., g n are linearly independent.
This finishes the proof.

S2 Proof of Theorem 3.3
In order to prove the theorem, the following lemma is provided, which can be found in Proposition 10.28 of .
Lemma S2.1. Suppose Ψ is a symmetric and positive definite kernel on Ω. Then the integral operator T maps L 2 (Ω) continuously into the reproducing kernel Hilbert space It is the adjoint of the embedding operator of the reproducing kernel Hilbert space Now we are ready to prove Theorem 3.3. For any u n = (u 1 , ..., u n ) T ∈ R n , it follows that where the third equality is by Lemma S2.1.
Sinceû n := K −1 n k n (g) minimizes the MSPE, it also minimizes T g − n j=1 u j T g j 2 N Ψ (Ω) , where K n and k n (g) are as in (2). By taking minimum on both sides of (S2.1), we obtain which finishes the proof.

S3 Proof of Corollary 3.4
In order to prove this corollary, we need the following lemma, which states the asymptotic rates of the eigenvalues of K. Lemma S3.1 is implied by the proof of Lemma 18 of Tuo and Wang (2020). In the rest of the Supplementary Materials, we will use the following notation. For two positive sequences a n and b n , we write a n b n if, for some constants C, C > 0, C ≤ a n /b n ≤ C . For notational simplicity, we will use C, C , C 1 , C 2 , ... to denote the constants, of which the values can change from line to line.
Proof of Corollary 3.4. By (5) where the last equality holds by Theorem 10.29 of . Recall that g k = φ k for k = 1, ..., n. Now take u k = g, φ k L 2 (Ω) for k = 1, ..., n. It follows from Theorem 3.3 and This indicates that the MSPE depends on the tail behavior of the summation ∞ j=1 λ j g, φ j Because Ψ is a Matérn kernel function defined in (9), Lemma S3.1 implies λ j j − 2ν d . Then, an explicit convergence rate can be obtained via This finishes the proof of (11).
Next, we consider the case g ∈ N Ψ (Ω). Theorem 10.29 of  yields that Then an alternative way to bound (S3.4) is by This finishes the proof of (12), and thus the proof of Corollary 3.4.

S4 Proof of Corollary 3.5
The following lemma is used in the proof.
Proof of Corollary 3.5. For any u n = (u 1 , ..., u n ) T , by (S3.3), we have where the last equality holds because φ j 's are orthogonal basis in L 2 (Ω). Therefore, we can take u n = R −1 n g n , where g n = (g(x 1 ), ..., g(x n )) T . Then the term, sup x∈Ω |g(x) − n k=1 u k Ψ(x k , x)|, becomes the prediction error of the radial basis function interpolation, which is well established in the literature .
For the completeness of this proof, we include the proof of an upper bound here. Let v k (x) = (R −1 n r n (x)) k for k = 1, . . . , n. For any x ∈ Ω, the reproducing property implies that where the inequality is by the Cauchy-Schwarz inequality. A bound on Ψ(x, x)−r n (x) T R −1 n r n (x) can be obtained via Lemma S4.1, which gives where h Xn,Ω is the fill distance of X n . Combining (8) with (S4.6), (S4.7) and (S4.8), we obtain that which finishes the proof.

S5 Proof of Proposition 3.6
Without loss of generality, we assume γ = 1. It suffices to show that K n = (K(g j , g k )) jk is positive definite for any g 1 , ..., g n ∈ V , which can be done by showing that there exist n points a 1 , ..., a n ∈ R m with m < ∞ such that a j − a k 2 = g j − g k L 2 (Ω) . Let V n = span({g 1 , ..., g n }), which is the linear space spanned by g 1 , ..., g n . Clearly, V n is a finite dimensional space. Let φ 1 , ..., φ m be an orthogonal basis of V n with respect to L 2 (Ω), with m ≤ n, which can be found via the Gram-Schmidt process. Thus, given this basis, each g j can be written as a jk φ k with a j = (a j1 , ..., a jm ) T ∈ R m . Then it can be verified that a j − a k 2 = g j − g k L 2 (Ω) .

S6 Proof of Theorem 3.7
We first provide a characterization on the function class V . Since Φ is a positive definite function, we can apply Mercer's theorem to Φ and obtain where λ Φ,1 ≥ λ Φ,2 ≥ ... > 0 are the eigenvalues and {φ k } k∈N are the eigenfunctions, and the summation is uniformly and absolutely convergent. Because g ∈ V ⊂ N Φ (Ω), by Theorem 10.29 of , the summation Definition S6.1. A set Ω ⊂ R d is said to satisfy an interior cone condition if there exists an angle α ∈ (0, π/2) and a radius r > 0 such that for every x ∈ Ω, a unit vector ξ(x) exists such that the cone is contained in Ω.
We need the following lemma, which is Theorem 11.8 of ; also see Theorem 3.14 of . Lemma S6.1 ensures the existence of the local polynomial reproduction.
Lemma S6.1 (Local polynomial reproduction). Let l ∈ N 0 and π l (R d ) be the set of dvariate polynomials with absolute degree no more than l. Suppose Ω ⊂ R d is bounded and satisfies an interior cone condition. Define C 1 = 2, C 2 = 16(1 + sin(α)) 2 l 2 3 sin 2 (α) where α and r are defined in Definition S6.1. Then for all X n = {x 1 , . . . , x n } ⊂ Ω with h Xn,Ω ≤ c 0 and every x ∈ Ω there exist numbersũ 1 (x), . . . ,ũ n (x) with The following lemma provides an upper bound of the MSPE usingũ j (x), which can be found from the proof of Theorem 11.9 and (11.6) of .
Let Ω be a compact region satisfying an interior cone condition. Then for h Xn,Ω ≤ h 0 , where p ∈ π l/2 (R).