\section{Introduction}\label{sec.1}
There is a growing demand in various fields for flexible models that are able to accurately characterize complex relations among a vector of binary responses $ {\bf y} =(y_1, \ldots, y_n)^{\intercal}$ and a set of predictors $ {\bf X} =( {\bf x} _{1}, \ldots,  {\bf x} _{n})^{\intercal}$, with $y_i \in \{0; 1\}$ and $ {\bf x} _{i}=(x_{i1}, \ldots, x_{iq})^{\intercal} \in \Bbb{R}^q$, for  each $i=1, \ldots, n$. Common solutions address this goal by replacing the  linear predictor $ {\bf X}  {\boldsymbol \beta} =( {\bf x} _1^{\intercal} {\boldsymbol \beta} , \ldots,  {\bf x} _n^{\intercal} {\boldsymbol \beta} )^{\intercal} \in  \Bbb{R}^n$ in the generalized linear model for $ {\bf y} $ \citep{nelder1972}, with a more flexible vector $ {\bf f} ( {\bf X} )=[f( {\bf x} _1), \ldots,f( {\bf x} _n)]^{\intercal} \in  \Bbb{R}^n$ that accounts for complex non--linear relations between the response and the predictors, thus enhancing predictive power. Notable examples of this approach within the Bayesian setting define $ {\bf f} ( {\bf X} )$ via additive trees \citep{chipman_2010}, Bayesian \textsc{p}--splines \citep{brezger2006} and Gaussian processes (\textsc{gp}) \citep{ras_2006}, among others. Motivated by the success of \textsc{gp}s for classification  \citep[e.g.,][]{neal1999, opper2000, de2005, chu2005, girolami2006,ras_2006,choudhuri2007,riihimaki2013}, we focus  on deriving improved methods to evaluate the predictive probabilities for the latter class of models under the probit link. In particular, we assume that the responses $y_i$, $i=1, \ldots, n$ are conditionally independent realizations from Bernoulli variables with probability parameters $\Phi[f( {\bf x} _i)]$, $i=1, \ldots, n$, where $\Phi[f( {\bf x} )]$ is the cumulative distribution function of a standard Gaussian  evaluated at $f( {\bf x} )$, whereas $f( {\bf x} )$   is assigned a \textsc{gp} prior with mean function $m( {\bf x} )=\mathbb{E}[f( {\bf x} )]$ and covariance kernel $K( {\bf x} , {\bf x} ')=\mathbb{E}\{[f( {\bf x} )-m( {\bf x} )] [f( {\bf x} ')-m( {\bf x} ')]\}$. Exploiting standard \textsc{gp}'s properties \citep[][]{ras_2006} and assuming no overlap in $ {\bf x} _1, \ldots,  {\bf x} _n$, such assumptions lead to the model
\begin{eqnarray}
p[ {\bf y}  \mid  {\bf f} ( {\bf X} )]= \prod\nolimits_{i=1}^n \Phi[f( {\bf x} _i)]^{y_i}\{1-\Phi[f( {\bf x} _i)] \}^{1-y_i}, \qquad \mbox{with} \ \ p[ {\bf f} ( {\bf X} )]= \phi_n[ {\bf f} ( {\bf X} )- {\boldsymbol \xi} ;  {\boldsymbol \Omega} ],
\label{eq1}
\end{eqnarray}
where $\phi_n[ {\bf f} ( {\bf X} )- {\boldsymbol \xi} ;  {\boldsymbol \Omega} ]$ denotes the density of a multivariate Gaussian distribution $\mbox{N}_n( {\boldsymbol \xi} ,  {\boldsymbol \Omega} )$ for $ {\bf f} ( {\bf X} )$, with mean vector $ {\boldsymbol \xi} =[m( {\bf x} _1), \ldots, m( {\bf x} _n)]^{\intercal}$, and $n \times n$ covariance matrix $ {\boldsymbol \Omega} $ having entries $ {\boldsymbol \Omega} _{[ii']}=K( {\bf x} _i, {\bf x} _{i'})$, for every $i=1, \ldots, n$ and $i'=1, \ldots, n$. Model \eqref{eq1} has attracted considerable interest due to its flexibility and its direct connection with binary discrete choice models based on Gaussian latent utilities  $z_i=f( {\bf x} _i)+ \varepsilon_i$, with $ \varepsilon_i \sim \mbox{N}(0,1)$, independently for $i=1, \ldots, n$  \citep{Albert_1993}. In fact, $p[y_i=1 \mid f( {\bf x} _i)]=\Phi[f( {\bf x} _i)]=p[z_i>0 \mid f( {\bf x} _i)]$. In such settings, a main goal of inference is to evaluate the predictive probabilities of new responses $y_{n+1}$, defined as
\begin{eqnarray}
\begingroup\textstyle
p(y_{n+1}= 1 \mid   {\bf y} )=1-p(y_{n+1}= 0\mid  {\bf y} )={\int} \Phi[f( {\bf x} _{n+1})] \cdot \left[ {\int} p[f( {\bf x} _{n+1}),  {\bf f} ( {\bf X} ) \mid   {\bf y} ] \mbox{d}  {\bf f} ( {\bf X} ) \right]\mbox{d} f( {\bf x} _{n+1}),
\endgroup
\label{eq2}
\end{eqnarray}
where $p[f( {\bf x} _{n+1}),  {\bf f} ( {\bf X} ) \mid   {\bf y} ]$ is the joint posterior density of $[f( {\bf x} _{n+1}), {\bf f} ( {\bf X} ) ]$ under model \eqref{eq1}, which seems not available in closed form due to the apparent absence  of conjugacy between the probit likelihood and the multivariate Gaussian prior for $[f( {\bf x} _{n+1}), {\bf f} ( {\bf X} ) ]$ under model \eqref{eq1}. This issue has motivated extensive research to compute predictive probabilities in probit models with multivariate Gaussian priors either via  Monte Carlo methods relying on samples from $p[f( {\bf x} _{n+1}),  {\bf f} ( {\bf X} ) 	\mid   {\bf y} ]$ \citep{neal1999,Albert_1993,de2005, holmes_2006,choudhuri2007,pakman_2014, hoff_2014,Durante2018} or by deriving tractable approximations of $p[f( {\bf x} _{n+1}),  {\bf f} ( {\bf X} ) 	\mid   {\bf y} ]$ \citep{consonni_2007,girolami2006, chu2005,ras_2006,riihimaki2013} that allow simple evaluation of \eqref{eq2}. These methods provide state--of--the--art solutions in small--to--moderate dimensional settings, but tend to become rapidly inaccurate or computationally impractical in higher dimension \citep{chopin_2017,Johndrow2018,Durante2018, fasano2019asymptotically}. This issue is inherently found in probit \textsc{gp}s where, by definition, there are $p \approx n$ parameters in $ {\bf f} ( {\bf X} )$, with the sample size $n$ being typically large in most studies.

In this article we aim to cover the above gap by providing novel closed--form expressions for the predictive probabilities in probit  \textsc{gp}s along with improved methods to evaluate the involved quantities in high dimensions. More specifically, in Section \ref{sec.2.1} we first derive a closed--form expression for the marginal likelihood $p( {\bf y} )$ under model \eqref{eq1}, and then exploit this result to  show that $p(y_{n+1}= 1\mid  {\bf y} )$ can be expressed as the ratio between cumulative distribution functions of multivariate Gaussians with dimensions $n+1$ and $n$, respectively. To overcome the known issues associated with the evaluation of these two quantities in high dimensions \citep[][]{chop_2011,botev_2017,cao2019,cao2020} we introduce an error--reduction technique for computing ratios of Gaussian cumulative distribution functions that builds on the tile--low--rank   method in \citet{cao2020} and substantially reduces the computational time of state--of--the--art strategies such as minimax tilting methods \citep{botev_2017} and Hamiltonian Monte Carlo samplers (\textsc{stan}) \citep{hoff_2014}, without affecting accuracy. To further improve the scalability of the methods presented in Section \ref{sec.2.1}, we show in Section \ref{sec.2.2} that  $p(y_{n+1}= 1\mid  {\bf y} )$ has also an alternative representation based on functionals of multivariate truncated normals, and we address the intractability of such variables in high dimensions by proposing a  variational approximation based on univariate truncated normals which allows accurate and efficient evaluation of predictive probabilities in high dimensions, thus scaling to settings where available strategies are unfeasible. Such results are also related to the conditional distribution of the \textsc{gp} given the binary responses which we show to coincide with a unified skew--normal  (\textsc{sun}) \citep{arellano_2006} by adapting recent results in \citet{Durante2018}  on classical Bayesian probit regression. The magnitude of the improvements provided by the new methods presented in Sections \ref{sec.2.1}--\ref{sec.2.2} relative to state--of--the--art competitors is illustrated in simulations in Section \ref{sec.3} and in an environmental application to Saudi Arabia windspeed in Section \ref{sec.4}. Section \ref{sec.5} contains concluding remarks, whereas all the proofs can be found in the Appendix.

\section{Computation of Predictive Probabilities in Probit Gaussian Processes}
\label{sec.2}
In Sections \ref{sec.2.1} and \ref{sec.2.2} we present novel expressions for the predictive probabilities in probit Gaussian processes along with improved methods to evaluate efficiently the involved quantities in high dimensions.

\subsection{Evaluation via multivariate Gaussian probability ratios}
\label{sec.2.1}
To introduce the closed--form expression for $p(y_{n+1}= 1\mid  {\bf y} )$ based on ratios of multivariate Gaussian cumulative distribution functions, first notice that by leveraging known properties of Gaussian variables, the probit likelihood in \eqref{eq1} can be re--expressed as $p[ {\bf y}  \mid  {\bf f} ( {\bf X} )]= \prod\nolimits_{i=1}^n \Phi[f( {\bf x} _i)]^{y_i}\{1-\Phi[f( {\bf x} _i)] \}^{1-y_i}= \prod_{i=1}^n \Phi[(2y_i-1)f( {\bf x} _i)]=\Phi_n[ {\bf D}  {\bf f} ( {\bf X} );  {\bf I} _n]$, where $ {\bf D} =\mbox{diag}[(2y_1-1), \ldots, (2y_n-1)]$, and $\Phi_n[ {\bf D}  {\bf f} ( {\bf X} );  {\bf I} _n]$ is the cumulative distribution function of a zero--mean $n$--variate Gaussian with identity covariance matrix $ {\bf I} _n$, evaluated at $ {\bf D}  {\bf f} ( {\bf X} )$. Leveraging this  form and adapting Lemma 7.1 in \citet{azzalini_2013} to our setting, we can easily express the marginal likelihood under model \eqref{eq1} as
\begin{eqnarray}
\begingroup\textstyle
p( {\bf y} )={\int} \Phi_n[ {\bf D}  {\bf f} ( {\bf X} );  {\bf I} _n] \phi_n[ {\bf f} ( {\bf X} )- {\boldsymbol \xi} ;  {\boldsymbol \Omega} ]\mbox{d}  {\bf f} ( {\bf X} ) = \Phi_n( {\bf D}  {\boldsymbol \xi} ;  {\bf I} _n+ {\bf D}  {\boldsymbol \Omega}  {\bf D} ^{\intercal}).
\endgroup
\label{eq3}
\end{eqnarray}
Equation \eqref{eq3} provides a closed--form expression that can be useful to estimate fixed parameters in the \textsc{gp} kernel via maximization of $p( {\bf y} )$ and, as shown in Proposition \ref{prop1},  is also a key to obtain closed--form expressions for $p(y_{n+1}= 1|  {\bf y} )$.

\begin{Proposition}
Under \eqref{eq1}, the predictive probability for a new response $y_{n+1} \in \{0;1\}$ with predictor $ {\bf x} _{n+1} \in \Bbb{R}^q$ is
\begin{eqnarray}
p(y_{n+1}= 1 \mid  {\bf y} )=1-p(y_{n+1}= 0 \mid  {\bf y} )=\frac{\Phi_{n+1}( {\bf D} ^* {\boldsymbol \xi} ^*;  {\bf I} _{n+1}+ {\bf D} ^* {\boldsymbol \Omega} ^* {\bf D} ^{*\intercal})}{\Phi_n( {\bf D}  {\boldsymbol \xi} ;  {\bf I} _n+ {\bf D}  {\boldsymbol \Omega}  {\bf D} ^{\intercal})},
\label{eq4}
\end{eqnarray}
where $ {\bf D} ^*=\mbox{\normalfont diag}[(2y_1-1), \ldots, (2y_n-1), 1]$, $ {\boldsymbol \xi} ^*=[ {\boldsymbol \xi} ^{\intercal},m( {\bf x} _{n+1})]^{\intercal}$, and $ {\boldsymbol \Omega} ^*$ is obtained by adding an additional row and column to $ {\boldsymbol \Omega} $ defined as $ {\boldsymbol \Omega} _{[n+1, \cdot]}^{*\intercal}= {\boldsymbol \Omega} ^*_{[\cdot, n+1]}=[K( {\bf x} _{n+1}, {\bf x} _1), \ldots, K( {\bf x} _{n+1}, {\bf x} _n), K( {\bf x} _{n+1}, {\bf x} _{n+1})]^{\intercal}$.
\label{prop1}
\end{Proposition}

To prove Proposition \ref{prop1}, it is sufficient to notice that, by the Bayes rule, $p(y_{n+1}= 1 \mid  {\bf y} )=p(y_{n+1}= 1,  {\bf y} )/p( {\bf y} )$ where $p(y_{n+1}= 1,  {\bf y} )$ and $p( {\bf y} )$ are the marginal likelihoods of $(y_{n+1}= 1,  {\bf y} )$ and $ {\bf y} $, respectively, under model \eqref{eq1}. Replacing such quantities with their closed--form expression in \eqref{eq3}, leads to \eqref{eq4}. See the Appendix for a more detailed proof which also includes additional clarifications on equation \eqref{eq3}.

Evaluation of \eqref{eq4} requires the calculation of cumulative distribution functions of multivariate Gaussians, which is known to be a challenging task in high dimensions \citep[e.g.,][]{genz_1992,chop_2011,botev_2017,genton2018,cao2019,cao2020}. Recent advances via minimax tilting methods \citep{botev_2017} allow accurate evaluation of such quantities, but face an increased computational cost which makes such strategies rapidly impractical as $n$ grows. A more scalable solution can be found in the separation--of--variable (\textsc{sov}) algorithm originally introduced by \citet{genz_1992} and subsequently improved in terms of scalability by \citet{cao2020}. Such a routine decomposes the generic multivariate Gaussian probability $\Phi_n( {\bf a} , {\bf b} ; {\boldsymbol \Sigma} )$ as
\begin{equation}
\begingroup\textstyle
\Phi_n( {\bf a} , {\bf b} ; {\boldsymbol \Sigma} )=\int^{ {\bf b} }_{ {\bf a} } \phi_n( {\bf u} ;  {\boldsymbol \Sigma} )\mbox{d} {\bf u} =(e_1-d_1) \int^1_{0}(e_2-d_2) \cdots  \int^1_{0}(e_n-d_n)  \int^1_{0}\mbox{d}  {\bf w} =\mathbb{E}_{ {\bf w} }[(e_1-d_1) \cdots (e_n-d_n)],
\label{eq5}
\endgroup
\end{equation}
with $d_i = \Phi(\{a_i - \sum_{j=1}^{i-1}l_{ij} \Phi^{-1}[d_j + w_j(e_j - d_j)]\}/l_{ii})$ and $e_i = \Phi(\{b_i - \sum_{j=1}^{i-1}l_{ij} \Phi^{-1}[d_j + w_j(e_j - d_j)]\}/l_{ii})$ for each $i=1, \ldots, n$, where  $l_{ij}$ is the $(ij)$-th coefficient in the lower Cholesky factor of $ {\boldsymbol \Sigma} $, and $ {\bf w} =(w_1, \ldots, w_{n-1})^{\intercal}$ denotes a vector with uniform entries $w_j  \sim \mbox{U}(0,1)$, for $j=1, \ldots, n-1$. Such a decomposition transforms the integration region into the unit hypercube, thus allowing the evaluation of  $\Phi_n( {\bf a} , {\bf b} ; {\boldsymbol \Sigma} )$ via functionals of uniform densities.   To further improve the quality of the above estimator, more recent implementations \citep{trinh2015} combine \eqref{eq5} with a univariate reordering  preconditioner that rearranges the integration variables and produces the corresponding Cholesky factor simultaneously at the same $O(n^3)$ cost of the Cholesky factorization. This ordering strategy processes the integration variables from left to right iteratively and, at each step, it switches the original integration variable with the one having the narrowest conditional integration limits on its right side. Positioning left the integration variables with narrower integration limits is shown in \cite{trinh2015} and \cite{cao2020} to improve the Monte Carlo convergence rate of \eqref{eq5}, whose integrand is evaluated $R$ times --- corresponding to the Monte Carlo sample size --- each of which has a cost of $O(n^2)$. Such costs allow the implementation of this strategy in settings with $n \leq 1{,}000$, thus motivating more scalable options in high dimensions. \citet{cao2020} address this issue via a tile--low--rank representation for $ {\boldsymbol \Sigma} $ that reduces the cost of the \textsc{sov} algorithm by substituting the dense matrix--vector multiplication with the low--rank matrix--vector multiplication. A compatible block--reordering is also introduced in place of the univariate reordering to improve the convergence rate at the same cost as the low--rank Cholesky factorization. Specifically, the block--reordering orders integration variables on the block level based on crude estimates of the block--wise marginal probabilities. Both the block--reordering and the tile--low--rank  version of the \textsc{sov} algorithm reach their optimal complexities of $O(n^{5/2})$ and $O(n^{3/2})$, respectively, when the block size in the tile--low--rank representation is $n^{1/2}$, thus reducing the computational complexity of the classical  \textsc{sov} algorithm by $n^{1/2}$ and allowing implementation in  tens of thousands of dimensions.

Although the above techniques can be effectively implemented to evaluate multivariate Gaussian probabilities as in \eqref{eq3}, the calculation of ratios among such quantities as in \eqref{eq4} has typically high accuracy requirements. Unfortunately, as discussed in \citet{botev_2017} and \citet{cao2020}, the estimation errors of tail multivariate Gaussian probabilities, that also include the cumulative distribution function, can be as large as the probability estimates themselves when $n$ is in hundreds to thousands of dimensions, thus producing unreliable ratio estimates. To address this issue, we propose an error--reduction technique that avoids computing  the numerator and the denominator in \eqref{eq4}  separately, but combines their evaluation under the tile--low--rank representation. Indeed, as is clear from Proposition \ref{prop1}, the denominator in \eqref{eq4} is the marginalization of the numerator over the last integration variable. Therefore, keeping the general notation of the \textsc{sov} algorithm and leveraging equation \eqref{eq5}, expression \eqref{eq4} can be re--written in the general form
\begin{eqnarray}
\frac{\Phi_{n+1}( {\bf a} , {\bf b} ; {\boldsymbol \Sigma} )}{\Phi_{n}( {\bf a} _{-(n+1)}, {\bf b} _{-(n+1)}; {\boldsymbol \Sigma} _{-(n+1)})}=\frac{\mathbb{E}_{ {\bf w} }[(e_1-d_1)\cdots (e_{n}-d_{n})\cdot(e_{n+1}-d_{n+1})]}{\mathbb{E}_{ {\bf w} _{-n}}[(e_1-d_1)\cdots (e_{n}-d_{n})]},
\label{eq6}
\end{eqnarray}
where $e_i$ and $d_i$ are defined as in equation \eqref{eq5} for $i=1, \ldots, n+1$, whereas $ {\bf a} _{-(n+1)}$, $ {\bf b} _{-(n+1)}$ and $ {\bf w} _{-n}$ are obtained by removing the $(n+1)$-th entry in $ {\bf a} $ and $ {\bf b} $, and the $n$-th entry in $ {\bf w} $, respectively. Similarly, $ {\boldsymbol \Sigma} _{-(n+1)}$ coincides with $ {\boldsymbol \Sigma} $ without the $(n+1)$-th row and column. As is clear from \eqref{eq6}, the quantities $(e_1-d_1), \ldots, (e_{n}-d_{n})$ are the  same deterministic functions of $ {\bf w} $ both in the numerator and in the denominator, and hence, using the same set of Monte Carlo samples $ {\bf w} $ in the $n$--dimensional hypercube for estimating the two expectations might significantly reduce the estimation error of their ratio. In particular, our proposed ratio estimator is
\begin{eqnarray}
\hat{p}(y_{n+1}= 1 \mid  {\bf y} )=\frac{\sum_{r=1}^R[e_1( {\bf w} ^{(r)}_{-n})-d_1( {\bf w} ^{(r)}_{-n})] \cdots [e_n( {\bf w} ^{(r)}_{-n})-d_n( {\bf w} ^{(r)}_{-n})]\cdot [e_{n+1}( {\bf w} ^{(r)})-d_{n+1}( {\bf w} ^{(r)})]}{\sum_{r=1}^R[e_1( {\bf w} ^{(r)}_{-n})-d_1( {\bf w} ^{(r)}_{-n})] \cdots [e_n( {\bf w} ^{(r)}_{-n})-d_n( {\bf w} ^{(r)}_{-n})]},
\label{eq7}
\end{eqnarray}
where the generic $e_i( {\bf w} ^{(r)})$ and $d_i( {\bf w} ^{(r)})$ denote the values of $e_i$ and $d_i$ in \eqref{eq5} evaluated at the Monte Carlo  sample $ {\bf w} ^{(r)}$ of $ {\bf w} $. This estimator is asymptotically unbiased because the numerator and the denominator converge to $\mathbb{E}_{ {\bf w} }[(e_1-d_1)\cdots (e_{n}-d_{n}) \cdot (e_{n+1}-d_{n+1})]$ and $\mathbb{E}_{ {\bf w} _{-n}}[(e_1-d_1)\cdots (e_{n}-d_{n})]$, respectively, and hence equation \eqref{eq7} converges to \eqref{eq6} in probability. Moreover, equation~\eqref{eq7} is guaranteed to be in $(0, 1)$, thus producing an estimator whose variance is always smaller than $0.25$. This is not the case when  the numerator and the denominator in \eqref{eq4} are estimated separately. Indeed, as discussed in \citet{botev_2017} and  \citet{cao2020}, when $n$ is high the estimation errors of the two cumulative distribution functions in \eqref{eq4} are often as large as the estimates themselves, thus producing estimated ratios possibly outside of the range $(0, 1)$ and with high variance.

The pseudo--code for evaluating~\eqref{eq4} via the estimator outlined in~\eqref{eq7} is provided in Algorithm~\ref{algo:1}. In step~{\bf[b]}, the block reordering produces a new variable order that is used to reorder the integration limits $ {\bf a} $ and $ {\bf b} $, whereas in step~{\bf[c]} the inverse matrices of the diagonal blocks of $ {\bf L} $ computed in step~{\bf[b]} are recycled to maximize efficiency. Also the quantities in step~{\bf[d.1]} do not need to be re--evaluated every time a new prediction is required since they only depend on the observed training data, and hence such quantities can be pre--computed and stored separately.

\begin{algorithm*}[t]
	\caption{Compute \eqref{eq4} via the estimator \eqref{eq7}}
	\vspace{5pt}
{\bf [a]} Define $ {\bf a} =-\infty$, $ {\bf b} = {\bf D} ^* {\boldsymbol \xi} ^*$, $ {\boldsymbol \Sigma} =  {\bf I} _{n+1}+ {\bf D} ^* {\boldsymbol \Omega} ^* {\bf D} ^{*\intercal}$, $ {\bf a} _{-(n+1)}=-\infty$, $ {\bf b} _{-(n+1)}= {\bf D}  {\boldsymbol \xi} $, $ {\boldsymbol \Sigma} _{-(n+1)}=  {\bf I} _{n}+ {\bf D}  {\boldsymbol \Omega}  {\bf D} ^{\intercal}$ and let $ {\bf w} ^{(r)}, \ldots,  {\bf w} ^{(R)}$ denote uniform samples from the $n$--dimensional unit hypercurbe.\\
\vspace{5pt}
		{\bf [b]} Apply {\bf block reordering} \citep{cao2020} to ($ {\bf a} _{-(n+1)}$, $ {\bf b} _{-(n+1)}$, $ {\boldsymbol \Sigma} _{-(n+1)}$), which produces the tile--low--rank Cholesky factor $ {\bf L} $ after variable reordering, and reorder $ {\bf a} _{-(n+1)}$ and $ {\bf b} _{-(n+1)}$ accordingly.\\
\vspace{5pt}
		{\bf [c]} Compute $ {\bf L} _{[n+1,1:n+1]}$ using $ {\boldsymbol \Sigma} $ and $ {\bf L} $.\\
\vspace{5pt}
{\bf [d]} Obtain the quantities required to evaluate equation \eqref{eq7}. \\
		\For{r=1,\ldots,R}{
		{\bf [d.1]} Compute the differences $e_1( {\bf w} ^{(r)}_{-n})-d_1( {\bf w} ^{(r)}_{-n}), \ldots, e_n( {\bf w} ^{(r)}_{-n})-d_n( {\bf w} ^{(r)}_{-n})$ by applying the tile--low--rank variant of equation~\eqref{eq5} \citep{cao2020} to $( {\bf a} ,  {\bf b} ,  {\bf L} )$. Store also the vector $ {\bf v} ^{(r)}=[\Phi^{-1}\{d_1( {\bf w} ^{(r)}_{-n})+w^{(r)}_1[e_1( {\bf w} ^{(r)}_{-n})-d_1( {\bf w} ^{(r)}_{-n})]\}, \ldots, \Phi^{-1}\{d_n( {\bf w} ^{(r)}_{-n})+w^{(r)}_n[e_n( {\bf w} ^{(r)}_{-n})-d_n( {\bf w} ^{(r)}_{-n})]\}]^{\intercal}$. \\
	\vspace{3pt}
	{\bf [d.2]} Set $e_{n+1}( {\bf w} ^{(r)})-d_{n+1}( {\bf w} ^{(r)})=\Phi\left[\frac{b_{n+1}{-}\  {\bf L} _{[n+1,1:n]} {\bf v} ^{(r)}}{l_{n+1,n+1}}\right]-\Phi\left[\frac{a_{n+1}{-}\  {\bf L} _{[n+1,1:n]} {\bf v} ^{(r)}}{l_{n+1,n+1}}\right]$.
	}
		\vspace{3pt}
	{\bf [e]}  Estimate~\eqref{eq4} via Monte Carlo as in~\eqref{eq7} using the quantities computed in step {\bf [d]}.
			\vspace{2pt}
		\label{algo:1}
\end{algorithm*}

\subsection{Evaluation via functionals of multivariate truncated normals}
\label{sec.2.2}
\vspace{3pt}
The methods presented in Section \ref{sec.2.1} allow substantial improvements in terms of accuracy and scalability in the evaluation of the predictive probabilities under probit \textsc{gp}s. However, the cost of the tile--low--rank version of the \textsc{sov} algorithm might still be too expensive in high dimensional settings with a large $n$. To further reduce computational times, we derive an alternative expression for $p(y_{n+1}= 1 |  {\bf y} )$ relying on functionals of multivariate truncated normals which are then approximated via mean--field variational Bayes \citep[e.g.,][]{blei2017} to facilitate simple evaluation of   $p(y_{n+1}= 1 \mid  {\bf y} )$ using Monte Carlo samples from univariate truncated normals.

To derive this alternative expression, we shall first notice that the joint posterior $p[f( {\bf x} _{n+1}),  {\bf f} ( {\bf X} ) \mid   {\bf y} ]$ in \eqref{eq2} can be factorized as $p[f( {\bf x} _{n+1}) \mid  {\bf f} ( {\bf X} )]\cdot p[ {\bf f} ( {\bf X} ) \mid  {\bf y} ]$, provided that $f( {\bf x} _{n+1})$ does not appear in the likelihood for $ {\bf y} $, which is true because there is no overlap among predictors. Exploiting the well--known properties of \textsc{gp}s \citep[][]{ras_2006}, the first factor $p[f( {\bf x} _{n+1}) \mid  {\bf f} ( {\bf X} )]$ in the above expression can be easily derived by applying the closure under conditioning property of multivariate Gaussians, thus obtaining the univariate normal density
\begin{eqnarray}
\begin{split}
& p[f( {\bf x} _{n+1}) \mid  {\bf f} ( {\bf X} )] \\
&=\phi(f( {\bf x} _{n+1}) -[m( {\bf x} _{n+1})+ {\boldsymbol \Omega} ^*_{[n+1,1:n]} {\boldsymbol \Omega} ^{-1}( {\bf f} ( {\bf X} )- {\boldsymbol \xi} )]; K( {\bf x} _{n+1}, {\bf x} _{n+1})- {\boldsymbol \Omega} ^*_{[n+1,1:n]} {\boldsymbol \Omega} ^{-1} {\boldsymbol \Omega} ^*_{[1:n,n+1]} ),\\
&=\phi(f( {\bf x} _{n+1}) -[\mu_{x_{n+1}}+ {\bf H} _{x_{n+1}} {\bf f} ( {\bf X} )]; \sigma_{x_{n+1}}^2 ),
\end{split}
\label{eq8}
\end{eqnarray}
where $ {\bf H} _{x_{n+1}}{=}\  {\boldsymbol \Omega} ^*_{[n+1,1:n]} {\boldsymbol \Omega} ^{-1}$, $\mu_{x_{n+1}}{=}\ m( {\bf x} _{n+1})- {\bf H} _{x_{n+1}} {\boldsymbol \xi} $ and $ \sigma_{x_{n+1}}^2{=}\ K( {\bf x} _{n+1}, {\bf x} _{n+1})- {\boldsymbol \Omega} ^*_{[n+1,1:n]} {\boldsymbol \Omega} ^{-1} {\boldsymbol \Omega} ^*_{[1:n,n+1]} $, whereas the other quantities are defined as in equation \eqref{eq1} and  \eqref{eq4}. Adapting the recent conjugacy results for probit models with Gaussian priors in  \citet{Durante2018} to the \textsc{gp} setting, it is also possible to express $p[ {\bf f} ( {\bf X} ) \mid  {\bf y} ]$ as the density of the unified skew--normal (\textsc{sun}) \citep{arellano_2006} $\mbox{SUN}_{n,n}[ {\boldsymbol \xi} , {\boldsymbol \Omega} ,\bar{ {\boldsymbol \Omega} } {\boldsymbol \omega}  {\bf D} ^{\intercal} {\bf s} ^{-1}, {\bf s} ^{-1} {\bf D}  {\boldsymbol \xi} ,  {\bf s} ^{-1}( {\bf D}  {\boldsymbol \Omega}  {\bf D} ^{\intercal}{+} {\bf I} _n) {\bf s} ^{-1}]$, with $ {\bf s} =[( {\bf D}  {\boldsymbol \Omega}  {\bf D} ^{\intercal}{+} {\bf I} _n)\odot  {\bf I} _n]^{1/2}$, $\bar{ {\boldsymbol \Omega} }= {\boldsymbol \omega} ^{-1} {\boldsymbol \Omega}  {\boldsymbol \omega} ^{-1}$ and $ {\boldsymbol \omega} =( {\boldsymbol \Omega}  \odot  {\bf I} _n)^{1/2}$. Indeed, recalling the results discussed in Sections \ref{sec.1}--\ref{sec.2.1} and applying the Bayes rule, we have that $p[ {\bf f} ( {\bf X} ) \mid  {\bf y} ] \propto p[ {\bf f} ( {\bf X} )] \cdot p[ {\bf y}  \mid  {\bf f} ( {\bf X} )]=\phi_n[ {\bf f} ( {\bf X} )- {\boldsymbol \xi} ;  {\boldsymbol \Omega} ]\cdot \Phi_n[ {\bf D}  {\bf f} ( {\bf X} );  {\bf I} _n]$ which coincides with the kernel of a \textsc{sun} density as discussed in the proof of Theorem 1 in  \citet{Durante2018}. Such a class of random variables introduces asymmetric shapes in Gaussian densities via a skewness--inducing mechanism driven by the cumulative distribution function of an $n$--variate Gaussian with a full-rank covariance matrix. Hence, the evaluation of $ p[ {\bf f} ( {\bf X} ) \mid  {\bf y} ]$ still requires calculation of multivariate Gaussian probabilities, leading to the same issues discussed in Section \ref{sec.2.1}; see \citet{arellano_2006}, \citet{azzalini_2013} and \citet{Durante2018} for an in--depth discussion on the properties of \textsc{sun} variables for posterior inference.

A possibility to address the above issue is to leverage the discrete-choice interpretation of the probit \textsc{gp} introduced in Section \ref{sec.1}. Under this alternative representation, model \eqref{eq1} can be equivalently re--expressed as $y_i=1(z_i>0)$, with $[z_i \mid  {\bf f} ( {\bf x} _i)] \sim \mbox{N}[ {\bf f} ( {\bf x} _i),1]$, independently for $i=1, \ldots, n$, and $ {\bf f} ( {\bf X} )=[f( {\bf x} _1), \ldots,f( {\bf x} _n)]^{\intercal} \sim \mbox{N}_n( {\boldsymbol \xi} ,  {\boldsymbol \Omega} )$. Adapting the results in \citet{holmes_2006} to our \textsc{gp} setting, the joint posterior $p[ {\bf f} ( {\bf X} ),  {\bf z}  \mid  {\bf y} ]$ of $ {\bf f} ( {\bf X} )$ and the augmented data $ {\bf z} =(z_1, \ldots, z_n)^{\intercal}$, factorizes as $p[ {\bf f} ( {\bf X} ) \mid   {\bf z} ]\cdot p( {\bf z}  \mid  {\bf y} )$, with
\begin{eqnarray}
\begin{split}
p[ {\bf f} ( {\bf X} ) \mid  {\bf z} ]& =\phi_n[ {\bf f} ( {\bf X} )-( {\boldsymbol \Omega} ^{-1}+ {\bf I} _n)^{-1}( {\boldsymbol \Omega} ^{-1} {\boldsymbol \xi} + {\bf z} ); ( {\boldsymbol \Omega} ^{-1}+ {\bf I} _n)^{-1}]=\phi_n( {\bf f} ( {\bf X} )-[ {\boldsymbol \mu} _{ {\bf X} }+ {\boldsymbol \Sigma} _{ {\bf X} } {\bf z} ];  {\boldsymbol \Sigma} _{ {\bf X} }), \\
p( {\bf z}  \mid  {\bf y} )& \propto \phi_n[ {\bf z}  -  {\boldsymbol \xi} ;  {\bf I} _n+ {\boldsymbol \Omega} ]\prod\nolimits_{i=1}^n1[(2y_i-1)z_i>0]=\phi_n[ {\bf z}  -  {\boldsymbol \xi} ;  {\boldsymbol \Sigma} _{ {\bf z} }]\prod\nolimits_{i=1}^n1[(2y_i-1)z_i>0],
\end{split}
\label{eq9}
\end{eqnarray}
where $ {\boldsymbol \Sigma} _{ {\bf X} }=( {\boldsymbol \Omega} ^{-1}+ {\bf I} _n)^{-1}$, $ {\boldsymbol \mu} _{ {\bf X} }= {\boldsymbol \Sigma} _{ {\bf X} } {\boldsymbol \Omega} ^{-1} {\boldsymbol \xi} $ and $ {\boldsymbol \Sigma} _{ {\bf z} }= {\bf I} _n+ {\boldsymbol \Omega} $. Therefore, the joint posterior  density $p[ {\bf f} ( {\bf X} ) \mid   {\bf z} ]\cdot p( {\bf z}  \mid  {\bf y} )$ factorizes as the product of a Gaussian for $p[ {\bf f} ( {\bf X} ) \mid  {\bf z} ]$ and a multivariate truncated normal for $p( {\bf z}  \mid  {\bf y} )$ obtained via component--wise truncation of $\mbox{N}( {\boldsymbol \xi} , {\boldsymbol \Sigma} _{ {\bf z} })$ above or below $0$, depending on whether $y_i=1$ or $y_i=0$, respectively, for each $i=1, \ldots, n$. As shown in Proposition \ref{prop2}, by combining equations \eqref{eq8}--\eqref{eq9} with Lemma 7.1 in \citet{azzalini_2013}, it is possible to obtain an alternative expression for $p(y_{n+1}= 1 |  {\bf y} )$ based on functionals of multivariate truncated normals. See the Appendix for a detailed proof.
\begin{Proposition}
Under \eqref{eq1}, the predictive probability for a new response $y_{n+1} \in \{0;1\}$ with predictor $ {\bf x} _{n+1} \in \Bbb{R}^q$ is
\begin{eqnarray}
\begin{split}
&p(y_{n+1}= 1 \mid  {\bf y} )=1-p(y_{n+1}=0\mid  {\bf y} )= \mathbb{E}_{ {\bf z}  \mid  {\bf y} }[\mathbb{E}_{ {\bf f} ( {\bf X} ) \mid  {\bf z} }(\mathbb{E}_{f( {\bf x} _{n+1}) \mid  {\bf f} ( {\bf X} )}\{\Phi[f( {\bf x} _{n+1})]\})],\\
&= \mathbb{E}_{ {\bf z}  \mid  {\bf y} }[\mathbb{E}_{ {\bf f} ( {\bf X} ) \mid  {\bf z} }[\Phi(\mu_{x_{n+1}}{+} {\bf H} _{x_{n+1}} {\bf f} ( {\bf X} );1{+}\sigma_{x_{n+1}}^2 )]],\\
&= \mathbb{E}_{ {\bf z}  \mid  {\bf y} }[\Phi(\mu_{x_{n+1}}{+}  {\bf H} _{x_{n+1}} [ {\boldsymbol \mu} _{ {\bf X} }{+}\  {\boldsymbol \Sigma} _{ {\bf X} } {\bf z} ];1{+}\sigma_{x_{n+1}}^2{+}  {\bf H} _{x_{n+1}}  {\boldsymbol \Sigma} _{ {\bf X} } {\bf H} ^{\intercal}_{x_{n+1}})], \ \ \
\label{eq10}
\end{split}
\end{eqnarray}
where the different quantities in \eqref{eq10} are defined as in equations \eqref{eq8} and \eqref{eq9}, whereas $\mathbb{E}_{ {\bf z}  \mid  {\bf y} }[\cdot]$ denotes the expectation with respect to the multivariate truncated normal density $p( {\bf z}  \mid  {\bf y} )$ in \eqref{eq9}.
\label{prop2}
\end{Proposition}

Leveraging Proposition \ref{prop2} it is possible to evaluate $p(y_{n+1}= 1 \mid  {\bf y} )$ via Monte Carlo methods based on independent samples from the multivariate truncated normal with density as in \eqref{eq9}, thus producing the estimate $\hat{p}(y_{n+1}= 1 \mid  {\bf y} )=1-\hat{p}(y_{n+1}= 0 \mid  {\bf y} )=\sum_{r=1}^R\Phi(\mu_{x_{n+1}}+  {\bf H} _{x_{n+1}} [ {\boldsymbol \mu} _{ {\bf X} }+  {\boldsymbol \Sigma} _{ {\bf X} } {\bf z} ^{(r)}];1+\sigma_{x_{n+1}}^2+  {\bf H} _{x_{n+1}}  {\boldsymbol \Sigma} _{ {\bf X} } {\bf H} ^{\intercal}_{x_{n+1}})/R$, where $ {\bf z} ^{(1)}, \ldots,  {\bf z} ^{(R)}$ denote independent and identically distributed samples from $p( {\bf z}  \mid  {\bf y} )$. Unfortunately, sampling from multivariate truncated normals in settings where $n$ is larger than a few hundreds raises the same computational issues discussed in Section \ref{sec.2.1}, i.e., the evaluation of multivariate Gaussian cumulative distribution functions \citep{holmes_2006,botev_2017, pakman_2014,Durante2018,fasano2019asymptotically}.

\begin{algorithm*}[t]
	\caption{Compute \eqref{eq10} via Monte Carlo based on a mean--field approximation of $p( {\bf z}  \mid  {\bf y} )$}
\vspace{3pt}
{\bf [a]}  Initialize $ {\bf z} ^{(0)}{=}\ [0, \ldots, 0]^{\intercal}$.	\\
\vspace{5pt}
{\bf [b]} Apply the \textsc{cavi} algorithm to obtain the optimal mean--field approximation $q^*( {\bf z} )=\prod_{i=1}^n q^*(z_i)$ for $p( {\bf z}  |  {\bf y} )$.	\\
\vspace{1pt}
		\For{t=1 until convergence}{
		\For{i=1, \ldots, n}{
		Set the univariate truncated normal approximating density $q^{(t)}(z_i)$ for $z_i$ at step $t$ equal to
		\begin{eqnarray*}
		q^{(t)}(z_i) \propto \phi\{z_i-[\xi_i+ {\bf H} _{z_i}( {\bf z} ^{(t-1)}_{-i}- {\boldsymbol \xi} _{-i})]; \sigma_{z_i}^{2}\}1[(2y_i-1)z_i>0]
		\end{eqnarray*}
		 with $ {\bf z} ^{(t-1)}_{-i}=[\mathbb{E}_{q^{(t)}(z_1)}(z_1), \ldots, \mathbb{E}_{q^{(t)}(z_{i-1})}(z_{i-1}),\mathbb{E}_{q^{(t-1)}(z_{i+1})}(z_{i+1}), \ldots, \mathbb{E}_{q^{(t-1)}(z_{n})}(z_{n})]^{\intercal}$.
		}
		\vspace{3pt}
		{\bf Output:} $q^*( {\bf z} )=\prod_{i=1}^n q^*(z_i)$, where each $q^*(z_i)$ is a univariate truncated normal.
	}
		\vspace{5pt}
	{\bf [c]} Estimate  \eqref{eq10} via Monte Carlo as in \eqref{eq13}.
			\vspace{2pt}
		\label{algo:2}
\end{algorithm*}

To address the above issue, we propose to replace the intractable sampling density $p( {\bf z}  \mid  {\bf y} )$ with a mean--field variational approximation $q^*( {\bf z} )=\prod_{i=1}^n q^*(z_i)$ that factorizes over its marginals $q^*(z_1), \ldots, q^*(z_n)$. In this way, the Monte Carlo estimate for $p(y_{n+1}= 1 \mid  {\bf y} )$ can be obtained by sampling $R$ times from $n$ independent univariate approximate densities $q^*(z_1), \ldots, q^*(z_n)$ instead of the exact but intractable joint density $p( {\bf z}  \mid  {\bf y} )$. Recalling the classical mean--field variational Bayes framework \citep[e.g.,][]{blei2017}, the optimal approximating density $q^*( {\bf z} )$ is the one that minimizes the Kullback--Leibler (\textsc{kl}) divergence $\textsc{kl}[q( {\bf z} ) \| p( {\bf z}  \mid  {\bf y} )]=\mathbb{E}_{q( {\bf z} )}\{\log [q( {\bf z} )/p( {\bf z}  \mid  {\bf y} )]\}$ \citep{kullback1951} to $p( {\bf z}  \mid  {\bf y} )$ among all the densities in the mean--field variational family $\mathcal{Q}=\{q( {\bf z} ): q( {\bf z} )=\prod_{i=1}^n q(z_i) \}$. The solution of such a minimization problem is, typically, not available in closed--form but can be obtained via coordinate ascent variational inference (\textsc{cavi}) algorithms \citep[][]{bishop2006,blei2017} that iteratively minimize the \textsc{kl} with respect to one component $q(z_i)$ at a time, keeping fixed  the others at their most recent estimate $ {\bf q} ^{(t-1)}( {\bf z} _{-i})=[q^{(t)}(z_1), \ldots, q^{(t)}(z_{i-1}), q^{(t-1)}(z_{i+1}), \ldots, q^{(t-1)}(z_n)]$. Recalling \citet{bishop2006}, this is accomplished via the updates
\begin{eqnarray}
q^{(t)}(z_i) \propto \exp\{ \mathbb{E}_{ {\bf q} ^{(t-1)}( {\bf z} _{-i})}[\log[p(z_i \mid  {\bf z} _{-i}, {\bf y} )]]\}, \quad \mbox{for each } i=1, \ldots, n,
\label{eq11}
\end{eqnarray}
at iteration $t$, until convergence. In \eqref{eq11}, the quantity  $p(z_i \mid  {\bf z} _{-i}, {\bf y} )$ denotes the full conditional density of $z_i$. Due to the closure under conditioning property of the multivariate truncated normal \citep{horrace2005}, such a quantity can be derived explicitly from $p( {\bf z}  \mid  {\bf y} )$ in \eqref{eq9} and coincides with the density of a univariate truncated normal. In particular, we can express  $p(z_i \mid  {\bf z} _{-i}, {\bf y} )$ as
\begin{eqnarray}
p(z_i \mid  {\bf z} _{-i}, {\bf y} )\propto\phi\{z_i-[\xi_i+ {\bf H} _{z_i}( {\bf z} _{-i}- {\boldsymbol \xi} _{-i})]; \sigma_{z_i}^{2}\}1[(2y_i-1)z_i>0],
\label{eq12}
\end{eqnarray}
where $ {\boldsymbol \xi} _{-i}$ denotes the vector $ {\boldsymbol \xi} $ without the $i$th entry, $ {\bf H} _{z_i}= {\boldsymbol \Sigma} _{ {\bf z} [i,-i]} {\boldsymbol \Sigma} ^{-1}_{ {\bf z} [-i,-i]}$, and $ \sigma_{z_i}^{2}= {\boldsymbol \Sigma} _{ {\bf z} [i,i]}- {\boldsymbol \Sigma} _{ {\bf z} [i,-i]} {\boldsymbol \Sigma} ^{-1}_{ {\bf z} [-i,-i]} {\boldsymbol \Sigma} _{ {\bf z} [-i,i]}$. The density in \eqref{eq12} has an exponential kernel which is linear in $ {\bf z} _{-i}$ and, hence, replacing the expression for $p(z_i \mid  {\bf z} _{-i}, {\bf y} )$ in the \textsc{cavi} updates reported in \eqref{eq11}, it follows that also $q^{(t)}(z_i)$ has a univariate truncated normal density as in  \eqref{eq12} with $ {\bf z} _{-i}$ replaced by $ {\bf z} _{-i}^{(t-1)}=[\mathbb{E}_{q^{(t)}(z_1)}(z_1), \ldots, \mathbb{E}_{q^{(t)}(z_{i-1})}(z_{i-1}),\mathbb{E}_{q^{(t-1)}(z_{i+1})}(z_{i+1}), \ldots, \mathbb{E}_{q^{(t-1)}(z_{n})}(z_{n})]^{\intercal}$. Each term in $ {\bf z} ^{(t-1)}_{-i}$ is the expected value of a univariate   truncated normal which is available in closed--form, thus producing a simple \textsc{cavi} relying on closed--form updates, as outlined in Algorithm \ref{algo:2}.

Once the optimal univariate truncated normal approximating densities $q^*(z_1), \ldots, q^*(z_n)$ are available, equation \eqref{eq10} can be easily evaluated via Monte Carlo by letting
\begin{eqnarray}
		\hat{p}(y_{n+1}= 1 \mid  {\bf y} )=\frac{1}{R}\sum\nolimits_{r=1}^R\Phi(\mu_{x_{n+1}}{+} \  {\bf H} _{x_{n+1}} [ {\boldsymbol \mu} _{ {\bf X} }+  {\boldsymbol \Sigma} _{ {\bf X} } {\bf z} ^{*(r)}];1+\sigma_{x_{n+1}}^2{+} \  {\bf H} _{x_{n+1}}  {\boldsymbol \Sigma} _{ {\bf X} } {\bf H} ^{\intercal}_{x_{n+1}}), \quad
		\label{eq13}
		\end{eqnarray}
		with $ {\bf z} ^{*(r)}=[z_1^{*(r)}, \ldots, z_n^{*(r)}]^{\intercal}$, for $r=1,\ldots, R$, where each $z_i^{*(r)}$ can be efficiently sampled from the corresponding univariate truncated normal approximating density $q^*(z_i)$, independently for $i=1, \ldots, n$ and $r=1, \ldots, R$. Unlike for the multivariate case, sampling from independent univariate truncated normals can be effectively done in standard statistical softwares, thus avoiding issues in large $n$ settings. As outlined in the simulation studies in Section \ref{sec.3}, such an approximate strategy massively reduces computational times without affecting accuracy.

\section{Simulation Studies}
\label{sec.3}
In this section, we study the gains in accuracy and computational scalability of the methods developed in Section \ref{sec.2} relative to state--of--the--art competitors, which include Monte Carlo inference under the widely--used \textsc{stan} implementation (see \texttt{R} package  \texttt{rstan}) of the Hamiltonian no--u--turn sampler \citep{hoff_2014}, and minimax tilting (see \texttt{R} package  \texttt{TruncatedNormal}) methods \citep{botev_2017} to evaluate the multivariate Gaussian cumulative distribution functions involved in the predictive probability expressed in equation \eqref{eq4}.

\begin{figure}[t]
	\centering
		\includegraphics[scale=.58]{F1_1.pdf}
		\includegraphics[scale=.58]{F1_2.pdf}
	\caption{Simulated probabilities $\Phi[f( {\bf x} )]$ on the $100 \times 100$ grid $\mathcal{G}=\{ {\bf x} =(x_1,x_2): x_1 \in (1/100,2/100, \ldots, 100/100),  x_2 \in (1/100,2/100, \ldots, 100/100)\}$ in the unit square, where $f( {\bf x} )$ is a zero mean \textsc{gp} with squared exponential covariance kernel. White circles denote the 100 unknown locations distributed randomly (left) and on a grid (right), used for prediction.}
	\label{fig:F1}
\vspace{-10pt}
\end{figure}

To evaluate performance in high--dimensional settings, we generate the binary responses on the $100 \times 100$ unit grid $\mathcal{G}=\{ {\bf x} =(x_1,x_2): x_1 \in (1/100,2/100, \ldots, 100/100),  x_2 \in (1/100,2/100, \ldots, 100/100)\}$ with equally--spaced predictors, thus obtaining $n=10{,}000$ non--overlapping configurations. At these locations, we simulate the responses $y_1, \ldots, y_{10{,}000}$ from independent Bernoulli distributions with probability parameters $\Phi[f_0( {\bf x} _1)], \ldots, \Phi[f_0( {\bf x} _{10{,}000})]$ displayed in Figure \ref{fig:F1}, where $ {\bf f} _0( {\bf X} )=[f_0( {\bf x} _1), \ldots, f_0( {\bf x} _{10{,}000})]^{\intercal}$ is a sample from a \textsc{gp} with mean function $m( {\bf x} )=0$ and squared exponential covariance kernel $K( {\bf x} , {\bf x} ')=\exp(-\alpha\cdot \|  {\bf x} - {\bf x} ' \|^2)$, where $\alpha=30$. Such a kernel is frequently used in machine learning, thus we focus on this choice for the simulation study. To assess  performance in estimating the predictive probabilities, we also simulate the probability parameters and the associated binary responses for $100$ out--of--sample units under two scenarios. As outlined in Figure \ref{fig:F1}, the first one relies on randomly distributed locations, whereas the second focuses on a grid structure. We also compare performance in lower--dimensional problems with $n \in \{15^2,25^2,50^2 \}$ obtained by selecting a $n^{1/2} \times n^{1/2}$ sub--grid of $\mathcal{G}$ with equally spaced configurations between $0$ and $1$, along with their associated probability parameters and simulated binary responses.

\begin{table}
\caption{Computational efficiency and predictive performance, at varying sample size $n$, of \textsc{stan} \citep{hoff_2014}, \textsc{tn} \citep{botev_2017}, 	\textsc{tlr} (Section \ref{sec.2.1}) and \textsc{vb} (Section \ref{sec.2.2}), when the units to be predicted are distributed either randomly [random] or on a grid [grid]. \textsc{time}: cost in seconds to compute 100 predictive probabilities. \textsc{mse}: mean squared error between the estimated predictive probabilities and the true ones. \textsc{auc}: area under the \textsc{roc} when predicting the out--of--sample binary responses with the estimated predictive probabilities.}\label{tab:T1}
\addtolength{\tabcolsep}{-1pt}
	\centering
	\resizebox{13cm}{!}{
\begin{tabular}{ll|cccc}
Method& Perfomance measures&$n=250$  & $n=625$ & $n=2{,}500$& $n=10{,}000$ \\
\midrule
		&\textsc{time} {\footnotesize{[seconds]}} & 56{,}273 & --- & ---& ---\\
		&\textsc{mse} {\footnotesize{[random]}} & 0.013 & --- & ---& --- \\
	\textsc{stan}		&\textsc{mse} {\footnotesize{[grid]}} & 0.011 & --- & ---& --- \\
	&\textsc{auc}  {\footnotesize{[random]}} & 0.745 & --- & ---& --- \\
				&\textsc{auc}  {\footnotesize{[grid]}} & 0.820 & --- & ---& --- \\
\midrule
		&\textsc{time}  {\footnotesize{[seconds]}} & 349 & 2{,}339 & --- & ---  \\
		&\textsc{mse}  {\footnotesize{[random]}} & 0.012 & 0.011 & --- & ---  \\
	\textsc{tn}		&\textsc{mse}  {\footnotesize{[grid]}} & 0.010 & 0.013 & --- & ---  \\
	&\textsc{auc}   {\footnotesize{[random]}} & 0.731 & 0.687 & --- & ---  \\
		&\textsc{auc}   {\footnotesize{[grid]}} & 0.834 & 0.776 & --- & ---  \\
\midrule
		&\textsc{time}  {\footnotesize{[seconds]}} & 108 & 402 & 2{,}509& 18{,}709 \\
		&\textsc{mse}  {\footnotesize{[random]}} & 0.013 & 0.016 & 0.004& 0.002 \\
	\textsc{tlr}	&\textsc{mse}  {\footnotesize{[grid]}} & 0.011 & 0.022 & 0.005& 0.001\\
		&\textsc{auc}   {\footnotesize{[random]}} & 0.731 & 0.696 & 0.742& 0.758 \\
		&\textsc{auc}   {\footnotesize{[grid]}} & 0.836 & 0.746 & 0.846& 0.836 \\
\midrule
		&\textsc{time}  {\footnotesize{[seconds]}} & 2 & 6 & 115& 4{,}305 \\
		&\textsc{mse}  {\footnotesize{[random]}} & 0.012 & 0.011 &  0.002& 0.001 \\
	\textsc{vb}	&\textsc{mse}  {\footnotesize{[grid]}} & 0.011 & 0.017 & 0.004& 0.001 \\
	&\textsc{auc}   {\footnotesize{[random]}} & 0.736 & 0.697 & 0.743& 0.771 \\
		&\textsc{auc}   {\footnotesize{[grid]}} & 0.832 & 0.752 & 0.840& 0.828 \\
\bottomrule
\end{tabular}}
\end{table}

Table \ref{tab:T1} summarizes the performance in terms of accuracy and computational scalability of the different methods analyzed, at varying $n$ and under the two different scenarios considered for prediction. In reporting the results, we set a conservative computational budget of one day and display the performance measures only for those routines  with a running time below this computational budget. Monte Carlo inference via \textsc{stan} \citep{hoff_2014} relies on the \texttt{rstan} package applied to model \eqref{eq1} in order to obtain posterior samples from $ {\bf f} ( {\bf X} )$ which are then used to compute the predictive probabilities at the unknown locations via ordinary kriging. Such evaluations rely on $10{,}000$ \textsc{mcmc} samples after a burnin of $10{,}000$, setting $\alpha=30$. In evaluating the performance of minimax tilting (\textsc{tn}) \citep{botev_2017}, we  compute the numerator and the denominator in \eqref{eq4} separately via the \texttt{R} package  \texttt{TruncatedNormal}, using the default settings. Since the denominator is shared among all predictions in \eqref{eq4}, the average cost per prediction under the package  \texttt{TruncatedNormal} is close to that of computing one multivariate Gaussian cumulative distribution function. Equation \eqref{eq4} is also evaluated under the novel methods proposed in Section \ref{sec.2.1} (\textsc{tlr}) and summarized in Algorithm \ref{algo:1}, which can be implemented via simple adaptations of the \texttt{R} package  \texttt{tlrmvnmvt} \citep{cao2020}. In implementing such a routine, we set the block size to $n^{1/2}$, the truncation level to $10^{-4}$ and $R=20{,}000$. In evaluating the predictive probabilities under \textsc{tn} and \textsc{tlr}, we avoid setting $\alpha$ equal to the true value $30$, but instead estimate such a quantity by applying the \texttt{R} packages \texttt{TruncatedNormal} and  \texttt{tlrmvnmvt} to maximize, with respect to $\alpha$, the marginal likelihood in \eqref{eq3} on a grid of $60$ equally spaced $\alpha$ values between $15$ and $45$. Results are comparable, although \texttt{TruncatedNormal} provides slightly less noisy estimates of \eqref{eq3} than \texttt{tlrmvnmvt}, at the cost of a substantially higher running time. The estimate of $\alpha$ provided by  \texttt{tlrmvnmvt} is also used in the implementation of the novel variational strategy (\textsc{vb}) presented in Section \ref{sec.2.2} and summarized in Algorithm \ref{algo:2}. Also in this case we consider $R=20{,}000$ Monte Carlo samples to evaluate \eqref{eq10}. Such values are generated from the optimal univariate truncated normal approximating densities produced by the \textsc{cavi} in Algorithm \ref{algo:2}, which can be implemented via minor adaptations of the source code in the GitHub repository \texttt{Probit-PFMVB} \citep{fasano2019asymptotically}. All computations were run on a 2.5 GHz Intel Core i7 CPU workstation, without multithreading.

As clarified in Table \ref{tab:T1}, the methods proposed in Sections \ref{sec.2.1}--\ref{sec.2.2} notably reduce the running times relative to state--of--the--art competitors, thus making prediction under probit \textsc{gp} computationally feasible in those high--dimensional settings that arise commonly in various applications. This is especially true for the \textsc{vb} solution proposed in Section \ref{sec.2.2} which is orders of magnitude faster than its competitors. Such a notable reduction in running times under \textsc{tlr} and \textsc{vb}  is crucially obtained at almost no costs in terms of  accuracy in out--of--sample prediction (see \textsc{auc}) and in the estimation of the predictive probabilities (see \textsc{mse}), when compared to competitors relying on \textsc{mcmc} samples  from the exact posterior (\textsc{stan}) or on exact evaluation of multivariate Gaussian cumulative distribution functions (\textsc{tn}).

\section{Saudi Arabia Windspeed Application}
\label{sec.4}
We conclude by applying the methods developed in Sections \ref{sec.2.1} and \ref{sec.2.2} to a real--world environmental application aimed at modeling whether the local windspeed exceeds a pre--specified working threshold for energy production in a given region of interest in Saudi Arabia.
Wind turbines for generating electricity typically have two windspeed thresholds, of which the lower controls when the blades of the turbine start to be in motion and the higher indicates if the turbine should be switched off to avoid strong wind damage. Here, the binary response $y_i$ is defined as whether the windspeed at the $i$-th location exceeds the lower threshold, thus allowing production of wind power, which is referred to as the working threshold of wind turbines.
This important application is motivated by the growing   domestic energy consumption in Saudi Arabia and by the attempt to reduce the reliance on fossil fuels, thereby leading to an increasing interest on renewable energy sources, including wind \citep[][]{shaahid2014,chen2018,tagle2019,giani2020}. The effective exploitation of such resources and the careful management of the energy stations require careful modeling and prediction at a fine spatial resolution of whether the local windspeed exceeds or not a given threshold for energy production. As we will discuss in the following, such a fine grid of observations commonly produces a sample size around tens of thousands units. This makes state--of--the--art algorithms for probit \textsc{gp} computationally unfeasible in such studies, thus motivating our scalable solutions presented in  Sections \ref{sec.2.1} and \ref{sec.2.2}.

\begin{figure}[t]
	\centering
		\includegraphics[scale=.57]{F2_1.pdf}
		\includegraphics[scale=.57]{F2_2.pdf}
	\caption{Heatmaps representing the windspeed at $140$ meters high (left) and a binary version $y$ of this measure defining whether the local windspeed is sufficiently high for energy production (dark gray: \textsc{yes}; light gray: \textsc{no}) based on the 4m/s threshold (right) on Jan 21st, 2014. The red area denotes the spatial region that is used for modeling and prediction.}
	\label{fig:F2}
\end{figure}

The windspeed dataset considered in this article is produced by the Weather Research and Forecasting (\textsc{wrf}) model \citep{yip2018statistical}, which constructs the weather system via partial differential equations on the mesoscale and demands strong computation capacity to serve meteorological applications \citep{skamarock2008}. The time resolution of our data is daily and we use the windspeed over the region of north--west Saudi Arabia on January 21st, 2014 for modeling and out--of--sample prediction. Such a region covers the wind farm at Dumat Al Jandal, which is the first wind farm in Saudi Arabia and currently under construction as well as the future smart city of NEOM, a strategic component of the Saudi 2030 Vision, where wind power is expected to be a key energy source. Moreover, the windspeed on January 21st, 2014 has high variability across this region, which makes the out--of--sample prediction task more challenging. As shown in Figures \ref{fig:F2} and \ref{fig:F3} the region under analysis is obtained by intersecting the Saudi Arabia territorial map with the rectangle ranging from \textsc{e}$34^\circ 30^\prime$ to \textsc{e}$43^\circ$ and from \textsc{n}$25^\circ$ to \textsc{n}$32^\circ$. Within this region we consider a fine grid of $n = 9{,}036$ equally--spaced locations $ {\bf x} _i=(x_{i1},x_{i2})^{\intercal}=(\texttt{long}_i,\texttt{lat}_i)^{\intercal}$ at which we monitor whether the windspeed is either above $(y_i=1)$ or below $(y_i=0)$ the working threshold of wind turbines for each $i=1, \ldots, 9{,}036$. Following \citet{chen2018}, such a threshold is set at 4 m/s. As for the simulation study in Section \ref{sec.3}, we monitor predictive performance at $100$ out--of--sample locations displayed in Figure \ref{fig:F3}, which are distributed either randomly or on a grid centered at the Dumat Al Jandal wind farm.

\begin{figure}[t]
	\centering
		\includegraphics[scale=.57]{F3_1.pdf}
		\includegraphics[scale=.57]{F3_2.pdf}
	\caption{For the the spatial region that is used for modeling and prediction, heatmaps  defining whether the local windspeed is sufficiently high for energy production (dark gray: \textsc{yes}; light gray: \textsc{no}) based on the 4m/s threshold on Jan 21st, 2014. Black circles denote the 100 unknown locations distributed randomly (left) and on a grid (right), used for prediction.}
	\label{fig:F3}
\end{figure}

Motivated by the results in the simulation study in Section \ref{sec.3}, we consider a probit \textsc{gp} with zero mean function and squared exponential covariance kernel $K( {\bf x} , {\bf x} ')=\exp(-\alpha\cdot \| {\bf x} - {\bf x} ' \|^2)$, where $\alpha$ is estimated via maximization of the marginal likelihood in \eqref{eq3} evaluated via the  \texttt{tlrmvnmvt} package on a grid of values in $[1,30]$.  The estimated  $\alpha$ is close to $17$, which has an effective range of $0.42$. This relatively short effective range is consistent with the abrupt changes of the binary responses. Recalling the results in Table \ref{tab:T1}, the calculation of the predictive probabilities is only performed under the methods presented in Sections \ref{sec.2.1} (\textsc{tlr}) and \ref{sec.2.2} (\textsc{vb}) since \textsc{stan} and \textsc{tn} would be computationally impractical in such a high--dimensional setting with $n = 9{,}036$. In implementing both methods we set $\alpha=\hat{\alpha}=17$ and consider the same settings as in the simulation study in Section \ref{sec.3}, thus obtaining running times that are comparable to the simulation scenario  with $n=10{,}000$. More specifically,  \textsc{tlr} and \textsc{vb}  require about 168 and 35 seconds per prediction, respectively, thus confirming the feasibility of such novel methods in high--dimensional settings. We shall also emphasize that these running times are mostly affected by the matrix pre--computation operations in Algorithms \ref{algo:1}--\ref{algo:2}, and hence could be carefully reduced via sparse matrix representations or careful algebraic operations. Out--of--sample predictive performance measured via the \textsc{auc} is similarly accurate under both methods, with a slightly increased improvement provided by \textsc{vb}. In particular, the \textsc{auc} for the random scenario is $0.968$ for \textsc{tlr} and $0.971$ for \textsc{vb}, whereas in the grid setting such a measure is  $0.881$ for \textsc{tlr} and $0.906$ for \textsc{vb}.

\section{Discussion}
\label{sec.5}
This article provides novel expressions for the predictive probabilities under probit models with \textsc{gp} priors, relying either on multivariate Gaussian cumulative distribution functions or on functionals of multivariate truncated normals, and proposes scalable computational strategies to evaluate such quantities in common high--dimensional settings, thus covering an important gap in the literature. As highlighted in the simulations studies in Section \ref{sec.3}, these computational gains are notable and do not sacrifice accuracy, thereby allowing tractable prediction under probit \textsc{gp} in applications that were previously computationally impractical, such as the Saudi Arabia windspeed study  in Section \ref{sec.4}.

The above results open up several avenues for future research. For instance, the methods in Section \ref{sec.2} can be adapted to any probit model with a multivariate Gaussian prior for the linear predictor. These include classical Bayesian probit regression, multivariate probit models and general additive representations relying on basis expansions. Extensions to categorical responses under a multinomial probit \textsc{gp} model or to more general priors such as unified skew--normals can also be explored by leveraging results in \citet{Durante2018}, \citet{fasano2020} and \citet{benavoli2020}.

\section*{Acknowledgments}

This publication is based upon work supported by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research (OSR) under Award No: OSR-2018-CRG7-3742.