# Application and Evaluation of Surrogate Models for Radiation Source Search
## 3 Surrogate Models
### Legendre Polynomials  
For polynomial interpolation and regression, we assume that the physical model responses can be expressed as  
\[u(q)=u(q,w)=\sum_{j=1}^{\infty}w_{j}\Psi_{j}(q), \tag{13}\]  
which we approximate to obtain the surrogate model  
\[u^{N}(q,w)=\sum_{j=1}^{N}w_{j}\Psi_{j}(q). \tag{14}\]  
We take \(\Psi_{j}(q)\) to be multivariate Legendre polynomials--i.e., products of univariate Legendre polynomials--which are defined as the solution to the differential equation  
\[\frac{d}{dq}\Big{[}(1-q^{2})\frac{d}{dq}\Psi_{j}(q)\Big{]}+j(j+1)\Psi_{j}(q)=0. \tag{15}\]  
The first three univariate Legendre polynomials are  
\[\Psi_{1}(q) =1, \tag{16}\] \[\Psi_{2}(q) =q,\] (17) \[\Psi_{3}(q) =\frac{1}{2}(3q^{2}-1) \tag{18}\]  
for \(q\in\hat{\Omega}=[-1,1]\). These polynomials are orthogonal on the interval \(\hat{\Omega}=[-1,1]\) with respect to the density \(\rho(q)=\frac{1}{2}\). We note that the roots of the univariate Legendre polynomials provide the Gauss-Legendre training points, \(q^{m},\ m=1,...,M\), that we use to construct the surrogate models. An important feature of polynomial approximations is that their accuracy is improved for functions with more regularity. We have spectral convergence for Legendre polynomial expansions of functions on \([-1,1]\) with error bounds  
\[||u(q,w)-u^{N}(q,w)||_{L^{2}[-1,1]}\leq C_{k}N^{-k}||u(q,w)||_{H^{k}([-1,1])}. \tag{19}\]  
Here, \(C_{k}\geq 0\) is a constant that may depend on \(k\) and \(u\in H^{k}([-1,1])\), where \(H^{k}([-1,1])\) is a Sobolev space of \(L^{2}\) functions with weak derivatives of all orders up to \(k\) in \(L^{2}\). We should observe the convergence rate \(CN^{-1}||u(q,w)||_{1}\), since our model is continuous, but has non-continuous derivatives [19]. We observe a convergence rate of \(rRMSE(u^{N})/rRMSE(u^{N/2})\approx 0.8\), where \(rRMSE(u^{N})\) means the relative root mean squared error defined in Equation (10) of the Legendre-polynomial-based surrogate model employing \(N\) basis functions. This is not the convergence rate of \(rRMSE(u^{N})/rRMSE(u^{N/2})\approx 0.5\) that we expected to observe, since  
\[\frac{CN^{-1}||u(q,w)||_{1}}{C(N/2)^{-1}||u(q,w)||_{1}}=\frac{N^{-1}}{(N/2)^{ -1}}=\frac{1}{2}. \tag{20}\]  
However, this rate is closer to \(0.5\) for certain values of \(N\), as can be seen in the convergence results in \(rRMSE\) of the Legendre polynomials that we compile in Table 1. We also bound \(N\) from above to address overfitting the noise in the physical model response. We consider two methods to solve for the coefficients \(w_{j}\). The first is based on discrete projection, which exploits the orthogonality of the polynomials \(\Psi_{j}(q)\) to interpolate the function of interest at\(u(q^{m}),\ m=1,\ldots,M\); i.e., the physical model simulations defined in Equation (3). Multiplying both sides of Equation (13) by \(\Psi_{j}(q)\) and integrating over \(\hat{\Omega}\), yields  
\[w_{j}=\frac{1}{\gamma_{j}}\int_{\hat{\Omega}}u(q,w)\Psi_{j}(q)\rho(q)dq,\ \ j=1, \ldots,N, \tag{21}\]  
by exploiting the orthogonality relation  
\[\int_{\hat{\Omega}}\Psi_{j}(q)\Psi_{i}(q)\rho(q)dq= \begin{cases}0&\text{for $i \
eq j$}\\ \gamma_{j}&\text{for $i=j$.}\end{cases}  \tag{22}\]  
Here, \(\gamma_{j}=\frac{2}{2j+1}\) is a normalization constant and \(\rho(q)=\frac{1}{2}\) is the associated density over \(\hat{\Omega}=[-1,1]\). We approximate this integral by a quadrature rule, such as Gauss-Legendre quadrature, to obtain the approximate relation  
\[w_{j}\approx\frac{1}{\gamma_{j}}\sum_{r=1}^{R}u(q^{r})\Psi_{j}(q^{r})\rho(q^{r })\omega^{r},\ j=1,\ldots,N \tag{23}\]  
for the coefficients of the surrogate model. Here, \(q^{r}\) are the sampled inputs, \(\omega^{r}\) are the quadrature weights, and \(u(q)=\hat{u}_{i}(q)+B\), where \(\hat{u}_{i}(q)\) solves Equation (3) and we have dropped the dependence on \(w\). We employ Gaussian quadrature, since with only \(R+1\) training points, polynomials of degree less than or equal to \(2R+1\) can be integrated exactly. We set \(R=M\) and \(N=\frac{(p+K)!}{p!K!