\section{Introduction and simulation details}
\label{intro}

In all extensions of the Standard Model the transitions between pseudoscalar (PS) mesons can be parametrized in terms of three form factors, namely $f_+$, $f_0$ and $f_T$. Recently, using the gauge configurations generated by the European Twisted Mass Collaboration (ETMC) with $N_f = 2 + 1 + 1$ dynamical quarks \cite{Baron:2010bv,Baron:2011sf}, we have computed the three form factors relevant for the $D \to \pi(K) \ell \nu$ transitions, namely the vector and scalar form factors $f_+^{D \to \pi(K)}(q^2)$ and $f_0^{D \to \pi(K)}(q^2)$ in Ref.~\cite{Lubicz:2017syv}, and the tensor form factor $f_T^{D \to \pi(K)}(q^2)$ in Ref.~\cite{Lubicz:2018rfs} with $q^2 = (p_D - p_{\pi(K)})^2$. The form factors have been evaluated in the whole range of values of $q^2$ accessible by the experiments, i.e.~from $q^2 = 0$ up to $q^2_{\rm max} = (M_D - M_{\pi(K)})^2$. In this contribution we present the results of an exploratory study toward a lattice determination of the vector and scalar form factors for the $B \to \pi \ell \nu$ decays.

The analysis of Refs.~\cite{Lubicz:2017syv,Lubicz:2018rfs} has highlighted an important issue, which warrants further investigations. We found evidence of a remarkable breaking of Lorentz symmetry due to hypercubic effects for all the form factors and we presented a method to subtract the hypercubic artefacts and to recover the Lorentz-invariant form factors in the continuum limit.
These effects appear to be dependent on the difference between the parent and the child meson masses. We thus expect that they will play an even more important role in the determination of the form factors governing semileptonic $B$-meson decays into lighter mesons. It is therefore crucial to have them under control. For this reason we have started an investigation of the semileptonic transition $D^\prime \to \pi \ell \nu$, where $D^\prime$ is an unphysical heavy-light meson, built with a heavier-than-charm quark.

In our calculations quark momenta are injected on the lattice using non-periodic (twisted) boundary conditions and the matrix elements of the vector  current and the scalar density are determined for many kinematical conditions, in which parent and child mesons are  moving or at rest.

The gauge ensembles and the QCD simulations used in this work are the same adopted in Refs.~\cite{Lubicz:2017syv,Lubicz:2018rfs}, where the reader is referred to for details.
The simulations have been carried out at three different values of the inverse bare lattice coupling $\beta$, at different lattice volumes (with spacial sizes varying between $\simeq 2$ and $\simeq 3$ fm) and for pion masses ranging from $\simeq 210$ to $ \simeq 450$ MeV~\cite{Carrasco:2014cwa}.

\section{Form factors and hypercubic effects}
\label{formfactors}

Let us introduce the local bare vector and scalar operators, $V_\mu = \bar{c}^\prime \gamma_\mu d$ and $S = \bar{c}^\prime d$, where $c^\prime$ is a heavy quark with mass up to twice the physical charm quark mass determined in Ref.~\cite{Carrasco:2014cwa}. Since in our lattice setup we employ maximally twisted fermions, the vector and scalar operators renormalize multiplicatively~\cite{Frezzotti:2003ni}, i.e.~$\widehat{V}_\mu =  {\cal{Z}}_V\cdot V_\mu$ and $\widehat{S} = {\cal{Z}}_P\cdot S$, where $\widehat{V}_\mu$ and $\widehat{S}$ are the renormalized operators with ${\cal Z}_V$ and ${\cal Z}_P$ being the corresponding renormalization constants.

The matrix elements $\braket{\pi(p_\pi) | \widehat{V}_\mu | D^\prime(p_D^\prime)}\equiv\braket{\widehat{V}_\mu}$ and $\braket{\pi(p_\pi) | {S} | D^\prime(p_D^\prime)}\equiv\braket{{S}}$ can be decomposed into the Lorentz-invariant vector and scalar form factors $f_+(q^2)$ and $f_0(q^2)$ as
\begin{eqnarray}
       \label{eq:vector0L}
       \braket{\widehat{V}_0}  & = & (E_{D^\prime} + E_\pi) f_+(q^2) + (E_{D^\prime} - E_\pi) \frac{M_{D^\prime}^2 - M_\pi^2}{q^2} \left[ f_0(q^2) - f_+(q^2) \right] +
                                                       {\cal{O}}(a^2) ~ , \\
       \label{eq:vectoriL}
       \braket{\widehat{V}_{\rm sp}} & = &  \left\{ (p_{D^\prime} + p_\pi) f_+(q^2) + (p_{D^\prime} - p_\pi) \frac{M_{D^\prime}^2 - M_\pi^2}{q^2}
                                                                                \left[ f_0(q^2) - f_+(q^2) \right] \right\} + {\cal{O}}(a^2) ~ , \\
       \label{eq:scalarL}
       \braket{S} & = & \frac{M_{D^\prime}^2 - M_\pi^2}{\mu_h - \mu_{ud}} f_0(q^2) + {\cal{O}}(a^2) ~ ,
\end{eqnarray}
where $\mu_h$ ($\mu_{ud}$) is the heavy(light) quark bare mass, $\widehat{V}_0$ is the temporal component of the vector current, while $\widehat{V}_{\rm sp}$ is the average of the three spacial components of $\widehat{V}_\mu$ (we inject democratically distributed 3-momenta). The vector and scalar matrix elements can be extracted from the large time distance behavior of three ratios, $R_\mu$ ($\mu = 0, \rm{sp}$) and $R_S$, which are given by
\begin{eqnarray}
\label{eq:Rmu}
    R_\mu(t,\vec{p}_{D^\prime}, \vec{p}_\pi) &\equiv& 4\, p_{{D^\prime} \mu}\, p_{\pi \mu}\, \frac{C^{{D^\prime}\pi}_{V_\mu}(t, t^\prime, \vec{p}_{D^\prime}, \vec{p}_\pi) \,
        C^{\pi {D^\prime}}_{V_\mu}(t, t^\prime, \vec{p}_\pi, \vec{p}_{D^\prime})} {C^{\pi \pi}_{V_\mu}(t, t^\prime, \vec{p}_\pi, \vec{p}_\pi) \,
        C^{{D^\prime}{D^\prime}}_{V_\mu}(t, t^\prime, \vec{p}_{D^\prime}, \vec{p}_{D^\prime})} ~ _{\overrightarrow{t\gg a, \, (t^\prime-t)\gg a}} ~\lvert\braket{\widehat{V}_\mu} \rvert^2 \,,\\[2mm]
    \label{eq:RS}
    R_S(t,\vec{p}_{D^\prime}, \vec{p}_\pi) &\equiv& 4\, E_{D^\prime}\, E_\pi\,\frac{C^{{D^\prime}\pi}_S(t, t^\prime,\vec{p}_{D^\prime}, \vec{p}_\pi)\, C^{\pi {D^\prime}}_S(t, t^\prime,\vec{p}_\pi, \vec{p}_{D^\prime})} {\widetilde{C}_2^{D^\prime} \left(t^\prime,\vec{p}_{{D^\prime}} \right)\, \widetilde{C}_2^\pi \left(t^\prime,\vec{p}_{\pi} \right)} \,~ _{\overrightarrow{t\gg a, \, (t^\prime-t)\gg a}}~\lvert\braket{S} \rvert^2  \,.
\end{eqnarray}
In Eqs.~(\ref{eq:Rmu})-(\ref{eq:RS}) $C^{D^\prime\pi}_{\Gamma}$ ($\Gamma=V_\mu,\,S$) and $\widetilde{C}_2^M$ are respectively the 3-point correlation function between the $D^\prime$ and the $\pi$ mesons and the 2-point correlation function for the meson $M$ in which the backward signal has been cancelled (see Ref.~\cite{Lubicz:2017syv}). At large time distances they are defined as
\begin{eqnarray}
\label{eq:C3_larget}
        C^{{D^\prime}\pi}_{\widehat{\Gamma}}\left( t,\, t^\prime,\, \vec{p}_{D^\prime},\, \vec{p}_\pi \right) && ~ _{\overrightarrow{t\gg a\,, \, (t^\prime-t)\gg a}} ~
              Z_\pi Z_{D^\prime}^* \braket{\pi(p_\pi)|\widehat{\Gamma}|{D^\prime}(p_{D^\prime})}\, e^{-E_{D^\prime} t} \, e^{-E_\pi (t^\prime - t)} / (4E_\pi E_{D^\prime}) ~ , \\
\label{eq:C2_tilde_larget}
        \widetilde{C}_2^{M} \left(t,\,\vec{p}_{M}\right) && ~ _{\overrightarrow{t \gg a}} ~ Z_M \, e^{-E_M t} / (2 E_M) ~ ,
\end{eqnarray}
where $E_{D^\prime}$ and $E_{\pi}$ are the energies of the $D^\prime$ and $\pi$ mesons, while $Z_D^\prime$ and $Z_\pi$ are the matrix elements $\braket{0\lvert\,P_5^{D^
\prime}(0)\,\rvert\,D^\prime(\vec{p}_{D^\prime}) }$ and $\braket{0\lvert\,P_5^\pi(0)\,\rvert\,\pi(\vec{p}_{\pi})}$, which depend on the meson momenta $\vec{p}_D^\prime$ and $\vec{p}_
\pi$ because of the use of smeared interpolating fields.
From the 2- and 3-point correlators we are able to extract the matrix elements $\braket{\widehat{V}_\mu}$ and $\braket{S}$, which allow us to determine $f_+(q^2)$ and $f_0(q^2)$ as the best-fit values of Eqs.~(\ref{eq:vector0L})-(\ref{eq:scalarL}).

The momentum dependencies of the vector and scalar form factors are illustrated in Fig.~\ref{fig:fishbone}, where different markers and colors correspond to different values of the child meson momentum for two different ensembles. The unphysical $D^\prime$-meson has a mass of $\simeq 2.7$ GeV corresponding to a heavy-quark mass $m_h \simeq 2.13$ GeV.
If the Lorentz-covariant decomposition~(\ref{eq:vector0L})-(\ref{eq:scalarL}) were adequate to describe the lattice data, the extracted form factors would depend only on the squared 4-momentum transfer $q^2$ (and on meson masses).
This is not the case and an extra dependence on the value of the meson 3-momenta is clearly visible in Fig.~\ref{fig:fishbone}.
\begin{figure}[htb!]
\centering
\includegraphics[width=6.65cm,clip]{fzero_uncorrected_vs_q2_B25_32_0h3.pdf}
\includegraphics[width=6.65cm,clip]{fplus_uncorrected_vs_q2_B75_32_0h3.pdf}
\vspace{-0.25cm}
\caption{\it \footnotesize Momentum dependence of the scalar form factor for the ensemble B25.32 (left panel) and of the vector form factor for the ensemble B75.32 (right panel), corresponding respectively to $M_\pi \approx 250$ MeV and $M_\pi \approx 400$ MeV~\protect\cite{Carrasco:2014cwa}. In both cases the unphysical $D^\prime$-meson has a mass of $\simeq 2.7$ GeV corresponding to a heavy-quark mass $m_h \simeq 2.13$ GeV. Different markers and colors distinguish different values of the pion momentum.}
\label{fig:fishbone}
\end{figure}
The origin of this effect is that the lattice breaks Lorentz symmetry and it is invariant only under discrete hypercubic rotations. Therefore, the form factors determined form Eqs.~(\ref{eq:vector0L})-(\ref{eq:scalarL}) depend also on hypercubic invariants.

The upper panels of Fig.~\ref{fig:elastic_vs_nonelastic} shows how the hypercubic effects on the vector and scalar form factors change with the mass difference ($M_{D^\prime} - M_\pi$). In the lower panels the case of the transition between two nearly degenerate PS mesons close to the physical D-meson is shown. Within the statistical uncertainties the corresponding form factors exhibit no evidence of hypercubic effects.
\begin{figure}[htb!]
\centering
\makebox[\textwidth][c]{
\includegraphics[width=6.65cm,clip]{fzero_uncorrected_various_mh.pdf}
\includegraphics[width=6.65cm,clip]{fplus_uncorrected_various_mh.pdf}
}
\makebox[\textwidth][c]{
\includegraphics[width=6.65cm,clip]{fzero_uncorrected_elastic.pdf}
\includegraphics[width=6.65cm,clip]{fplus_uncorrected_elastic.pdf}
}
\vspace{-0.75cm}
\caption{\it \footnotesize Upper panels: the vector and scalar form factors for a transition between two non-degenerate PS mesons for three different values of the heavy-quark mass. Lower panels: the vector and scalar form factors for transitions between two nearly degenerate PS mesons close to the physical D-meson.}
\label{fig:elastic_vs_nonelastic}
\end{figure}

\section{Global fit}
\label{globalfit}

A possible way to describe hypercubic artifacts affecting the form factors $f_+$ and $f_0$ is to address them directly on the vector and scalar matrix elements (\ref{eq:vector0L})-(\ref{eq:scalarL}).
The analysis closely follows the strategy presented in Ref.~\cite{Lubicz:2017syv} for the vector and scalar $D \to \pi(K)$ form factors. We refer the interested reader to the more detailed discussion presented in Refs.~\cite{Lubicz:2017syv,Lubicz:2018rfs}.
To illustrate the procedure we consider the following decomposition of the vector matrix elements
\begin{equation}
    \braket{\pi(p_\pi) | \widehat{V}_\mu^E | {D^\prime}(p_{D^\prime})} = \braket{\widehat{V}_\mu^E}_{\rm Lor} + \braket{\widehat{V}_\mu^E}_{\rm hyp} ~ ,
    \label{eq:vector_decomposition}
\end{equation}
where $\braket{\widehat{V}_\mu^E}_{\rm Lor}$ is the Lorentz-covariant term
\begin{equation}
    \braket{\widehat{V}_\mu^E}_{\rm Lor} = P_\mu^E~ f_+(q^2, a^2) + q_\mu^E ~ \frac{M_{D^\prime}^2 - M_{\pi}^2}{q^2} \left[ f_0(q^2, a^2) - f_+(q^2, a^2) \right] ~ ,
    \label{eq:vector_Lorentz}
\end{equation}
while $\braket{\widehat{V}_\mu^E}_{\rm hyp}$ is given by
\begin{equation}
    \braket{\widehat{V}_\mu^E}_{\rm hyp} = a^2 \left[ \left( q_\mu^E \right)^3 ~ H_1+ \left( q_\mu^E \right)^2 P_\mu^E ~ H_2 +
                                                                    q_\mu^E \left( P_\mu^E \right)^2 ~ H_3 + \left( P_\mu^E \right)^3 ~ H_4 \right] ~ .
    \label{eq:vector_hypercubic}
\end{equation}
The suffix $E$ emphasizes that all the relations are written in the Euclidean space, e.g.~$q_\mu^E = \left( \vec{q}, ~ q_4 \right) = \left( \vec{q}, ~ iq_0 \right)$ so that $\sum_\mu q_\mu^E q_\mu^E = - q^2$.
The new structures, allowed by the lattice hypercubic symmetry, introduce the new form factors $H_i$ ($i = 1, ..., 4$).
For them we have adopted a polynomial expansion $H_i(z) = d_0^i + d_1^i z + d_2^i z^2$ given in terms of the $z$ variable~\cite{Boyd:1995cf}.
The coefficients $d_{0,1,2}^i$, which depend on the light-quark mass, are free parameters in the fitting procedure.

As for the Lorentz-invariant terms $\braket{\widehat{V}_\mu^E}_{\rm Lor}$ and $\braket{S}_{\rm Lor}$, the form factors $f_{+,0}(q^2, a^2)$ can be parametrized by the modified z-expansion of Ref.~\cite{Bourrely:2008za}, viz.
 \begin{eqnarray}
   \label{eq:z-exp_f+}
    f_+^{{D^\prime} \to \pi}(q^2, a^2) & = & \frac{f^{{D^\prime} \to \pi}(0, a^2) + c_+(a^2)\, (z - z_0)
                                                               \left(1 + \frac{z + z_0}{2} \right)}{1 - q^2 \left( 1 + P_+ a^2 \right) / \left( M_{D^\prime}^2 + \Delta^2 \right)} ~ ,\\
   \label{eq:z-exp_f0}
    f_0^{{D^\prime} \to \pi}(q^2, a^2) & = & \frac{f^{{D^\prime} \to \pi}(0, a^2) + c_0(a^2)\, (z - z_0)
                                                               \left(1 + \frac{z + z_0}{2} \right)}{1 - q^2 / M_S^2} ~ ,
 \end{eqnarray}
where for $c_{+,0}$ we have adopted a linear dependence in $a^2$, $M_{D^\prime}$ is the PS mass calculated on the lattice, $\Delta^2 \equiv M_{D*}^2 - M_D^2$ is taken from experiments and $M_S$ is left as a free parameter.
For the vector form factor at zero 4-momentum transfer, $f^{{D^\prime} \to \pi}(0, a^2)$, we use the following Ansatz
\begin{eqnarray}
    \label{eq:ChLim}
    f^{{D^\prime} \to \pi}(0, a^2)\!\!&=&\!\! F_+ \left[ 1 + - \frac{3}{4} \left( 1 + 3 \widehat{g}^2 \right) \xi_\ell \log\xi_\ell + b_1\, \xi_\ell + b_2 \, \xi^2_\ell + b_3 \, a^2 \right] ~,
\end{eqnarray}
where the coefficients $F_+$, $b_1$, $b_2$ and $b_3$ are treated as free parameters.
In Fig.~\ref{fig:corrected} we show the same form factors as in Fig.~\ref{fig:fishbone} after the hypercubic contributions determined by the global fit have been subtracted. It can be seen that the vector and scalar form factors depend now only on the $4-$momentum transfer $q^2$.
\begin{figure}[htb!]
\centering
\includegraphics[width=6.65cm,clip]{fzero_corrected_vs_q2_B25_32_0h3.pdf}
\includegraphics[width=6.65cm,clip]{fplus_corrected_vs_q2_B75_32_0h3.pdf}
\vspace{-0.25cm}
\caption{\it \footnotesize The same as in Fig.~\ref{fig:fishbone} after subtracting the hypercubic contributions as determined by the global fit.}
\label{fig:corrected}
\end{figure}

The impact of the hypercubic effects at fixed lattice spacing is shown in Fig.~\ref{fig:HMaRcomparison_latt}, where we compare the scalar form factor $f_0(q^2)$ obtained by fitting the data from Eqs.~(\ref{eq:vector0L})-(\ref{eq:scalarL}) corresponding only to the heavy meson at rest with the one obtained using all the available data after subtracting the hypercubic effects determined from the global fit.
\begin{figure}[htb!]
\centering
\makebox[\textwidth][c]{
\includegraphics[width=6.65cm,clip]{fit_band_fzero_vs_HMaR_B55_32_0h2.pdf}
\includegraphics[width=6.65cm,clip]{fit_band_fzero_vs_HMaR_B55_32_0h4.pdf}
}
\makebox[\textwidth][c]{
\includegraphics[width=6.65cm,clip]{fit_band_fzero_vs_HMaR_D20_48_0h2.pdf}
\includegraphics[width=6.65cm,clip]{fit_band_fzero_vs_HMaR_D20_48_0h4.pdf}
}
\vspace{-0.75cm}
\caption{\it \footnotesize Momentum dependence of the scalar form factor $f_0(q^2)$ corresponding to $m_h \simeq 1.85$ GeV (left panels) and $m_h \simeq 2.50$ GeV (right panels) in the case of the gauge ensembles B55.32 and D20.48~\protect\cite{Carrasco:2014cwa}. The filled bands correspond to the global fit where all the available data are included after subtraction of the hypercubic effects. The dashed bands correspond to the fit of the data corresponding only to the heavy meson at rest without taking into account hypercubic corrections. The pion mass and the lattice spacing are respectively equal to $M_\pi \simeq 375$ MeV and $a \simeq 0.0815$ fm for the ensemble B55.32, and $M_\pi \simeq 250$ MeV and $a \simeq 0.0620$ fm in the case of the ensemble D20.48.}
\label{fig:HMaRcomparison_latt}
\end{figure}

We have performed the extrapolation to the physical pion point and to the continuum and infinite volume limits using two different fitting procedures: ~ i) including all available data taking into account the hypercubic corrections (see Eq.~(\ref{eq:vector_hypercubic})), and ~ ii) considering only the data with the heavy meson at rest neglecting hypercubic effects. In Fig.~\ref{fig:HMaRcomparison_cont} the (percentage) difference between the corresponding (correlated) results is shown and it indicates that the subtraction of hypercubic effects is required for obtaining properly the continuum limit.
\begin{figure}[htb!]
\centering
\makebox[\textwidth][c]{
\includegraphics[width=7.10cm,clip]{fplus_diff_q2_h2345_ver3.pdf}
\includegraphics[width=7.10cm,clip]{fzero_diff_q2_h2345_ver3.pdf}
}
\vspace{-0.75cm}
\caption{\it \footnotesize Percentage difference between the form factors extrapolated to the physical point obtained using either all available data (corrected for hypercubic effects) or using only data with the heavy meson at rest (uncorrected for hypercubic effects). Purple, red, green and grey bands correspond respectively to $m_h = 1.84, 2.13, 2.47$ and $2.87$ GeV. The short-dashed lines identify the extrapolation to values of $q^2$ not covered directly by the lattice data.}
\label{fig:HMaRcomparison_cont}
\end{figure}

\section{Conclusions}
\label{conclusions}

In this contribution we have presented the preliminary results of our study of the form factors describing the semileptonic transition between a heavy-light pseudoscalar meson, with masses above the physical D-meson one, and the pion.
Our analysis is based on the gauge configurations produced by ETMC with $N_f = 2 + 1 + 1$ flavors of dynamical quarks at three different values of the lattice spacing and with pion masses as small as $210$ MeV.
We have simulated a series of heavy-quark masses in the range $m_c^{phys} < m_h < 2m_c^{phys}$.

The Lorentz symmetry breaking due to hypercubic effects is clearly observed in the data as already found in Refs.~\cite{Lubicz:2017syv,Lubicz:2018rfs} in the case of the semileptonic form factors of the $D \to \pi(K) \ell \nu$ transitions.  The impact of the hypercubic contributions, included in the decomposition of the current matrix elements in terms of additional form factors, remains significant as the decaying meson mass increases and it should be taken into account for performing properly the continuum limit.

Therefore, we expect that hypercubic effects will play an important role in the lattice determination of the form factors governing the $B \to \pi \ell \nu$ semileptonic decays.
Further investigations in this direction are in progress.

\bibliographystyle{JHEP}


\section*{Bayesian Numerical Methods}

Numerical analysis is concerned with the approximation of typically high or infinite-dimensional mathematical quantities using discretisations of the space on which these are defined. Examples include integrands, and as a consequence their corresponding integrals, or the solutions of differential equations. Different discretisation schemes lead to different numerical algorithms, whose stability and convergence properties need to be carefully assessed. Such numerical methods have undeniably played, and continue to play, an important role in the implementation of statistical methods, from the computation of intractable integrals in Bayesian statistics to the computation of estimators using optimisation routines in frequentist statistics.

In essence, a Bayesian numerical method models the quantity of interest as a latent variable. A prior model is posited on this quantity, then the discretisation of the space is used as data to turn the prior distribution into a posterior distribution. The prior models used are often constructed over functions spaces, and the literature on Bayesian numerical methods therefore draws much inspiration from the work on Bayesian nonparametrics, such as Gaussian processes.

Since their inception in the seventies and eighties, Bayesian numerical methods have initiated much debate as to the foundations of computation and approximation of unknown functions. In this regard, Bayesian numerical methods have initiated fresh considerations of what it means to ``know a function'' \citep{Diaconis1988}. Several important questions have been raised - many of which are highlighted by the discussants of this paper. Are Bayesians simply reinventing the bread and butter of numerical analysis using a different scientific language? How should a prior model on such an abstract quantity be selected? What interpretation should be given to the posterior?

These types of questions have a different emphasis depending on the priorities of the varying scientific communities that are engaged in this debate.  We hope that the present paper and contributed comments can bring some focus to this discussion, and help clarify the advantages and disadvantages of the approach. The answer from the discussants to our first central question can be summarised as ``Yes, but...", which we believe is representative of the overall thoughts of the various research communities represented in this field. There are of course many caveats with these methods, and they are not reasonably expected to supplant numerical analysis, but we will argue that some of their advantages can be significant for certain application areas.

\section*{The Quantification of Uncertainty in Numerical Problems}

One motivation for the use of Bayesian methods in numerical analysis is the possibility of quantifying uncertainty emanating from discretisation error. Indeed, the posterior provides an entire distribution over the solution of the numerical problem which reflects any uncertainty remaining due to the fact that we have access to a finite amount of data, or in other words, that our observations of the quantity of interest are limited to the discretisation of the space. This is particularly appealing since this distribution can be used for several tasks, including to decide whether to refine the discretisation, and if this is of interest, to construct an experimental-design strategy to do so.

However, there are several technical challenges remaining to guarantee that our quantification of uncertainty is valid. First, we must justify our choice of prior model using knowledge of the function derived from domain expertise. This can clearly be a difficult task. Second, we also need to verify the consistency of the method, a problem which, unlike in finite-dimensional spaces, can be technically challenging \citep{Diaconis1986,Owhadi2015Brittleness}. Finally, even though the confidence intervals are valid from a Bayesian viewpoint (in the sense that they reflect our remaining uncertainty given our prior and the data), we may want to verify that the posterior is calibrated from a frequentist viewpoint; once again a non-trivial matter \citep{Cox1993}.

In the paper, we illustrated these issues for an algorithm called Bayesian quadrature/cubature, which tackles numerical integration by modelling the integrand using a prior (usually a Gaussian process), and using $n$ function evaluations to obtain a posterior on this integrand, which itself induces a posterior on the value of the integral. Our paper studied the consistency problem for this algorithm, and provided asymptotic rates. It also studied the calibration of the algorithm on some synthetic problems, where it was demonstrated that the posterior had good frequentist coverage for large $n$ (but no theory for this was provided).

In his comment, Owen argues that Monte Carlo methods are `unreasonably effective' thanks to the confidence intervals provided by the central limit theorem (CLT), which can be estimated at a faster rate in $n$ than the integral of interest. However, we point out that the CLT only provides confidence intervals which are valid in an asymptotic regime where $n \rightarrow \infty$. On the other hand, worst-case error used in the quasi-Monte Carlo (QMC) literature are also valid for finite $n$, but as pointed out by Owen, Hickernell and Jagadeeswaran, these are too pessimistic for most problems of interest. In comparison, contingent on a good methodology for eliciting a prior (a non-trivial task), the Bayesian approach can provide non-pessimistic uncertainty quantification for finite $n$.

Recent developments in the theory of Bayesian nonparametrics also demonstrate that certain Bayesian models can be well-calibrated in a frequentist sense for large classes of target functions \cite{Szabo2015,Sniekers2015}. Interestingly, Stein and Hung pointed out further work on calibration in the computer experiments literature, where it was demonstrated that the data collection mechanism was of importance for good calibration \citep{Zhu2006}. Note that a common assumption in these is that of noisy observations, which is not common in numerical analysis where functions are usually evaluated without noise (although the noisy case may still be relevant for certain applications). Nevertheless, these papers give us hope that similar calibration results could also be derived for Bayesian numerical methods.

\section*{Bayesian Priors, Posteriors and the Difficulty of Working with Measures}

The main premise required to justify the use of Bayesian methods is that we are able to come up with a good prior model which represents all our information prior to observing any data. This is clearly somewhat unreasonable, and it is debatable to what extent it is possible to construct such priors for the particular task of solving numerical problems. This issue is highlighted by Owen's claim that credible intervals obtained from Bayesian quadrature are `not anybody's belief'. Similar criticisms are of course also valid for most Bayesian numerical methods, and it is clear that much research is needed to help in the construction of such priors.

In this paper, we studied the case of Gaussian process priors, which is by far the most common choice for Bayesian numerical methods. The main reason for this choice is that working with general measures defined on infinite-dimensional domains is a notably challenging problem and requires the use of advanced sampling methods \citep{Cotter2013}. On the other hand, working with Gaussian processes is much simpler since linear operations preserve Gaussianity, and so boil down to linear algebra. As pointed out by Owen, Hickernell and Jagadeeswaran, this may still usually require $O(n^3)$ operations, but many approximation schemes for Gaussian processes exist. Some cheaper exact schemes based on particular choices of point sets or covariance functions have also recently been developed for the specific case of Bayesian numerical methods \cite{Karvonen2017symmetric,Jagadeeswaran2018}. Furthermore, we believe that the approach based on Fourier transforms suggested by Stein and Hung could also be fruitful.

Gaussian priors also lead to Bayesian numerical methods which are closely related to traditional methods. Indeed, connections between the posterior variance of Bayesian methods using Gaussian measures and worst-case results in reproducing kernel Hilbert spaces (RKHS) are well known and studied in the information-based complexity literature; see \cite{Ritter2000}. One issue which often leads to much confusion when making use of this connection, and which is highlighted by Stein, Hung, Hickernell and Jagadeeswaran, is the fact that draws from the Gaussian process are not in the space reproduced by the covariance function of the Gaussian process. The connection between Gaussian processes and RKHSs is nonetheless deep and offers opportunities to transfer results between the two frameworks, as highlighted by \cite{Kanagawa2018}. This may be extremely relevant to problems in numerical analysis as demonstrated by the Brownian bridge example of Stein and Hung, and by the recent work of \citep{Kanagawa2016} on misspecified Gaussian process priors.

For this reason, a central question is the following: How should we pick a prior for such an abstract quantity as an integrand? Clearly certain choices can be made based on prior knowledge of properties of the function, such as smoothness or periodicity, but we often still need to set certain amplitude or lengthscale parameters. These can be set using objective priors, as demonstrated in the paper and extended by Hickernell and Jagadeeswaran, or inferred using empirical Bayes approaches. These choices are clearly not easy to make and can have significant influence on the resulting algorithm.

The other option is of course to consider more complex priors. We agree with Stein and Hung that non-stationary priors could be a way forward, and are currently working in this direction. Another alternative to construct expressive priors would be the recent work on deep Gaussian process priors; see \cite{Dunlop2017}. Finally, we could move away from Gaussianity, allowing us more flexibility in the choice of functions on which to assign mass \citep{Cockayne2017BPNM}. Unfortunately, many of these alternatives come at the cost of additional computational requirements which may be prohibitive or require resorting to approximations which may not faithfully quantify uncertainty.

\section*{So, Should we be Bayesians for Numerical Problems?}

Regardless of the choice of prior, what is clear is that working in the space of measures is significantly harder than working in many of the function spaces commonly used in numerical analysis. As a consequence, the uncertainty quantification provided by Bayesian methods will usually come at a cost, and one should weigh this cost with the potential benefits before making a decision on whether to use a Bayesian numerical method.

A first application area where this may be useful is the field of inverse problems or computer experiments. The models studied here are often computationally expensive and based on systems of differential equations. The most popular algorithms for inference rely on repeated solution of numerical tasks such as interpolation, integration, linear algebra and optimisation. As pointed out by Stein and Hung, many of the issues facing Bayesian numerical methods, such as the choice of a prior, are also common in this literature.

Machine learning also offers many opportunities for the deployment of Bayesian numerical analysis since it is highly reliant on numerical methods and there has been a recent focus on better explaining the workings and failure modes of algorithms (that is, to provide interpretability and quantification of uncertainty). The difficulty of selecting priors can also be somewhat alleviated since we have an excellent source of prior information: the source code of the problem.

On the other hand, there is clearly no point providing uncertainty quantification for applications where uncertainty can be made negligible without excessive amounts of computation. For integration, this will usually be the case when the integrand and sampling mechanism are cheap. In this case, we would argue that being Bayesian about the numerical problem may not be necessary, or desirable.

\section*{Some Promising Research Directions}

With these application areas in mind, we conclude by highlighting areas which we believe could benefit the Bayesian approach to numerical analysis.

The first is the ability to jointly model the solution of several quantities of interest, for which data can be obtained simultaneously or sequentially. This could for example be the solutions of several related integrals, as was considered in \cite{Xi2018MultiOutput}, but could also include other problems such as solving differential equations, linear systems, optimisation  or a mixture of these problems. In such cases, data from one of the problems could be leveraged to improve the accuracy of the algorithms approximating the other problems. Such information sharing is extremely natural within a Bayesian framework and currently underexplored in the numerical analysis literature.

Another advantage which comes out of the Bayesian viewpoint is that of chains of computations and the ability to condition on all sources of uncertainty. Take for example some of the complex Bayesian inverse problems present in the computer experiments literature. There, it is common to have to sequentially solve numerical problems in both the forward and inverse problem. Here, reducing the discretisation error in each of these problems to a negligible level may not always be possible due to computational cost. Being Bayesian about numerical problems allow us to condition jointly on the uncertainty coming from the inverse problem and the uncertainty inherent due to numerical error using Bayes' theorem. We can hence obtain posteriors which are reflect our combined uncertainty due to all sources of discretisation \cite{Cockayne2017BPNM,Oates2017hydrocyclones}. The numerical analysis of large scientific codes is not straightforward in general, but Bayesian numerical analysis provides a coherent overarching framework.

\subsubsection*{Acknowledgments}

FXB was supported by the EPSRC grants [EP/L016710/1, EP/R018413/1]. MG was supported by the EPSRC grants [EP/J016934/3, EP/K034154/1, EP/P020720/1, EP/R018413/1], an EPSRC Established Career Fellowship, the EU grant [EU/259348]. FXB, CJO, MG \& DS were all supported by The Alan Turing Institute under grant [EP/N510129/1] and/or the Lloyds Register Foundation Programme on Data-Centric Engineering.

\bibliographystyle{plainnat}