Gaussian Process KRR Kernel Ridge Regression RKHS Reproducing Kernel Hilbert Space LRR Linear Ridge Regression OPU Optical Processing Unit RFF Random Fourier Features RBF Radial Basis Function DFA Direct Feedback Alignment MSE Mean Squared Error CE Cross Entropy MAC Method of Auxiliary Coordinates POVI Particle Optimization Variational Inference MLP Multilayer Perceptron MCMC Markov Chain Monte-Carlo SGHMC Stochastic Gradient Hamiltonian Monte Carlo

Keywords:

27 3.7 Lowering the Cost of REINFORCE 28 3.8 Experiments 29 3.8.1 Experiments on the UCI datasets 29 3.8.2 Experiments on image data 31 3.9 Discussion 34 4 Bayesian approach to classification 36 4.1 Introduction 36 4.2 Related Work 37 4.3 Variational Bootstrap for Classification Models 38 4.3.1 Considered Models 38 4.3.2 Variational Bootstrap for Neural Network Regression 39 4.3.3 Variational Bootstrap for Random Fourier Features 39 4.3.4 Dirichlet Label Transformation 40 4.3.5 Classification with Variational Bootstrap 41 vii contents viii 4.3.6 Experiments: Toy dataset 42 4.3.7 Experiments: UCI Datasets 43 4.3.8 Discussion of the approach 45 4.4 Improving the Dirichlet-transformation 46 4.4.1 Motivation 46 4.4.2 Cumulative Density Function Transformation 47 4.4.3 Experiments 49 4.4.4 Discussion of the approach 50 Local Gaussian Process regression 52 5.1 Introduction 52 5.2 Related work 53 5.3 Gaussian Processes, Kernel Ridge Regression, and Localization 54 5.3.1 Gaussian Process Regression 54 5.3.2 Locally Smoothed Gaussian Process Regression 55 5.3.3 Local Kernel Ridge Regression 57 5.4 Experiments 58 5.4.1 Toy dataset 58 5.4.2 UCI datasets 59 5.5 Discussion of the method 61 Imposing Functional Priors on Bayesian Neural Networks 62 6.1 Introduction 62 6.2 A review of types of prior in BNN 63 6.3 Imposing Functional Priors on BNNs 64 6.4 Gaussian Process priors 66 6.4.1 Mini-batching 66 6.4.2 Hyper-parameter optimization 68 6.4.3 Classification 68 6.4.4 Connections with Kernel Ridge Regression 68 6.5 Experimental evaluation of MCMC for GP priors 70 6.5.1 Toy regression dataset 70 6.5.2 UCI regression datasets 71 6.5.3 Toy classification dataset 73 6.5.4 UCI classification datasets 73 6.6 Limitations of the approach 74 6.7 Discussion of the MCMC inference for functional priors 75 6.8 Variational Bootstrap for BNNs with functional priors 76 6.8.1 Linear model case 77 6.8.2 Approximation of the Hilbert space norm 78 6.8.3 Emprical evaluation 79 6.8.4 Discussion of the FVBI method 79 Conclusions 80 7.1 Summary of contributions and open questions 80 7.2 Concluding remarks 82

A B S T R A C T

Deep learning played a significant role in establishing machine learning as a must-have instrument in several modern industries and sciences. However, the use of deep learning poses several new challenges. On the one hand, deep learning requires significant computational power for training models and applying them. Modern hardware comes to its limits when it is necessary to provide access to deep learning for millions of users. The enormous energy consumption of such kind of applications poses additional difficulties. Another problem with modern deep learning is its inability to estimate the uncertainty of the predictions, which creates significant obstacles to deploying deep learning in risk-sensitive applications. This thesis presents four projects to address these problems.

• We propose an approach making use of Optical Processing Units for deep learning applications. Optical Processing Units (OPUs) are computing devices that perform random projections of input data by exploiting the physical phenomenon of scattering a light source through a diffusive medium. This operation can be employed to construct kernel ridge regression models and serve as a component in deep neural networks. This novel hardware requires less energy and time to perform computations compared to classical hardware. However, OPUs require the input data to be binary and performs a non-differentiable operation, which poses problems when using this device on non-binary datasets or when integrating it into models trained with backpropagation. We overcome this difficulty by considering OPU as a black box and employing the REINFORCE gradient estimator. Then we integrate the OPU into a deep model, where the first layers serve as a binary encoder. REINFORCE gradient estimator allows us to calculate the gradient of the loss function with respect to the weights of the binary encoder and optimize these together with the parameters of kernel ridge regression with gradient-based optimization Naturally, our approach allows for the integration of OPUs into deep models trained with backpropagation to bring the advantages of this hardware to the domain of deep learning.

• We address the problem of uncertainty estimates for classification with Bayesian inference. We introduce techniques for shallow and deep models that can decrease the cost of Bayesian inference.

Carrying out Bayesian inference over parameters of statistical models is intractable when the likelihood and the prior are non-conjugate. Variational bootstrap provides a way to obtain samples from the posterior distribution over model parameters, where each sample is the solution of a task where the labels are perturbed. For Bayesian linear regression with a Gaussian likelihood, variational bootstrap yields samples from the exact posterior, whereas for nonlinear models with a Gaussian likelihood, some guarantees of approaching the true posterior can be established. In this part of the thesis, we extend variational bootstrap to the Bernoulli likelihood to tackle classification tasks. We use a transformation of the labels, which allows us to turn the classification task into a regression one. Then we apply variational bootstrap to obtain samples from an approximate posterior distribution over the parameters of the model. Variational bootstrap allows us to employ advanced gradient optimization techniques, which provide fast convergence. We provide experimental evidence that the proposed approach allows us to achieve classification accuracy and uncertainty estimation comparable with Markov Chain Monte Carlo (MCMC) methods at a fraction of the cost.

• Classical non-parametric Bayesian linear models, such as Gaussian Processes, are hard to scale to large datasets. We developed a novel framework to accelerate Gaussian process regression (GPR). We considered localization kernels at each data point to down-weigh the contributions from other data points that are far away, and we derived the GPR model stemming from the application of such localization operation. Through a set of experiments, we demonstrated the competitive performance of the proposed approach compared to full GPR, other localized models, and deep Gaussian processes. Crucially, these performances were obtained with considerable speedups compared to standard global GPR due to the sparsification effect of the Gram matrix induced by the localization operation.

• Specifying sensible priors for Bayesian neural networks (BNNs) is key to obtaining good predictive performance and sound predictive uncertainties. However, this is generally difficult because of the complex way prior distributions induce distributions over the functions that BNNs can represent. Switching the focus from the prior over the weights to such functional priors allows for the reasoning on what meaningful prior information should be incorporated. We propose enforcing such meaningful functional priors through Gaussian processes (GPs), which we view as a form of implicit prior over the weights. We employ scalable Markov chain Monte Carlo (MCMC) to obtain samples from the posterior distribution over BNN weights. Unlike previous approaches, our proposal does not require the modification of the original BNN model, it does not require any expensive preliminary optimization, and it can use any MCMC techniques and any functional prior that can be expressed in closed form. We illustrate the effectiveness of our approach with an extensive experimental campaign.

A C K N O W L E D G M E N T S

During my PhD journey, I was fortunate to have met and been influenced by many remarkable people who have inspired and supported me in my work. I am now grateful for the opportunity to express my gratitude to everyone who made this journey possible. First and foremost, I would like to express my immense gratitude to Maurizio Filippone for his unrelenting support, guidance, and patience throughout my PhD. I am genuinely thankful for the countless hours he spent reviewing drafts, providing suggestions, and offering constructive criticism.

I would also like to extend my gratitude to the 3IA Interdisciplinary Institute for Artificial Intelligence, and through them, to the French government for their confidence and financial support during my PhD.

My sincere thanks go to EURECOM for providing the administrative support and infrastructure that made my research possible. I am grateful for the research-friendly environment they provided, which made my dissertation work an enjoyable experience.

I am thankful to all of my colleagues, especially those I worked closely with. In particular, I am grateful to Dimitrios Milios for the numerous inspiring discussions and advice that allowed many of the ideas presented in this work to emerge. I want to thank Davit Gogolashvili for his valuable ideas and collaboration on Chapter 5 of this work. I am also grateful to Jonas Wacker and Simone Rossi for answering my endless questions and introducing me to the environment. Finally, I want to express my gratitude to Ba-Hien Tran for the daily conversations that provided a much-needed break from the intensity of research. The plot of the activation functions g -1 1 , g -1 2 for µ = 12.74, σ = 10 48 Figure 4. 5 The histograms of predictions in the latent space (right) and corresponding transformed values in the space of class probabilities (left). The value of the quadratic regularization term with respect to the number of points used for evaluation of the function. Average classification accuracy for UCI classification datasets 74 Table 6. 4 Average test NLL for UCI classification datasets 74

C O N T E N T S

x Table 6.5

RMSE for GP regression (GPR), the proposed method (FBVI) and WDGPi-G 79

L I S T I N G S A C R O N Y M S 1 I N T R O D U C T I O N
Recent advances in machine learning allowed to completely reshape a lot of fields of technology, starting from accessible consumer-level computer vision systems [START_REF] He | Deep Residual Learning for Image Recognition[END_REF], dialog [START_REF] Openai | [END_REF], translation systems [START_REF] Vaswani | Attention is All you Need[END_REF] to complicated systems for controlling and improving production [START_REF] Li | Deep learning for smart industry: Efficient manufacture inspection system with fog computing[END_REF], and financial risk assessment [START_REF] Nhi | Deep learning for decision making and the optimization of socially responsible investments and portfolio[END_REF]. In addition, modern machine learning techniques are used to drive scientific discoveries in various fields like experimental physics [START_REF] Derkach | Cherenkov Detectors Fast Simulation Using Neural Networks[END_REF] and mathematics [START_REF] Lample | Deep Learning For Symbolic Mathematics[END_REF], medicine [START_REF] Khalid | Skin Cancer Classification using Deep Learning and Transfer Learning[END_REF], chemistry [3], and many others [START_REF] Raghu | A survey of deep learning for scientific discovery[END_REF].

One of the fascinating things about these advances is that machine learning as a field existed for a long time, but its applications were scarce. We can track the origins of machine learning to the middle of the 20th century. Arthur Samuel was one of the first who started to use the term "machine learning" [START_REF] Samuel | Some Studies in Machine Learning Using the Game of Checkers[END_REF]. His work in the late 40s introduces a self-learning computer program for playing checkers and introduces the first machine learning algorithms. Later, Rosenblatt combined Samuel's ideas with the model of brain cell interaction, proposed by [START_REF] Donald | The first stage of perception: growth of the assembly[END_REF], to introduce a perceptron model [START_REF] Rosenblatt | The perceptron: a probabilistic model for information storage and organization in the brain[END_REF]. At that time, this method mainly had theoretical value as a step to understanding the principles of information processing in the brain. Still, the practical applicability of the method was limited to some image processing problems [START_REF] Albert E Murray | Perceptron applications in photo interpretation[END_REF]. At the time, it seemed that the future of machine learning was in more sophisticated techniques from mathematical statistics.

One such statistical method that impacted machine learning was the classification algorithm proposed by [START_REF] Fisher | The use of multiple measurements in taxonomic problems[END_REF]. It played a notable role in applied machine learning in the second half of the 20th century for developing various computer vision systems [START_REF] Turk | Eigenfaces for recognition[END_REF], [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF].

In the 70s, Vapnik [START_REF] Vladimir | Estimation of Dependencies Based on Empirical Data[END_REF] and Chervonenkis [START_REF] Vapnik | Theory of pattern recognition[END_REF] developed a powerful theoretical foundation for a new type of machine learning algorithms called Support Vector Machines. These algorithms proved themselves very powerful and successful in regression, classification, and clustering problems. Among the aforementioned methods, SVMs probably have had the most significant impact in industry at the time and anticipated the deep learning revolution. They were used for natural language processing [START_REF] Joachim | Text categorization with support vector machines[END_REF], various tasks in computer vision [START_REF] Maity | Supervised classification of RADARSAT-2 polarimetric data for different land features[END_REF], problems in biology [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF], and other sciences [START_REF] Hastie | Discussion of "Support vector machines with applications" by Javier Moguerza and Alberto Munoz[END_REF].

The vanilla versions of the aforementioned methods are linear with respect to input features. Thus, their area of application was quite restricted. They owe their success on a practical level to the so-called kernel trick [START_REF] Gogolashvili | Locally Smoothed Gaussian Process Regression[END_REF], wich allowed SVMs [START_REF] Boser | A Training Algorithm for Optimal Margin Classifiers[END_REF], and other classical statistical methods, like Ridge Regression [START_REF] Schölkopf | Nonlinear Component Analysis as a Kernel Eigenvalue Problem[END_REF], to tackle complex non-linear problems. At the same time, kernel methods made possible the application of Gaussian Processes, proposed by [START_REF] Kac | An explicit representation of a stationary Gaussian process[END_REF] to tackle machine learning problems [START_REF] Opper | Gaussian processes for classification: Mean-field algorithms[END_REF], giving the linear models not only the flexibility to resolve non-linear problems but also an ability to estimate the uncertainty of the model. Thus, we can conclude that the advances of practical machine learning in the second half of the 20th century were mainly based on the solid foundations of mathematical statistics and their introduction to the Computer Science community. This strong basis enabled accurate predictions but also a theoretical analysis of model's properties, like statistical guarantees for the risk on unseen data [START_REF] Caponnetto | Optimal Rates for the Regularized Least-Squares Algorithm[END_REF], model interpretability [START_REF] Statnikov | Using SVM weight-based methods to identify causally relevant and non-causally relevant variables[END_REF], and uncertainty quantification [START_REF] Edward Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF]. The overwhelming success of deep learning at the beginning of the 21st century resulted from a completely different approach based on vast amounts of data, computational resources, and the simple Perceptron model proposed half a century ago. Extended with gradient descent optimization and backpropagation to train deeper models, it enabled striking versatility. Deep learning reduces pipeline creation time by automatizing the burden of feature engineering due to its ability to train machine learning pipelines in an "end-to-end" fashion.

Deep neural networks can solve much more complicated problems, but with a cost. Theoretical guarantees of deep models are an open question in machine learning, except for some exceptional cases [START_REF] Jacot | Neural Tangent Kernel: Convergence and Generalization in Neural Networks[END_REF]. Predictions of these models tend to be overconfident [START_REF] Guo | On Calibration of Modern Neural Networks[END_REF] and hard to interpret. The community is constantly looking for effective solutions to counter these problems, but cheap and reliable solutions are still lacking. Nowadays, it may seem that most of the advances of machine learning in the 20th century did not play any role in the current state of the art.

This work proposes four separate projects connected by the idea of reintroducing some advances of classical machine learning and merging them with the power of deep neural networks. In this chapter, we will discuss the motivation and contributions of each project.

machine learning at the speed of light

As mentioned above, the success of deep learning became possible due to significant advances in hardware. General purpose CPUs are not very efficient in matrix and tensor operations, which are a building block of all algorithms for inference and training of neural networks. Nowadays, GPUs and TPUs are absolutely necessary for research and practical applications of deep learning. The problem is that even this specialized hardware can face performance problems with applications that require high loads, such as requests from a significant amount of users or applications that require low latency.

One of the possible solutions to this problem is to use a completely different approach to perform computations. In Chapter 3 we consider a device that uses light to perform computations. It projects an input vector into a high-dimensional feature space using random projections. Then, it is possible to use these features to solve a linear regression or classification problem. From the mathematical point of view, this device allows solving a Kernel Ridge Regression (KRR) problem with a particular type of kernel: the OPU-kernel.

Utilization of the considered device for practical applications poses two main challenges. The first problem appears due to the fact that Optical Processing Units (OPUs) can work only with binary input vectors, and most real-life datasets, considered by machine learning, are not binary. It means that machine learning pipelines containing an OPU should also include a binarization mechanism for input data. The second challenge, posed by the OPU architecture, is due to the fact that OPUs performs an operation that is non-differentiable. This means that it is not possible to integrate OPUs into the deep learning pipelines for "end-to-end" training using standard backpropagation. This situation becomes especially problematic when it is necessary to work with complicated feature extraction procedures that should be trained with the rest of the pipeline, like, for example, image data, that requires using convolutions to extract shift-invariant features.

Chapter 3 is dedicated to creating a training procedure for resolving these problems. The main contributions of Chapter 3 are:

• a training method that allows building heterogeneous pipelines with the OPU and training them in an "end-to-end" fashion;

• an empirical demonstration that our approach allows obtaining binary embeddings that improve the performance of the pipeline on classification or regression problems compared to non-supervised binarization techniques;

• an empirical demonstration that the proposed procedure allows integrating the OPU as a part of complex deep learning pipelines to address machine learning problems that require learnable feature extraction procedures.

performing classification in a bayesian way

One of the common pitfalls of modern deep learning is overconfidence in predictions of deep models. Classification models take a vector as input and return a probability of belonging of this input vector to classes defined by the problem. But in practice, deep models assign a class label with high probability, even for input vectors that lie far away from data distribution and do not belong to any class of interest. We can consider a neural network trained to classify pictures of dogs and cats as a simple example. If we feed the model an image of a car, it will assign it to one of the classes with a probability close to 1, while the desired behavior would be 0.5 probability for each of the two classes. The problem of overconfidence in predictions can be addressed with a Bayesian approach. Instead of producing a single prediction, Bayesian Neural Network produces the distribution of class probabilities or samples from this distribution. Analyzing the entropy of this distribution, it is possible to reason about the model's confidence level.

However, this property of Bayesian neural networks comes with a cost. There are multiple approaches to performing Bayesian inference for Bayesian Neural Networks, such as Variational Inference (VI) or Markov Chain Monte-Carlo (MCMC). VI approaches, in most cases, are unable to recover the true posterior predictive

We consider VI and MCMC approaches in detail in Chapter 2 distributions, which is why they provide inadequate estimates of the uncertainty of the predictions. MCMC approaches provide good uncertainty estimates, but they require a lot of time to converge. Thus the main challenges for creating an inference procedure are obtaining reliable uncertainty estimates of the predictions and decreasing the time necessary for the inference procedure.

Chapter 4 is dedicated to addressing these challenges. In summary, in Chapter 4 we make the following contributions:

• we propose a novel algorithm for performing Bayesian inference for classification models;

• we demonstrated its applicability both for linear models and deep neural networks;

• we empirically demonstrated that the proposed method achieves a faster convergence rate than the MCMC approach while preserving comparable quality uncertainty estimations.

gaussian processes in the age of big data

Nowadays, it is common to rely on big amounts of data to train statistical models.

As mentioned, the success of deep learning owes much to a large amount of training data. And deep models can benefit from the big amount of data because of the vast amount of parameters that they have. If the big amount of data and parameters is the key to the machine learning method's success, can we extend the same principle to classical machine learning algorithms? Classical machine learning contains a category of methods based on this principle called non-parametric methods. This category includes many methods, like the aforementioned Gaussian Processes and other types of kernel machines. The general idea of all these methods is to grow the capacity of the model with the size of the training dataset. The problem is that most of these methods are not adapted to huge datasets of the deep learning era.

In Chapter 5 we consider one such method: Gaussian Process Regression. This algorithm possesses several nice properties, like flexibility, thanks to the use of kernels, which allows encoding various constraints on a modeling function, like smoothness, trend, or seasonality. In addition to this, Gaussian Processs (GPs) are Bayesian methods, which means, that they can provide information about the uncertainty of the predictions, and they have a natural mechanism of taking into account the information about measurement noise in training labels. However, classical GPs cannot deal with large data sets even using modern hardware. Another problem with this algorithm is that it uses that for most of the popular kernel functions, like Radial Basis Function (RBF) or Matérn, training points that lie far away from the input point of interest do not affect the prediction, while affecting the computational complexity of the algorithm.

We discuss GPs in detail in Section 2.2.1

There were multiple attempts to make GP scalable to large data sets, but most of these methods are based on distilling or synthesizing some small subset from training data and then making all predictions based on this small subset. We

We consider different approaches to build scalable GPs at Section 5.2 believe that this approach contradicts the philosophy that made deep learning successful. Indeed, modern neural networks very often have a larger number of parameters than the number of training samples. Thus, they do not perform a shrinking of any sort. In Chapter 5, we propose a method that extends the Gaussian Process Regression algorithm that is able to perform at scale while keeping all the training data available for predictions. To obtain this property, we leveraged the idea of localization. For each test data point, we use only the training points relevant to this specific region of space. From one point of view, that is how we avoid wasting computational resources on training points that do not affect this specific prediction. From another point of view, the test set, with points located in different regions of input space, is capable of leveraging the full training set.

The contributions of Chapter 5 are:

• the approach that allows performing Gaussian Process Regression at scale while effectively using computational resources and exploiting the full size of the training set;

• the empirical demonstration of the superiority of the proposed approach above classical GP Regression and other localization and sparsification techniques in terms of computational complexity and prediction performance.

how gaussian processes can improve bayesian neural networks

Bayesian neural networks have the same advantages as GP, like prediction uncertainty estimation and the ability to account for noise in training labels, but in addition they possess the expressiveness of deep models. The problem is that the Bayesian approach requires choosing the prior over the variable of interest. In the case of Bayesian Neural Network (BNN) it is the vector of neural network weights.

There are a lot of approaches to impose different kinds of prior distributions over the weights. Still, most of them cannot address the main challenge of imposing meaningful functional properties on the model's output by imposing prior distributions over the model's weights. The described situation is a particular case of

We discuss in detail priors for BNNs in Section 6.2 the problem, general for most deep learning techniques: lack of interpretability. But in the particular case of BNNs, it works the other way around. In standard deep learning, the functional properties of the model reflect the properties of a training set, but it is hard to interpret the intrinsic properties and structure of the weight. In BNNs, it is hard to interpret how the properties of the weights imposed by the prior affect the model's output.

In Chapter 6, we proposed a solution for this problem. Instead of manipulating the model's weights by imposing a prior, we decided to put prior distribution over the model's output. For this purpose, we need a distribution over functions. Such distribution can be provided by a GP. The idea of imposing the prior over BNNs is not novel, but we managed to propose a generalized approach that uses MCMC, standard for the BNN inference without any additional approximations. While our method does not provide exact posterior inference, we demonstrated empirically that it is competitive with other similar techniques. The main contributions of Chapter 6 are:

• A scalable approach for the inference of BNNs with a GP prior using Markov Chain Mote-Carlo

• Empirical demonstration of the effectiveness of the proposed approach compared to other methods of imposing a functional prior.

B A C K G R O U N D
2.1 linear methods in machine learning

Linear Ridge Regression

Let X = x 1 , . . . , x n be a set of input vectors x ∈ R d and let y = y 1 , . . . , y n be a set of labels associated with the input vectors. Let's model the connection between the output y and the input vector x in the following way:

y i = f (x i , w) + ε, ε ∼ N (0, σ noise ), (2.1)
where f (ϕ(x), w) is a deterministic function with parameters w and ϕ(•) is some function that extracts features from x. In this section, we will consider f to be a linear function

f (x, w) = w ⊤ ϕ(x). (2.2)
Given the Gaussian assumption about the noise ε, we can see that the likelihood of observing a response y given the input x is Gaussian.

p(y|x) = N (y| f (ϕ(x), w), σ noise ). (2.3)
Now we can find a maximum likelihood estimate of the parameters of the model w.

It is easier to minimize a log-likelihood instead of maximizing the likelihood itself.

We will obtain an equivalent solution because the logarithm is a monotonically increasing function

ŵ = argmin w 1 2σ 2 noise n ∑ i=1 (y i -w ⊤ ϕ(x i )) 2 + const. (2.4)
It is possible to introduce an L2 regularization term to the optimization objective that will penalize large values in w. Thus, we can write the optimization objective of Linear Ridge Regression (LRR):

ŵ = argmin w 1 2 n ∑ i=1 (y i -w ⊤ ϕ(x i )) 2 + λ 2 ∥w∥ 2 2 , (2.5)
where parameter λ controls the level of regularization. In this formulation, it reflects the assumption regarding the variance of the observation noise λ = σ 2 noise .

Kernel Ridge Regression

Now we can consider a KRR algorithm, which is an extension of LRR. Kernel Ridge Regression (KRR) is a statistical model which constructs a functional relationship between the inputs and the labels which belongs to the so-called Reproducing Kernel Hilbert Space (RKHS). The properties of such functions, such as smoothness, are characterized by choice of a so-called kernel function k(•, •) : R d × R d → R [START_REF] Murphy | Machine learning: a probabilistic perspective[END_REF], which is a positive semi-definite function of pairs of input points returning a scalar. The reproducing property of kernel functions is ⟨k(x, •), k(y, •)⟩ = k(x, y).

According to Mercer theorem, positive definiteness of kernel functions implies that we can express k(x i , x j ) = ϕ(x i ) ⊤ ϕ(x j ) for some set of (possibly infinite) basis functions ϕ(•).

In order to derive the conventional formulation of kernel ridge regression, it is useful to start from linear regression, where a set of model parameters w is introduced to express a linear relationship between inputs and labels. Then, one introduces the following optimization problem:

ŵ = argmin w 1 2 n ∑ i=1 (y i -w ⊤ x i ) 2 + λ 2 ||w|| 2 2 .
(2.6)

The objective function contains two terms; the first is a model fitting term, while the second is a regularization term, which prevents the weights from becoming too large. The solution to this optimization problem is available in closed form, given that the objective is quadratic with respect to the parameters, yielding: ŵ = (X ⊤ X + λI) -1 y.

(2.7)

Using standard algebraic manipulations involving the Woodbury identity, we can re-express the solution as:

ŵ = X ⊤ (XX ⊤ + λI) -1 y. (2.8)
While this is more costly than the previous expression in the common case where d < n (inversion of a n × n matrix rather than a d × d matrix), this formulation is useful to derive kernel ridge regression. Imagining to introduce basis functions ϕ(•) = (ϕ 1 (•), . . . , ϕ D (•)) ⊤ , we can solve this new optimization problem

ŵ = argmin w 1 2 n ∑ i=1 (y i -wϕ(x i )) + λ 2 ||w|| 2 , (2.9)
with solution

ŵ = Φ ⊤ (ΦΦ ⊤ + λI) -1 y. (2.10)
Evaluating the model at a given input x * yields:

ϕ(x * ) ⊤ ŵ = ϕ(x * ) ⊤ Φ ⊤ (ΦΦ ⊤ + λI) -1 y. (2.11)
In this expression, we recognize the scalar product of vectors of basis functions.

Then we can express these scalar products as a kernel function and obtain the following:

ϕ(x * ) ⊤ ŵ = k * (K + λI) -1 y, (2.12)
where

k * = [k(x 1 , x * ), . . . , k(x n , x * )] ⊤ and K ij = k(x i , x j ).
In practice, one first chooses a kernel function, and this induces a set of basis functions; the beauty of this formulation is that one never explicitly works with the set of basis functions, and all we need to use this model in practice is the evaluation of kernel functions among inputs.

The same method can be approached from a functional space point of view. In this case, it is necessary to solve an optimization problem in RKHS imposed by the kernel function k(•, •) with respect to a regularization constraint over the RKHS norm of the function.

f * (•) = argmin f ∈H 1 2σ 2 n n ∑ i=1 (y i -f (x i )) 2 + 1 2 || f || 2 H . (2.13)
The Representer theorem [START_REF] Mohri | Foundations of Machine Learning[END_REF] shows that the solution to this problem has the following form:

f * (•) = n ∑ i=1 α i k(•, x i ). (2.14)
We can find an expression for f * (•) by substituting this result into (2.13) and solving it with respect to α

f * (x * ) = k * (K + σ 2 n I) -1 y, (2.15)
which is equivalent to the expression (2.13) for σ 2 n = λ. The choice of the kernel for the KRR algorithm strongly depends on the specific problem. The most popular kernel is a RBF kernel, also known as the Gaussian kernel:

k(x, x ′ ) = α exp(- ∥x -x ′ ∥ 2 l
).

(2.16)

This kernel imposes a smoothness assumption on the model f (•), and it has two parameters α, l. Using these parameters, it is possible to control the amplitude and smoothness of the modeling function f (•).

Another popular choice for the kernel function is the Polynomial Kernel:

k(x, x ′ ) = (x ⊤ x ′ + c) d . (2.17)
This kernel enables modeling of linear and non-linear relationships between input and output variables, depending on the degree d. In fact, it maps the input data into a higher-dimensional space, and it is particularly useful when the amplitude of the input vector x has significance for the task, as is often the case with computer vision tasks, for example.

Random Features

One of the main limitations of kernel methods is the scalability to large datasets. The problem arises from the need to evaluate and perform algebraic operations with the so-called Gram matrix K. Because K is an n × n matrix, evaluating and storing K requires O(n 2 ) computations and storage, while any algebraic operations, such as factorization to handle the inverse of K + λI, requires O(n 3 ) operations. These prevent the applicability of kernel methods in their exact form to datasets of size beyond a few thousand. It is worth noting that some approaches have been proposed to solve algebraic operations in an iterative fashion and without the need to store K [START_REF] Filippone | Enabling scalable stochastic gradient-based inference for Gaussian processes by employing the Unbiased LInear System SolvEr (ULISSE)[END_REF], [START_REF] Cutajar | Preconditioning Kernel Matrices[END_REF], [START_REF] Wang | Exact Gaussian Processes on a Million Data Points[END_REF]. However, they still require O(n 2 ) computations for each iteration of their solvers. Furthermore, while the number of iterations of the solvers is much lower than n in practice, in the worst case, it can be O(n), leading to a worst-case complexity of O(n 3 ).

The literature offers a number of solutions to scale kernel methods to large data linearly in the number of data, such as Nyström approximations [START_REF] Fowlkes | Efficient spatiotemporal grouping using the Nystrom method[END_REF] and random features [START_REF] Rahimi | Weighted Sums of Random Kitchen Sinks: Replacing minimization with randomization in learning[END_REF]. In this section, we consider the random features approach and in Chapter 3 we will consider a hardware implementation of random features in Optical Processing Units (OPUs) [START_REF] Saade | Random projections through multiple optical scattering: Approximating Kernels at the speed of light[END_REF]. Random feature approximations form a class of approximations that attempt to construct a finite set of basis functions

ϕ(•) ∈ R D such that k(x i , x j ) ≈ ϕ(x i ) ⊤ ϕ(x j ).
(

2.18)

There are different ways to construct such sets of basis functions, depending on the kernel. For example, so-called random Fourier features are commonly employed to approximate the Gaussian kernel:

k(x i , x j ) = exp(-||x i -x j || 2 ). (2.19)
Appealing to Bochner's theorem [START_REF] Rahimi | Weighted Sums of Random Kitchen Sinks: Replacing minimization with randomization in learning[END_REF], this kernel, which is shift-invariant due to the dependence on τ = x ix j , admits an alternative expression as:

k(τ) = p(ω) exp(i2πωτ)dω, (2.20)
where p(ω) is a proper density function and i = √ -1. Interpreting this as an expectation under p(ω), it is possible to approximate the integral as an expectation using Monte Carlo:

k(τ) = 1 D ∑ r exp(i2πω(r)τ), (2.21)
with ω (r) ∼ p(ω). Furthermore, it is possible to use simple trigonometric identities to verify that the complex exponential can be broken down as a scalar product with terms depending on x i and x j :

k(x i , x j ) = 1 D ϕ(x i ) ⊤ ϕ(x j ), (2.22)
with ϕ r (x) = (sin(x ⊤ ω (r) ), cos(x ⊤ ω (r) )).

(2.23)

We refer the reader to [START_REF] Rahimi | Weighted Sums of Random Kitchen Sinks: Replacing minimization with randomization in learning[END_REF], [START_REF] Cho | Kernel Methods for Deep Learning[END_REF], [START_REF] Cutajar | Random Feature Expansions for Deep Gaussian Processes[END_REF], [START_REF] Wacker | Improved Random Features for Dot Product Kernels[END_REF] for random features derived from alternative integral representation to the Fourier transform.

bayesian methods

In classification and regression tasks, the objective is to tune the parameters w of a function f (x, w) so that it best explains a set of given observations. The Bayesian paradigm dictates that we specify a prior distribution p(w), which captures any prior knowledge about the model. For a dataset D{(x i , y i ) | i = 1 . . . n} of n input vectors x i and labels y i , it is possible to define the posterior distribution over w after observing the data D by means of Bayes rule:

p(w|D) = p(D|w)p(w) p(D|w)p(w)dw , ( 2.24)
where p(D|w) is known as the likelihood of the model, while the integral term in the denominator above is referred to as marginal likelihood or evidence.

Bayesian Linear Models for Regression

bayesian linear regression We can consider LRR from a Bayesian perspective. If we assume that the prior distribution over the parameters N (0, σ prior ) is Gaussian, we can easily verify that the maximum a posteriori estimate of the parameters w corresponds to the LRR optimization objective.

ŵ = argmin w 1 2 n ∑ i=1 (y i -w ⊤ ϕ(x i )) 2 + 1 2 σ 2 noise σ 2 prior ∥w∥ 2 2 .
(2.25)

But for the linear regression case, it is possible to go further. Due to the fact, that the Gaussian prior is conjugate to the Gaussian likelihood [START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF] it is possible to solve the integral in the denominator of (2.24) analytically and obtain an analytical expression for the posterior distribution of the parameters w. After observing n data points, the posterior distribution of w will be Gaussian with parameters:

µ w = σ -2 noise Σ w Φ ⊤ y, Σ -1 w = σ -2 prior I + σ -2 noise Φ ⊤ Φ.
(2.26)

Having the posterior distribution over the weights it is possible to write the posterior predictive distribution given the test point

x * p(y * |x * , X, y) = N (y * |µ ⊤ w ϕ(x * ), σ y * ), σ 2 y * = σ 2 noise + ϕ(x * ) ⊤ Σ w ϕ(x * ).
(2.27)

gaussian process A GP regression can be seen as a Bayesian equivalent of the KRR method. And as the KRR method, GP regression can be approached from the parametric or the functional space points of view. Let's consider a parametric approach first. We can re-express (2.27) using the kernel k(x, x ′ ) = σ 2 w ϕ(x) ⊤ ϕ(x ′ ). In this case, the mean of the predictive distribution for the input x * is given by the following expression:

µ ⊤ w ϕ(x * ) =(σ -2 noise (σ -2 prior I + σ -2 noise Φ ⊤ Φ) -1 Φ ⊤ y) ⊤ ϕ(x * ) = σ 2 prior ϕ(x * ) ⊤ Φ ⊤ (σ 2 noise I + σ 2 prior ΦΦ ⊤ ) -1 y = k * (K + σ 2 noise I) -1 y.
(2.28)

And the variance of the predictive distribution becomes:

σ 2 y * = σ 2 noise + ϕ(x * ) ⊤ (σ 2 noise I + σ 2 prior Φ ⊤ Φ) -1 ϕ(x * ) = σ 2 noise + ϕ(x * ) ⊤ (σ 2 prior I -σ 2 prior Φ ⊤ (σ 2 noise I + σ 2 prior ΦΦ ⊤ ) -1 σ 2 prior Φ)ϕ(x * ) = σ 2 noise + k(x * , x * ) -k ⊤ * (K + σ 2 noise I) -1 k * .
(2.29)

From the functional space point of view, GP is a distribution over functions.

Definition 1 A Gaussian Process is a collection of random variables, any finite number of which have a joint Gaussian distribution [START_REF] Edward Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF].

The distribution of a finite set of these variables is determined by the mean function m(x) and the kernel function k(•, •):

m(x) = E[ f (x)], k(x, x ′ ) = E[ f (x) f (x ′ )]. (2.30)
Let's consider a GP with a constant mean function m(x) = 0. We can denote the GP as

f (x) ∼ GP (0, k(x, x ′ )). (2.31)
To apply this concept to regression problems, we consider the model

y i = f (x i ) + ε, ε ∼ N (0, σ noise ). (2.32)
We can choose a GP with m(x) = 0 as a prior distribution over f (x). Given the assumption about the Gaussainity of the noise case, the prior distribution over the outputs is also Gaussian.

y ∼ N (0, K + σ 2 noise I), where K ij = k(x i , x j ).

(2.33)

Now, given the dataset D{(x i , y i ) | i = 1 . . . n}, it is possible to obtain a posterior distribution over f (•) at the test point x * . For a non-parametric model such as GP the posterior distribution plays a role of a predictive distribution of a parametric model. The posterior distribution over f (•) is also a GP:

f (x)|X, y ∼ GP ( m(x), k(x, x ′ )). (2.34)
Using the expression for a conditional Gaussian distribution [START_REF] Edward Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF], we can obtain the functions m(x) and k(x, x ′ ):

m(x) = k K + σ 2 noise I -1 y, k x, x ′ = k x, x ′ -k K + σ 2 noise I -1 k ′ , (2.35)
where

k ⊤ = [k(x, x 1 ), . . . , k(x, x n )] and k ′⊤ = [k(x ′ , x 1 ), . . . , k(x ′ , x n )].
M MAP = argmax M p(D|M)p(M) = argmax M p(D|M).
(2.37)

In this context, let's consider the models discussed in Section 2.2.1. The performance of the Bayesian method depends on the choice of the noise parameter σ noise and the prior parameter σ prior . We have shown their connection to the regularization parameter λ that controls regularization in frequentist linear methods in (2.25). Thus, we can conclude that an incorrect choice of these parameters will lead to the same problems typical for the incorrect choice of regularization. If σ prior is too big, it will lead to model overfitting. If this parameter is too small, we will underfit the data. The opposite is true for σ noise . Thus, we can claim that different values of these hyperparameters define different models, and we can use Bayesian model selection to optimize these parameters. Thus, we can write p(D|σ noise , σ prior ) = p(D|M).

To find optimal values for these parameters, we can consider the evidence of the model, i.e., the denominator of (2.24). p(D|σ noise , σ prior ) = p(D|w, σ noise )p(w|σ prior )dw.

(2.38)

We can solve the optimization problem (2.37) with respect to the parameters σ noise and σ prior . For linear models with a Gaussian prior and likelihood, we can obtain an analytical expression for (2.38) using the property for the conditional distribution in a linear-Gaussian model. According to this property, the distribution p(D|σ noise , σ prior ) = p(y|σ noise , σ prior ) will also be Gaussian. For simplicity, we consider the prior p(w) to be zero-centered.

p(y|σ noise , σ prior ) = N (y|0, σ 2 noise I + σ 2 prior ΦΦ ⊤ ).

(2.39)

Instead of maximizing the expression for p(y|σ noise , σ prior ) directly, we can maximize its logarithm. After some simple transformations, we can show that:

log p(y|σ noise , σ prior ) = - n 2 log σ 2 noise - d 2 log σ 2 prior - 1 2σ 2 noise ∥y -Φµ w ∥ 2 - 1 2σ 2 prior µ ⊤ w µ w - 1 2 log |Σ w | - n 2 log(2π), (2.40)
where µ w and Σ w are defined in (2.26). The expression for log p(y|σ noise , σ prior ) can be maximized with respect to σ 2 prior and σ 2 noise using a gradient-based optimizer to obtain the optimal values for these parameters. The main difficulty here is the computation of µ w due to the matrix inversion in (2.26), which has O(d 3 ) complexity. Given that we have to perform this computation for each step of the gradient-based optimizer, the overall complexity of obtaining the optimal values for σ 2 prior , σ 2 noise is quite high. We can use the same approach for optimizing the hyperparameters of GP regression. We can use the expression (2.39) directly by using the representation of the Gram matrix through basis functions K = ΦΦ ⊤ . The logarithm of the model evidence in this case becomes:

log p(y) = - 1 2 log |(K + σ 2 noise I)| - 1 2 y ⊤ (K + σ noise I) -1 y - n 2 log(2π).
(2.41)

The computational complexity for this expression is O(n 3 ). We can optimize log p(y) with respect to the noise parameter σ 2 noise and the parameters of the kernel k(•, •), such as the lengthscale or output variance in (2.16).

Bayesian Linear Models for Classification

To address classification problems in a Bayesian way, we have to use another type of likelihood p(D|w) in the expression (2.24). The standard choice is categorical distribution:

p(y|w) = n ∏ i 1[y i = j]a ij , (2.42)
where

a ij = exp(w ⊤ j x i ) ∑ c k=1 exp(w ⊤ k x i )
for c-class problem, and w j is a subset of w, that corresponds to the output j. With this type of likelihood and a Gaussian prior p(w), we can no longer leverage the property of conjugate distributions.

The easiest way to bypass this difficulty is to use the Laplace approximation. We can approximate the true posterior p(w|D) with an approximate posterior distribution q(w|D). As the approximate posterior q(wvect|D), we can choose a Gaussian distribution.

q(w|D) = N (w MAP , Σ Laplace ), (2.43)
where Σ Laplace = -(∇∇ log p(w|D)) -1 and w MAP is the mode of the posterior distribution p(w|D). The idea of this approach is simple; we match the mode of the approximate distribution with the mode of the true posterior distribution. The Hessian of the logposterior arises from the second-order Taylor approximation of the unnormalized log-posterior distribution at w MAP [START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF].

There are plenty of other approaches to approximate Bayesian linear models, like Expectation Propagation [START_REF] Opper | Gaussian processes for classification: Mean-field algorithms[END_REF], Evidence Lower Bound maximization, and Markov Chain Monte Carlo sampling, some of which we will consider later in this chapter.

In Chapter 4, we will consider another approach, based on a Dirichlet approximation, that allows replacing the classification problem with the regression one. The Laplace approximation is used in some experiments carried out in this work.

Bayesian Neural Networks

We can extend the Bayesian approach for linear models to non-linear ones, like neural networks. Let's consider a Multilayer Perceptron (MLP). We consider f to be defined as an MLP with L layers, each of which is given by the following formula:

f l (x) = 1 √ D l-1 W l φ( f l-1 (x)) + b l , l ∈ {1, ..., L + 1}, (2.44)
where W l and b l are the weights and the biases of the l-th layer, φ is some nonlinear function (i.e. ReLU), and D l denotes the dimension of the input for the corresponding layer. W l and b l are the stochastic variables. By dividing each layer by √ D l , we ensure that the variance of the output does not explode in the limit where the width of the network tends to infinity; this specification is known as Neural Tangent Kernel parameterization [START_REF] Jacot | Neural Tangent Kernel: Convergence and Generalization in Neural Networks[END_REF]. We shall use w to refer to the set of parameters: w := vec({W l , b l } L l=1 ) ∈ R m to be the corresponding vectorized form. The most popular choice of prior distribution p(w) for BNNs is the fully factorized Gaussian. We will consider the subject of prior for BNNs more thoroughly in Chapter 6.

In this work, we will use mostly Bayesian MLPs, but we can note that the extension of the Bayesian framework to other architectures is straightforward.

Inference techniques for Bayesian Neural Networks

Performing inference for a BNN poses the same challenge as Bayesian Logistic Regression. We cannot obtain the posterior distribution over the weights p(w|D) in closed form because the prior and the likelihood are not conjugate. The two most popular directions to resolve this issue are Variational Inference and Markov Chain Monte Carlo sampling.

Variational inference

The main idea of the Variational Inference approach is to choose a family of distributions and find such an element q(w) in this family that is closest to the true posterior distribution p(w|D) according to some metric. The most popular choice for this metric is the Kullback-Leibler divergence.

KL(p(x)||q(x)) = p(x) log

p(x) q(x) dx.

(2.45)

The classical approach to solving min q KL(q||p) is to maximize the Evidence Lower Bound (ELBO).

log p(w|D) ≥ L(p(w)|D)),

where L(p) = E w∼q(w) log p(w|D) q(w) .

(2.46)

It is possible to formulate an ELBO maximization objective as follows [START_REF] Kingma | Auto-Encoding Variational Bayes[END_REF]: min q -E w∼q(w) [log p(D|w] + KL(q(w)||p(w)).

(2.47)

The classical approach proposes to choose a parametric family of approximate distributions and solve this optimization problem with respect to the parameters of the approximate distribution θ. The problem of this approach is how to choose a rich family for the approximate distribution that can approximate the true posterior well. From the other point of view, the second term in the RHS of (2.47) is analytically available for a small set of priors and approximate posteriors, which makes this approach quite restrictive.

Recently so-called Particle Optimization Variational Inference (POVI) [START_REF] Liu | Stein Variational Gradient Descent: A General Purpose Bayesian Inference Algorithm[END_REF] was proposed as an alternative approach which does not rely on parametric approximations but instead chooses an empirical distribution of particles, where each particle represents a particular realization of the random variable of interest. Good examples of this approach were investigated in [START_REF] Liu | Stein Variational Gradient Descent: A General Purpose Bayesian Inference Algorithm[END_REF], [START_REF] Francesco | Repulsive Deep Ensembles are Bayesian[END_REF] and [START_REF] Milios | A Variational View on Bootstrap Ensembles as Bayesian Inference[END_REF]. We will discuss in detail this particular type of VI in Chapters 4 and 6.

Markov Chain Monte Carlo

MCMC represents a completely different approach to characterize the posterior over model parameters in the Bayesian treatment of statistical models. This approach aims to find a sequence of samples generated by a Markov chain. A Markov chain is a sequence of random variables where the distribution of each variable depends only on the previous element of the sequence. In the case of Bayesian inference, the Markov chain is constructed in such a way as to generate samples from the posterior distribution.

To explain the details of this approach, we consider a Metropolis-Hastings algorithm, which is a very basic MCMC algorithm with slow convergence, and we consider it only for illustration purposes. Our goal is to generate samples from the posterior distribution p(w|D). For this algorithm, it is necessary to introduce an auxiliary distribution q(w ′ |w) that depends on the previous sample in the sequence w, and it is referred to as the is called a proposal distribution. Starting from an initial point w 0 , we generate a new candidate w ′ ∼ q(w ′ |w) at each step. We accept w ′ and add it to the sequence with probability a a = min 1, p(w ′ |D)q(w ′ |w) p(w|D)q(w|w ′ ) .

(2.48)

If the candidate is not accepted, we add to the sequence old sample w to the sequence. It is easy to notice that for this algorithm, we do not need to compute the normalization constant in the denominator of (2.24) because it cancels out during the division operation in (2.48). Thus, the computation of the acceptance probability boils down to a = min 1, p(w ′ |D)q(w ′ |w) p(w|D)q(w|w ′ ) , (2.49) where p(w|D) is unnormalized posterior probability density p(w|D) = p(D|w)p(w).

(2.50)

As mentioned before, the Metropolis-Hastings algorithm has poor convergence properties, and it requires computing the likelihood function p(D|w) over the whole dataset, which could be expensive in the case of large datasets. In this work, we widely use Stochastic Gradient Hamiltonian Monte Carlo (SGHMC) because it has better convergence rates [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF] and it allows for mini-batching, i.e., Monte Carlo estimates of the log-likelihood using mini-batches of data.

The SGHMC algorithm is based on Hamiltonian Monte Carlo (HMC) [START_REF] Radford | MCMC using Hamiltonian dynamics[END_REF]. HMC considers the sampling procedure as a mechanical system with total energy, defined by the Hamiltonian, and position variables w and momentum variables r:

H(w, r) = U(w) + K(r), (2.51)
In this expression, the potential energy is given by

U(w) = -log p(w|D), (2.52)
and the kinetic energy is

K(r) = 1 2 r ⊤ M -1 r. (2.53)
In most cases, the mass matrix M is set to the identity matrix I. To generate a new sample for w, it is necessary to sample r ∼ N (0, M), and then to simulate Hamiltonian dynamics of this system for some time t.

   dw = M -1 rdt dr = -∇U(w)dt (2.54)
To simulate this dynamics, it is necessary to use discrete integrators, which introduce discretization errors. That is why it is necessary to use Metropolis-Hastings correction at the end of the Hamiltonian simulation, i.e., accept the resulting w ′ with probability a given by (2.49).

The HMC algorithm has much better convergence properties but it still requires computing the exact likelihood p(D|w), giving it a limited scalability. The SGHMC algorithm allows resolving this issue using the augmented Hamiltonian dynamics:

   dw = M -1 rdt dr = -∇U(w)dt -BM -1 rdt + N (0, 2Bdt), (2.55)
where B = 1 2 ϵV, with the discretization step ϵ of the numerical integrator and the covariance V, that comes from the mini-batched estimation of the gradient. The intuition behind this dynamics is the following. Utilization of the gradient estimates via mini-batches introduces some noise that, due to the central limit theorem, is assumed to have a Gaussian distribution. Because the dynamics have to preserve entropy, it is necessary to introduce a friction term -BM -1 r to the momentum dynamics.

A discretized version of this dynamics, implemented in practice, is the following:

w i = w i-1 -ϵ t M -1 r i-1 r i = r i-1 -ϵ t ∇ Ũ(w i ) -CM -1 r i-1 + N (0, 2(C -B)ϵ t ), (2.56)
where C is a user-defined friction parameter, Ũ(w i ) is a mini-batched estimate of the U((w i )). The empirical noise covariance B can be estimated using a method proposed in [START_REF] Tobias Springenberg | Bayesian Optimization with Robust Bayesian Neural Networks[END_REF].

I N P U T D ATA B I N A R I Z AT I O N F O R O P T I C A L P R O C E S S I N G U N I T S

introduction

Statistical models based on kernel methods offer powerful and theoretically wellunderstood tools for complex data modeling problems. The limitation of employing these kernel-based models in practice is that a naive implementation scales poorly with the size of the data set, and there has been a tremendous amount of work in the direction of mitigating this issue by introducing approximations.

In this context, Nyström approximations [START_REF] Ronald | Simple statistical gradient-following algorithms for connectionist reinforcement learning[END_REF] and random features [START_REF] Rahimi | Weighted Sums of Random Kitchen Sinks: Replacing minimization with randomization in learning[END_REF] are very popular techniques to scale kernel methods virtually to any number of data, thanks to mini-batch formulations [START_REF] Cutajar | Random Feature Expansions for Deep Gaussian Processes[END_REF], [START_REF] Hensman | Gaussian Processes for Big Data[END_REF].

The focus of this chapter is on random feature approximations, whereby kernelbased models are "linearized" by an equivalent linear model with a set of suitably constructed random basis functions. The motivation behind this work is to considerably accelerate the construction of random features while reducing power consumption by resorting to dedicated hardware, which we refer to as Optical Processing Unit (OPU).

OPUs are computing devices that perform random projections of input vectors by exploiting the physical phenomenon of scattering a light source through a diffusive medium [START_REF] Saade | Random projections through multiple optical scattering: Approximating Kernels at the speed of light[END_REF]. The random projection is then followed by a nonlinear operation, making the whole pipeline of computation exactly what is needed to construct random features to approximate kernel-based models. Crucially, OPUs offer the possibility to operate with a number of random features at the speed of light and with low power consumption, representing a unique solution to improve the scalability of kernel machines further. As an example, OPU-based random feature approximations have successfully been proposed to carry out approximate kernel ridge regression in [START_REF] Ohana | Kernel Computations from Large-Scale Random Features Obtained by Optical Processing Units[END_REF], [START_REF] Cappelli | Adversarial Robustness by Design Through Analog Computing And Synthetic Gradients[END_REF].

One limitation associated with working with OPUs is that, because of the hardware setup, input vectors need to be binarized. In addition, the random projection matrix characterizing the device is unknown and can only be retrieved through an expensive calibration procedure. In this paper, we propose a novel binarization strategy for OPUs, which is learned along with the regression/classification task in an end-to-end manner, meaning that the parameters of the binarization part are learned along with the kernel-based model parameters. In order to achieve this, we overcome the limitation that OPU projection matrices are unknown by employing the so-called REINFORCE gradient estimator, which allows us to treat the OPU as a black-box. Through experiments on several UCI classification/regression problems, we show that our proposal outperforms alternative unsupervised and supervised binarization techniques. We provide some analysis of the bounds on the objective functions of the proposed approaches, and also we consider several classes of kernels and image-based classification problems.

3.2 overview of the subject field

Binary Neural Networks

In neural networks, binarization is generally targeting intermediate layer activations, and it may also stem from the binarization of model parameters. In these cases, binarization is mostly introduced to reduce computational cost and memory consumption [START_REF] Haotong Qin | Binary neural networks: A survey[END_REF]. Neural networks with binary hidden layers find applications in binary autoencoders for hashing [START_REF] Miguel | Hashing with binary autoencoders[END_REF], data compression [START_REF] Tissier | Near-Lossless Binarization of Word Embeddings[END_REF], and hard attention mechanism [START_REF] Xu | Show, Attend and Tell: Neural Image Caption Generation with Visual Attention[END_REF]. The binarization of layer activations is obtained by a suitable choice of activation functions; for instance, the sign or Heaviside functions for the deterministic case, or the sigmoid or tanh functions combined with the Bernoulli distribution for the stochastic case [START_REF] Courbariaux | Binarized Neural Networks: Training Deep Neural Networks with Weights and Activations Constrained to +1 or -1[END_REF], [START_REF] Jorn | Probabilistic Binary Neural Networks[END_REF]. The most popular technique to propagate gradients through such activation functions is the so-called straight-through estimator (STE) [START_REF] Bengio | Estimating or Propagating Gradients Through Stochastic Neurons for Conditional Computation[END_REF]. This technique assumes that the gradient of non-differentiable activation functions is equal to 1. Of course, this is quite a crude assumption and it leads to biased estimations of gradients of estimated variables, but it works surprisingly well in practice. More recently, there have been proposals to replace the STE with another estimator through a relaxation technique, also known as the Gumbel Softmax-trick [START_REF] Jang | Categorical Reparameterization with Gumbel-Softmax[END_REF].

Also, different kinds of target propagation are used to learn suitable targets for each binary layer and then train the associated parameters with relaxation techniques, or combinatorial optimization [START_REF] Friesen | Deep Learning as a Mixed Convex-Combinatorial Optimization Problem[END_REF], [START_REF] Lee | Difference Target Propagation[END_REF], [START_REF] Carreira-Perpinan | Distributed optimization of deeply nested systems[END_REF].

Auxiliary models for binarization

Focusing on OPUs, currently, the standard approach to binarize data makes use of a binary autoencoder [START_REF] Lighton | Lightonml documentation[END_REF], [START_REF] Tissier | Near-Lossless Binarization of Word Embeddings[END_REF]. Such a binary autoencoder is trained independently from the OPU device, and it gives the possibility to perform the binarization operation by means of its encoder part. The autoencoder consists of a fullyconnected encoder and decoder. The hidden layer has a Heaviside activation function, so its output is binary. The training procedure updates the weights of the decoder with backpropagation, and the weights of the encoder are forced to be equal to the weights of the decoder in order to be able to reproduce the input.

This task could be performed with a more advanced auxiliary model for binarization. The current state-of-the-art for autoencoders with latent binary space was proposed in [START_REF] Fajtl | Latent Bernoulli Autoencoder[END_REF]. While this paper is mostly focused on sampling in the latent Bernoulli space and generating new data, the proposed model itself provides very competitive image reconstruction performance even compared with non-binary autoencoders. At its core, it still uses the STE approach for training, but it introduces a novel method of sampling latent representations at the test phase. Thus, it could be beneficial to use this model to input binarization for OPU-based models.

Optimization of heterogeneous models

In this work, we aim to develop a supervised binarization model which is learned together with the supervised learning task. That is, we aim to provide a training procedure for the heterogeneous model, consisting of the kernel ridge regression model approximated with random features and the binarization encoder before the OPU.

In the context of heterogeneous models, a general-purpose framework called Method of Auxiliary Coordinates (MAC) was proposed in [START_REF] Carreira-Perpinan | Distributed optimization of deeply nested systems[END_REF] with examples of application in [START_REF] Miguel | Hashing with binary autoencoders[END_REF] and [START_REF] Choromanska | Beyond Backprop: Online Alternating Minimization with Auxiliary Variables[END_REF]. The authors propose to introduce auxiliary variables into a deep neural network. These auxiliary variables are assigned the role of pre-activations for each layer, and they get replaced during the forward pass. The first step of the optimization targets the auxiliary variables, and, after this step, the parameters of each layer are optimized to regress on these variables, which take the role of layer-specific labels. This is very beneficial when some layers are discrete and vanilla backpropagation is not applicable. In [START_REF] Choromanska | Beyond Backprop: Online Alternating Minimization with Auxiliary Variables[END_REF], this approach is used to train a fully connected network with binary activation functions, using an STE to propagate a learning signal through the non-differentiable parts. Reference [START_REF] Miguel | Hashing with binary autoencoders[END_REF] is especially interesting because the authors illustrate how discrete binary layers can be optimized within a larger, non-binary model.

While splitting the optimization of the binarization and the model is a viable option, we still need a way to train each part individually. There is a wide variety of ways to obtain a solution for kernel ridge regression with the random feature approximation, so the most difficult point is how to optimize the part consisting of the binary encoder and the OPU, because it combines a non-differentiable binarization function with an implicit random projection, which is a black-box function, so it is also non-differentiable at each point. These make the STE from [START_REF] Choromanska | Beyond Backprop: Online Alternating Minimization with Auxiliary Variables[END_REF] inapplicable. Also, we found that the combinatorial approach used in [START_REF] Miguel | Hashing with binary autoencoders[END_REF] and [START_REF] Friesen | Deep Learning as a Mixed Convex-Combinatorial Optimization Problem[END_REF] is inapplicable to our case for two reasons. First, it is suitable only when the binary dimension is relatively small, which might be a limitation for a general solution. Second, the combinatorial approach combined with MAC converges in one iteration to poor local optima, and this happens because of the model setup, which is different from the ones in [START_REF] Miguel | Hashing with binary autoencoders[END_REF] and [START_REF] Friesen | Deep Learning as a Mixed Convex-Combinatorial Optimization Problem[END_REF].

Reinforcement learning approach

From a different point of view, it is possible to view our problem through the lenses of reinforcement learning, where it is necessary to propagate binary codes through the OPU instead of discrete actions through the black-box environment. Instead of maximizing the reward from the environment, we are trying to minimize the loss function. The classical algorithm to solve this problem is REINFORCE [START_REF] Ronald | Simple statistical gradient-following algorithms for connectionist reinforcement learning[END_REF]. This allows one to calculate the gradients of the reward with respect to the parameters of the policy that generates actions. The applicability of this method to other settings with black-box elements was shown in [START_REF] Ranganath | Black Box Variational Inference[END_REF]. There are various versions of this algorithm intended to reduce the variance of the gradient of the parameters. Very frequently, they are based on relaxations of the non-differentiable sampling procedure [START_REF] Tucker | REBAR: Low-variance, unbiased gradient estimates for discrete latent variable models[END_REF], or approximation of the black-box part of the model [START_REF] Grathwohl | Backpropagation through the Void: Optimizing control variates for black-box gradient estimation[END_REF]. It is also worth noting that there exist competitive alternatives to REINFORCE, such as the Augment-REINFORCE-Merge (ARM) [START_REF] Yin | ARM: Augment-REINFORCE-Merge Gradient for Stochastic Binary Networks[END_REF], later extended with variance reduction [START_REF] Kool | Buy 4 REINFORCE Samples, Get a Baseline for Free![END_REF] or relaxation [START_REF] Dong | DisARM: An Antithetic Gradient Estimator for Binary Latent Variables[END_REF]. However, the experimental results, reported in the publications dedicated to different modifications of ARM, show that the performance gain is marginal compared to a REINFORCE method with a proper variance reduction technique. And given the fact that the REINFORCE method has proven itself as a very reliable technique applied to a large variety of tasks and datasets, we prefer to adopt this in the optimization of OPU-based models.

Variance reduction for REINFORCE

The most popular approaches that are used for variance reduction of REINFORCElike estimators are conditioning [START_REF] Owen | Monte Carlo theory, methods and examples[END_REF] and control variates [START_REF] Paisley | Variational Bayesian Inference with Stochastic Search[END_REF].

The conditioning technique is based on the Rao-Blackwell theorem. The tech-That is why, in the literature, it is often referred to as Rao-Blackwellization nique is intended to reduce the variance of the estimator of the expectation

E[ f (x, y)] of the function f over the joint distribution x, y ∼ p(x, y). If it is possi- ble to obtain the conditional expectation E[ f (x, y)|y], it is possible to define an estimator f f (y) = E[ f (x, y)|y]. (3.1)
First, we can note that

E[ f (y)] = E[ f (x, y)]. (3.2)
It means that it is possible to replace the standard empirical estimator of the expectation

E[ f (x, y)] ≈ 1 n n ∑ i f (x, y), x, y ∼ p(x, y) (3.3)
with the conditional estimator

E[ f (x, y)] ≈ 1 n n ∑ i f (y), y ∼ p(y). (3.4)
We can easily compute the variance of the conditional estimator as:

var( f (y)) = var( f (x, y)) -E[( f (x, y) -f (x)) 2 ]. (3.5)
Thus, the main motivation to do this replacement is the fact that the variance of the one-sample conditional estimator is always smaller than the variance of the original one.

var 1 n n ∑ i=1 f (y) ≤ var 1 n n ∑ i=1 f (x, y) , x, y ∼ p(x, y). (3.6)
The main problem of this approach is that the straightforward analytical expression for E[ f (x, y)|y] sometimes is not available for the chosen model. The control variates technique provides another way to obtain an estimator with reduced variance. Unlike the previous method, this technique reduces the variance of the estimator of E[ f (x)], x ∼ p(x). So, it is applicable to functions that depend on one variable, and we don't have to compute any conditional expectation. For this approach, it is necessary to find an auxiliary function g that is somehow similar to the target function f . The nature of this similarity will be revealed later. Let's consider the following estimator:

f (x) = f (x) -a(g(x) -E[g(x)]). (3.7)
There are different types of estimators that use control variates. The estimator (3.7) is one of the most popular ones because it allows for the optimization of variance reduction with respect to the parameter a. The main requirement for this type of control variates is the availability of the analytical form of

E[g(x)]. Again, we can note that E[ f (x)] = E[ f (x)].
And the most important thing is that this estimator has a lower variance compared to the original one:

var( f (x)) = var( f (x)) + a 2 var(g(x)) -2acov( f (x), g(x)). (3.8)
From this expression, we can see that by the similarity measure between f and g is the implied covariance between these functions. Also, we can find the minimum of the expression (3.8) with respect to the parameter a and find the optimal value a * :

a * = cov( f (x), g(x)) var(g(x)) , x ∼ p(x). (3.9)
In practice, it is too expensive to compute cov( f (x), g(x)) and var(g(x)) over the whole distribution p(x), so we have to use the Monte Carlo approximation.

random features on optical processing units

In this section, we discuss OPUs in the context of random features. In the previous section, we discussed random features as a way to approximate models involving kernels; for OPUs, instead, the device produces random features (fast and with little power consumption) and the question that we aim to address here is how to use these to implement approximate kernel machines.

OPUs are computing devices that exploit the physical process of scattering of light to perform a random projection operation of a given vector. In particular, given a binary vector x i ∈ R d , OPUs perform multiplication by a random matrix R and apply the nonlinear activation function ∥ • ∥ 2 . In other words,

ϕ(x) = 1 √ D ∥Rx∥ 2 (3.10)
The matrix R ∈ C D×d is a complex Gaussian matrix with elements R ij ∼ CN (0, 1). Previous works have established that in the limit of an infinite number of random features, the equivalent kernel is the following [START_REF] Ohana | Kernel Computations from Large-Scale Random Features Obtained by Optical Processing Units[END_REF]:

k(x, y) ≈ ϕ(x)ϕ(y) D→∞ = ∥x∥ 2 ∥y∥ 2 + (x ⊤ y) 2 (3.11)
Therefore, when using OPUs for kernel ridge regression, we are implicitly working with this polynomial kernel.

Recently, a new version of OPUs has been proposed and developed in [START_REF] Dua | UCI Machine Learning Repository[END_REF], which allows one to perform linear random feature projections. To avoid confusion with the polynomial OPU random features, we will denote them ψ(x).

ψ(x) = CRx (3.12)
where C is a fixed constant. This novel type of OPU opens the possibility to approximate a wide variety of kernels by choosing an appropriate activation function [START_REF] Rahimi | Weighted Sums of Random Kitchen Sinks: Replacing minimization with randomization in learning[END_REF], [START_REF] Wacker | Improved Random Features for Dot Product Kernels[END_REF]. For example, it is possible to apply trigonometric activation functions to the outputs of the OPU:

ψ ′ (x) = sin(ψ(x)) cos(ψ(x)) (3.13)
This type of random features is called Random Fourier Features (RFF). It was proven in [START_REF] Rahimi | Weighted Sums of Random Kitchen Sinks: Replacing minimization with randomization in learning[END_REF] that this kind of random features allows approximating RBF kernels.

1 D D ∑ i=1 ψ ′ (x) ⊤ ψ ′ (y) = 1 D D ∑ r=1 sin(ψ(x)) cos(ψ(x)) ⊤ sin(ψ(y)) cos(ψ(y)) = E ω [cos(ω(x -y))] = k RBF (x, y) (3.14)
As mentioned before, an important aspect of OPUs is that their input should be binary; this project proposes a novel way to carry out a binarization of its input along with the kernel ridge regression task in an end-to-end fashion.

combinatorial optimization

In order to be able to implement kernel ridge regression on OPUs we need to binarize the inputs x i . We propose to do so by employing an encoder implemented as a neural network and parameterized by a set of weights W enc . A combination of the OPU-based regression with the binarization network is a heterogeneous model.

Let's look at the application of the MAC on the OPU-based regression model. In this case, we have to propagate a temporary learning objective for the binary encoder through the optical random projection layer. Let's denote the desired output of the binary encoder as s and auxiliary variables corresponding to the activations of the OPU layer ϕ(•) as c.

In Fig. 3.1, we plotted the beginning of the optimization. Here, we use a 2D plane only for illustration purposes. In general, the discussed vectors are not twodimensional. Vectors in the space of auxiliary variables are depicted in blue, and binary vectors that correspond to the output of the OPU are in red. Originally, we have an initial binary output of the encoder s (0) . We initialize the auxiliary variable as c (0) = ϕ(s (0) ). Auxiliary vector c (0) is depicted collinear to the vector s (0) only for demonstration purposes to emphasize the deterministic dependence of c (0) from s (0) . And further, we will follow the same scheme of showing deterministic dependence between vectors. In general, these vectors lay in different spaces with different dimensionality. At the first step, we have to solve regr are the initial weights of the regression model and µ is a regularization parameter. Let's note, that c (0+ 1 2 ) belongs to a span of vectors c (0) (contribution of regularization) and c * (0) (contribution of data), where

c (0+ 1 2 ) = argmin c [∥y -cW (0) regr ∥ 2 + µ∥c -ϕ(s (0) )∥ 2 ]. (3.15)
c * (0) = argmin c [∥y -cW (0) regr ∥ 2 ] (3.16)
. Then, we update the weights of the binary encoder.

W

(1)

regr = argmin W regr [∥y -c (0+ 1 2 ) W regr ∥ 2 ] (3.17)
At the end of the iteration, we perform a combinatorial search for the optimal s (1) among all possible candidates s [START_REF] Gogolashvili | Locally Smoothed Gaussian Process Regression[END_REF] candidate , which generates the closest to c (0+ 1 2 ) activations ϕ(s) (Fig. 3.2). It is possible to argue that an exhaustive combinatorial search is not a very practical solution, but for illustrative purposes, in this example, we assume that we have access to some discrete optimization procedure that gives us a solution for the following optimization objective:

s (1) = argmin s ∥c (0+ 1 2 ) -ϕ(s)∥ 2 . (3.18)
On the next iteration, the new vector c (1+ 1 2 ) will be a solution of (3.15), and it will lay within a span of the vector c ( * 1) = c (0+ 1 2 ) of the data contribution and the vector c (1) = ϕ(s (1) ) of the regularization contribution.

It means that the closest vector of activations generated by vectors in discrete space of binary encodings will be again c (1) . The same will happen for all next iterations (Fig. 3.3). This optimum is mostly defined by W (0) , s (0) , µ

This argumentation shows that the MAC approach is not applicable for the models that involve OPUs because it almost immediately gets stuck in a poor local optimum. That is why we have to switch to another approach to train the binary encoder for input data.

reinforce for kernel regression with binarized inputs

Another approach to implement kernel ridge regression on OPUs is using stochastic binarization of the inputs x i , and we propose to do so by employing an encoder, implemented as a neural network with a stochastic output, parameterized by a set of weights W enc . The encoder transforms the inputs to kernel-based models x i and turns them into a set of Bernoulli-distributed binary random variables z i .

In particular, we denote by f k (x, W enc ) the function implemented by the encoder which parameterizes the Bernoulli distribution associated with the kth element of the output, that is z k .

Recalling the random feature formulation of linear regression of Chapter 2, we propose the following approach to construct an approximate kernel-based model with binary inputs:

y = E z [w ⊤ regr ϕ(z)] + ϵ, ( 3.19)
where z ∼ Bernoulli( f (x, W enc )) and w regr are the parameters of the linearized regression model. Note how in this formulation, the binary vectors z are treated stochastically due to the expectation under the Bernoulli distribution induced by the encoder. The reason for this is that it allows us to employ the so-called REINFORCE gradient estimator, as we discuss next. REINFORCE, also known as the log-derivative trick or score function estimator, offers a way to estimate the gradient of the expectation of a non-differentiable function f (z) under the distribution of the input random vector variables z:

∇ θ E p(z;θ) f (z) = ∇ θ p(z; θ) f (z)dz = (3.20) ∇ θ p(z; θ) f (z)dz = p(z; θ) ∇ θ p(z; θ) p(z; θ) f (z)dz = (3.21) E p(z;θ) ∇ θ log p(z; θ) f (z) ≈ 1 M M ∑ i=1 ∇ θ log p(z; θ) f (z) (3.22)
where M is number of samples drown from p(z, θ). Applying REINFORCE to our approximate kernel-based model yields the following optimization objective:

min w regr ,W enc E Z∼Bernoulli( f (X,W enc )) [L(y, ϕ(Z)w regr )]+ λ enc ∥W enc ∥ 2 + λ regr ∥w regr ∥ 2 (3.23)
In this expression, we denoted by L(y, ỹ) the loss function associated with the task at hand and by Z the matrix that contains binary encoded variables for the whole training set X. We can optimize this objective by means of gradient-based techniques; for this, we require that we are able to compute the gradient of the objective with respect to all parameters. The gradient of the first term of the objective with respect to W enc , which is the most involved part, is:

∇ W enc E z∼q(z) [L(y, w ⊤ regr ϕ(z))] ≈ ≈ 1 M M ∑ i=1 L(y, w ⊤ regr ϕ(z i ))∇ W enc log q(z i ) (3.24)
while the derivatives of the other terms are straightforward to compute. With this derivation, we observe that it is then possible to jointly optimize all parameters, leading to what it is commonly referred to as an end-to-end approach. In the remainder of this paper, we refer to this method as End-to-End SE, where SE stands for Supervised Encoder.

variance reduction

A straightforward application of the expressison (3.24) for gradient optimization will lead to a slow convergence due to the variance of the REINFORCE gradient estimator. In this section, we will discuss the application of variance reduction techniques for the proposed OPU pipeline.

For the proposed approach, there is no straightforward way of applying the Rao-Blackwellization technique because, in our case, the stochastic activations of the encoder are mutually independent. It means that the conditional expectation of each dimension of the activation vector z with respect to all other dimensions is equal to its marginal distribution, and the conditional estimator becomes useless. However, the Rao-Blackwellization technique can still be applied, but it will require changing the original pipeline or replacing the original Bernoulli stochastic activations with auxiliary variables [START_REF] Yin | ARM: Augment-REINFORCE-Merge Gradient for Stochastic Binary Networks[END_REF]. It means that for this particular setting, we have to use control variates.

There are multiple ways to apply of the control variates technique.As a first attempt, we could apply the regression estimator (3.7), which is one of the most popular approaches. In this scenario, we have to find a control variate for h

(z) = ∇ W enc E z∼q(z) [L(y, w ⊤ regr ϕ(z))].
The obvious choice for the control variate is the score function:

g(z) = ∇ W enc log q(z) (3.25)
This choice of control variate has an analytically available form for the expectation E[h(z)] = 0. And it has a computational benefit because the score function ∇ W enc log q(z) should be computed anyway for the gradient estimator itself (3.24). Thus, the modified gradient estimator has the following form:

∇ W enc E z∼q(z) [L(y, w ⊤ regr ϕ(z))] ≈ ≈ 1 M M ∑ i=1 (L(y, w ⊤ regr ϕ(z i )) -a)∇ W enc log q(z i ). (3.26)
Then we have to compute an optimal value for the parameter a according to (3.9), which implies computation of the empirical variance var(h(z)) and the empirical covariance cov(h(z), g(z)). The key element, which is necessary for this computation, is the evaluation of the ∇ W enc log q(z i ) with respect to each sample z i ∼ q(z).

All these evaluations should be stored and then used for the computation of the variance and covariance. This observation shows that this gradient estimator is computationally inefficient for the proposed pipeline because we have to evaluate separately the gradients for each data point and for each binary representation of this data point sampled from the encoder. In order to reduce the variance of this estimator, we employ control variates [START_REF] Kool | Buy 4 REINFORCE Samples, Get a Baseline for Free![END_REF]. In this approach, we add a set of random variables to the estimator, such that these variables have zero means, so they do not alter the expectation of the gradient. The aim is to construct such variables so as to reduce the overall variance of the estimator:

∇ W enc E z∼q(z) [L(y, w ⊤ regr ϕ(z))] ≈ 1 M M ∑ i=1 ∇ W enc log q(z i )(L(y, w ⊤ regr ϕ(z i )) -v i )
where

v i = 1 M -1 ∑ i̸ =j L(y, w ⊤ regr ϕ(z j )) (3.27)

lowering the cost of reinforce

The estimation of the gradient of the End-to-End SE with respect to W enc can be expensive when the number of random features is large. This is due to the fact that this requires multiple samples to be passed from the encoder through the random projection and the approximate kernel ridge regression model. In this section, we propose a strategy to reduce the complexity of REINFORCE applied to our model, whereby we average the set of basis functions under the resampling of the binary variables as follows:

ỹ = w ⊤ regr E[ϕ(z)] + ϵ (3.28)
where z ∼ Bernoulli( f (x, W enc )).

With this new modeling assumption, the training is based on a modified optimization problem as follows:

min w regr ,W enc L(y, E Z∼Bernoulli( f (X,W enc )) [ϕ(Z)]w regr )+ λ enc ∥W enc ∥ 2 + λ regr ∥w regr ∥ 2 (3.29)
Again, we can perform gradient-based optimization. Focusing on the first term, which is the nontrivial one to differentiate in the objective, we obtain

∇ W enc L = dL d(E[ϕ(z]) ∇ W enc E[ϕ(z)] + 2λ enc ∑ ij W enc ij (3.30)
where ∇ W enc E[ϕ(z)] is calculated with the REINFORCE estimator. In the remainder of the paper, we will refer to this method as Isolated Supervised Encoder (SE).

Regarding the comparison of the End-to-End SE and Isolated SE, we can note the following relationship between these models in the case of regression problems. In the data term of the optimization objective (3.29) we can put an expectation over the whole matrix product of the random features map and the regression weights instead of an expectation over the random features only. Then we can move the expectation in such a way that it is taken over the whole term within the squared norm. We can do this because, within one gradient step iteration, neither y nor W enc are considered as random variables. In this case, the data term looks as follows: In turn, End-to-End SE has the following data term as part of its optimization objective (3.23):

∥y -E[ϕ(Z)]w regr ∥ 2 = ∥y -E[ϕ(Z)w regr ]∥ 2 = ∥E[y -ϕ(Z)w regr ]∥ 2 (3.31)
E[∥y -ϕ(Z)w regr ∥ 2 ] (3.32)
We can note that the squared loss is a convex function. Thus, we can apply Jensen's inequality to obtain the following expression:

E[∥y -ϕ(Z)w regr ∥ 2 ] ≥ ∥E[y -[ϕ(Z)]w regr ]∥ 2 = ∥y -E[ϕ(Z)]w regr ∥ 2 (3.33)
As a result, End-to-End SE optimizes the upper bound of the Isolated SE objective.

experiments

Experiments on the UCI datasets

We compared the performance of the proposed approaches for a non-linear OPU (End-to-End SE and Isolated SE) and a linear OPU that uses trigonometric activations (End-to-End SE with RFF) against a model based on unsupervised autoencoder proposed in [START_REF] Tissier | Near-Lossless Binarization of Word Embeddings[END_REF], an encoder trained with Direct Feedback Alignment (DFA) [START_REF] Nøkland | Direct Feedback Alignment Provides Learning in Deep Neural Networks[END_REF] and a KRR based on a RBF. Results are reported in Fig. 3.4 for several UCI regression and classification problems [START_REF] Dua | UCI Machine Learning Repository[END_REF]. We want to emphasize that the main competitors of the proposed methods are the ones based on unsupervised autoencoder and encoder trained by DFA, because kernel ridge regression is unable to work with large datasets, and OPU-based regression just approximates this method and is intended to replace it on large datasets.

For the KRR experiments, we used MSE as a loss function. To apply KRR to the classification problems, we replaced 0 and 1 in class labels with -1, 1 and solved a classification problem as a regression one using MSE loss as an optimization objective. For all other models, we used MSE loss for the regression problems and Cross Entropy (CE) loss for the classification problems. We used a logistic activation function on the last layer for the classification tasks.

For Isolated SE and End-to-End SE as an encoding function f (x, W enc ) providing parameters for the Bernoulli distribution, we chose a single linear layer with sigmoid activation.

f (x, W enc ) = σ(W ⊤ enc x) (3.34)
All hyperparameters for the DFA encoder, End-to-End SE, and Isolated SE models (size of binary embedding, learning rate, L2 regularization for the encoder and the regression layer) were chosen with a random search during cross-validation. Kernel parameters of KRR were tuned by random search with cross-validation. This poses computational challenges for the large datasets (MiniBoo, MoCap), so we resort to RFF approximations for these cases.

For the models involving random features (both Fourier and OPU-generated ones), we have tuned the variance of the distribution that generates these random features. Concretely, assuming that the elements of the R matrix generating the random projections are distributed through the standard Normal distribution, we can obtain a new random matrix R ′ by multiplying R by any variance, for instance:

ϕ ′ (x) = c|R ′ x| 2 = c| R α x| 2 = c 1 α 2 |Rx| 2 (3.35)
It is enough to multiply the output of the OPU by an additional parameter γ, such that γ 2 = 1 α 2 , and α optimize them with standard gradient descent. The parameter γ is not equivalent to the lengthscale parameter of the RBF kernel. In practice, it has an effect on the outputscale parameter of the RBF kernel, as it has simply a scaling effect on the kernel.

On the regression problems, both proposed methods outperformed their main competitors. On the classification problems, the DFA-based approach was better only on one dataset, and on all other datasets, the proposed methods performed better or equally well. Regarding the type of kernel approximated by the OPU, the experiments show that the linear OPU with trigonometric activations performs as good as OPU kernel for most of the datasets. It gives performance gains only for some classification problems. Considering the comparison between the proposed methods, we see that End-to-End SE is more stable and requires a significantly fewer number of samples from the encoder, although Isolated SE showed slightly better results on classification problems.

We considered including results obtained by running these models on the real OPU (Fig. 3.5). Unfortunately, the regression problems required such a large number of epochs that we could not perform the experiments in a reasonable amount of time.

We also tested the performance of our approach with respect to the number of samples required to employ REINFORCE. We found that End-to-End SE can achieve good results with a small number of samples from the encoder, and the increase of amount of samples does not seem to improve performance. Finally, we evaluated the effect of variance reduction on convergence speed and performance for the End-to-End SE model. In Fig. 3.6, we report results for one classification and one regression problem. The convergence curves indicate that the convergence speed benefits from the gradient variance reduction.

Experiments on image data

In this section, we evaluate an optical random feature regression approach for image classification tasks with several different binarization techniques, including the proposed methods.

The kernel generated by OPU (3.11) is an example of a polynomial kernel. Polynomial kernels, unlike more popular RBF kernels, take into account interaction between different feature dimensions. This property is especially important for image data because a relative alignment of pixels is crucial for image classification. Of course, when we are working with OPUs, the kernel takes into account an alignment of different dimensions of the binary embedding of an image instead of the pixels. The relative alignment of different dimensions of the binary embedding most probably does not contain exactly the same information as the mutual alignment of pixels. But until the binary encoder does not have disentanglement properties, the mutual interaction of dimensions of the binary embedding has to contain additional information about the image. That is why it is still important to use a kernel that is capable to take into account these relationships.

For the experiments on image data, we used two convolutional architectures of the binary encoder. The first architecture was inspired by the LeNet model (Table 3.1). We performed experiments on three classical image classification datasets (Table 3.2). We compared End-to-End SE with a model that used an autoencoder to train the encoder (AE) and a model that used direct feedback alignment (DFA) for the same purpose. The results are shown in Table 3.3. We used the same encoder architecture with the same hyperparameters for each binarization method. The size of the binary embedding was set to 400, except of the unsupervised AE method.

The reason we trained the unsupervised AE differently for these experiments is because we were using complex convolutional models, and it was hard to adapt the method proposed in [START_REF] Tissier | Near-Lossless Binarization of Word Embeddings[END_REF]. This binarization approach requires using exactly the same values of weights both for the encoder and for the decoder models. It means that this method requires building a decoder model that is symmetrical to the encoder model. Achieving this property for convolutional neural networks is difficult because, for the decoder, it is difficult to pick equivalent symmetric operations for convolutional and pooling layers in the encoder. Transposed convolutions and interpolation operations that are used in decoders for image data are suitable for training the autoencoder by end-to-end backpropagation. They learn operations that are not symmetric to convolutions and pooling layers of the encoder. The method proposed in [START_REF] Tissier | Near-Lossless Binarization of Word Embeddings[END_REF] assumes that only the decoder is trained, and the encoder copies weights from it, which is impossible to do for asymmetric operations in the decoder and encoder. That is why we trained the autoencoder model in a different way. The encoder of the AE model used an tanh activation function at the output. We used a parameter β that controls the steepness of the tanh function. We slowly increased the value of this parameter from β = 1 to β = 100 in the process of training. At the end of the training, the function tanh with a high value of the parameter β is almost equivalent to a shifted and scaled Heaviside function.

2h(x) -1 ≈ tanh(βx), β → ∞ (3.36)
The results of the DFA approach signify that this type of gradient update is not suitable for convolutional models. This observation is supported by other researchers [5,[START_REF] Launay | Principled Training of Neural Networks with Direct Feedback Alignment[END_REF] Unsupervised AE was able to provide acceptable accuracy for the MNIST dataset, but on CIFAR-10, its performance dropped significantly. A possible explanation is that simple convolutional AE is unable to extract a reasonable binary representation of complex images. The AE model used in the experiments was able to reconstruct simple images from the MNIST dataset. But the reconstruction quality of the same model was much worse for the CIFAR-10 dataset. When the dimensionality of the binary embedding was equal to 400, the AE model was unable to generate any sensible images. Thus, we had to increase the size of the binary embedding to 1024. But the reconstructed images were very blurry even with this modification.

Because of the poor performance of the unsupervised AE baseline, we decided to add another baseline to the comparison. For this experiment, we used a RESNETbased convolutional network as the encoder. LBAE approach proposed in [START_REF] Fajtl | Latent Bernoulli Autoencoder[END_REF] implements an autoencoder with a latent binary space. The training procedure of this method is based on the straight-through gradient estimator. The architecture of the binary encoder is represented in Table 3.4 and Table 3.5. This encoder used leaky relu as an activation function.

Table 3.6 contains the results of the comparison between the LBAE-based and End-to-End SE-based encoders in terms of classification error.

Both binarization approaches used the same resnet-based architecture of the encoder network. We dropped the DFA approach from the comparison because of its poor performance.

The results of this experiment showed interesting properties of the unsupervised approach for training the binary encoder. The LBAE-based encoder with a deeper network performed worse than the simpler autoencoder with tanh annealing in terms of classification error. At the same time, the LBAE approach was better in the As with the UCI data, we evaluated the effect of variance reduction.

As we can see, variance reduction plays a crucial role in image classification. Without this technique, the proposed method is unable to train the model for cifar10.

discussion

Recent advances in alternatives to transistor-based hardware are bringing a new wave of sustainable computing solutions for machine learning [START_REF] Xia | Roadmap on emerging hardware and technology for machine learning[END_REF]. This paper focuses on optical-based computing through OPUs [START_REF] Saade | Random projections through multiple optical scattering: Approximating Kernels at the speed of light[END_REF], which performs random- ized projections of binary input vectors at the speed of light with low energy consumption. In this paper, we considered these randomized computations to implement kernel-based models for regression and classification tasks through random feature approximations. In particular, we proposed a novel strategy to binarize the inputs of the given task to be able to employ OPUs inspired by reinforcement learning. The proposed strategy uses an encoder to map the inputs to a set of binary variables and employs the REINFORCE gradient estimator to estimate its parameters jointly with the parameters of the kernel-based model. We also explored ways to reduce the variance of the gradient estimator and accelerate convergence, which is critical in a number of challenging modeling tasks such as image classification. Through a series of experiments, we showed that our proposal outperforms competitors based on unsupervised binarization and those that do not employ gradient information. We are currently investigating our approach in the context of other kernel-based models, such as Gaussian Processes [START_REF] Edward Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF], and their extension to deep models, such as Deep Kernel Learning [START_REF] Andrew G Wilson | Stochastic Variational Deep Kernel Learning[END_REF] and Deep Gaussian processes [START_REF] Cutajar | Random Feature Expansions for Deep Gaussian Processes[END_REF].

B AY E S I A N A P P R O A C H T O C L A S S I F I C AT I O N

introduction

The Bayesian treatment of statistical models is desirable in applications where quantification of uncertainty is a primary requirement. For many classes of models, this is analytically intractable, and one needs to resort to approximations. Given a statistical model with parameters on which a prior distribution is assumed, such approximations yield an approximation to the posterior distribution over these parameters, either in closed form or in the form of samples. Popular approaches include the Laplace Approximation, Variational Inference, and Markov chain Monte Carlo.

In this work, we focus on an alternative approach called variational bootstrap [START_REF] Milios | A Variational View on Bootstrap Ensembles as Bayesian Inference[END_REF]. Variational bootstrap works by producing a set of replicas of the data set with perturbed labels. Then, each of these perturbed problems is solved by maximuma-posteriori-type optimization. Perhaps not surprisingly, in the Bayesian linear regression case, it is possible to introduce perturbations in a way such that the set of solutions to the perturbed problems is distributed exactly as the posterior over the parameters. One remarkable property of this approach is that it transforms the problem of characterizing the posterior distribution over model parameters into a set of easily parallelizable optimization problems.

Milios et al. [START_REF] Milios | A Variational View on Bootstrap Ensembles as Bayesian Inference[END_REF] provide an extension of this result to the case of regression with deep neural networks featuring ReLU activations, where certain theoretical guarantees are given for the minimization of the KL divergence between the approximate and the true posterior.

The main limitation of this approach is that no guarantees are provided for non-Gaussian likelihoods, such as Bernoulli or Multinomial, which are associated with classification problems.

In this work, we broaden the scope of variational bootstrap by extending it to classification problems. We propose to transform Bernoulli/Multinomial distributed labels to a latent representation with Gaussian noise. We can then apply a model with a Gaussian likelihood to solve a regression problem within this latent space. To carry out this transformation, we follow the method in [START_REF] Milios | Dirichlet-based Gaussian Processes for Largescale Calibrated Classification[END_REF], where classification labels are interpreted as the output of a Dirichlet distribution. We thus propose the combination of variational bootstrap with Dirichlet-based classification; this allows one to obtain reliable uncertainty estimates for the classification model at a lower cost than other Bayesian inference methods, and it enables easy parallelization.

We study the proposed extension to the nonlinear case on deep neural networks for classification with ReLU activations. The transformation of the labels from discrete to continuous allows us to borrow the results obtained in the Gaussian likelihood case. As a result, we obtain a method for which we have guarantees that the optimization of the perturbed problems yields an improvement of the approximation to the posterior over model parameters. In the experiments, we showcase results on various data sets, demonstrating that this is a competitive approach to characterize the posterior over model parameters. Crucially, our proposal is extremely easy to parallelize, and we view this as a considerable advantage compared to alternatives to obtain the posterior over model parameters.

The proposed extension is also relevant for Bayesian logistic regression applied to large datasets with a large number of features. An interesting application of this approach that we showcase in the paper is the approximation of Gaussian process classifiers with random features [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF], which turns the model into a Bayesian logistic regression model. The larger the set of random features, the better the approximation, but this has a negative impact on the computational cost. The proposed approach yields a very practical and effective method to overcome these difficulties.

related work

In recent years, neural networks (NN) have gained popularity due to their effectiveness on a broad variety of tasks, including image classification, [START_REF] He | Deep Residual Learning for Image Recognition[END_REF], natural language processing [START_REF] Devlin | BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding[END_REF], and many others [START_REF] Bourilkov | Machine and Deep Learning Applications in Particle Physics[END_REF][START_REF] Gawehn | Deep learning in drug discovery[END_REF][START_REF] Nguyen | Deep learning methods in transportation domain: a review[END_REF][START_REF] Purwins | Deep Learning for Audio Signal Processing[END_REF]. One of the most persistent challenges is that NNs tend to make overconfident decisions [START_REF] Guo | On Calibration of Modern Neural Networks[END_REF][START_REF] Kendall | What Uncertainties Do We Need in Bayesian Deep Learning for Computer Vision?[END_REF][START_REF] Kurakin | Adversarial examples in the physical world[END_REF]. In the literature, overconfident decisions are treated by introducing uncertainty to the outputs of a NN. In practice, it means that the model produces a predictive distribution of outputs for each input vector. This practice offers ways to quantify the uncertainty of predictions, as high predictive variance implies a lack of confidence. Methods that provide uncertainty quantification for NNs include deep ensembles [START_REF] Balaji Lakshminarayanan | Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles[END_REF], Monte Carlo dropout [START_REF] Gal | Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning[END_REF], and Bayesian Neural Networks (BNN) [START_REF] David | Bayesian neural networks and density networks[END_REF].

In this work, our focus is on BNNs. In contrast with the deterministic networks, BNNs consider parameters to be random variables that are associated with some prior distribution. Instead of a training procedure based on optimization, this treatment requires characterizing a posterior distribution over their weights given observations by means of Bayes theorem. The main challenge is that this procedure is intractable for nonlinear models, which has motivated the development of approximations using, for instance, Variational Inference (VI) [START_REF] Blundell | Weight Uncertainty in Neural Network[END_REF][START_REF] Gal | Bayesian Convolutional Neural Networks with Bernoulli Approximate Variational Inference[END_REF]. These techniques, while being relatively cheap in terms of computational resources, require selecting a family of distributions that is used for approximation of the posterior distribution. The choice of a family of approximate posteriors has a crucial effect on the performance of the model and on the reliability of uncertainty estimations. Because different models and tasks require different families of approximate posterior distributions, there are no reliable heuristics for this task. As an alternative to approximate inference, it is possible to use Markov Chain Monte Carlo (MCMC) techniques, which allows one to get samples from the true posterior distribution of the weights. The recently proposed Stochastic Gradient Hamiltonian Monte-Carlo (SGHMC) method [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF] allows obtaining samples from a true posterior distribution over the weights for large-scale problems using mini-batching. Despite a significant reduction of the time complexity compared to other MCMC methods, SGHMC may still be slow to converge.

Another family of Bayesian inference techniques is particle optimization variational inference (POVI) [START_REF] Francesco | Repulsive Deep Ensembles are Bayesian[END_REF][START_REF] Liu | Stein Variational Gradient Descent: A General Purpose Bayesian Inference Algorithm[END_REF][START_REF] Milios | A Variational View on Bootstrap Ensembles as Bayesian Inference[END_REF]. Methods from this family use particles that are different instances of the model, which are optimized independently. The resulting set of optimized particles serves as an approximate posterior distribution that is more flexible than VI solutions. In this work we focus on variational bootstrap [START_REF] Milios | A Variational View on Bootstrap Ensembles as Bayesian Inference[END_REF], which provides a theoretical connection between POVI and parametric bootstrap [START_REF] Efron | Bayesian Inference and the Parametric Bootstrap[END_REF] for regression tasks, and our contribution is to extend this method to classification problems.

variational bootstrap for classification models

Considered Models

bayesian neural networks In this chapter, we will use Bayesian MLP described in Section 2.2.4, expression (2.44). The method that we propose is quite general and does not depend on the choice of the BNN architecture. It can easily be applied to convolutional and recurrent networks.

random fourier features approximation of gaussian processes

We shall also examine the case of linear models, as they can serve as scalable approximations to another class of Bayesian models, namely Gaussian processes [START_REF] Edward Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF]. Following the random Fourier features (RFF) approximation [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF], we consider ϕ(x i ) ∈ R D×1 to be the projection of an input point x i ∈ R d×1 onto a feature space of D trigonometric basis functions. Then Φ ∈ R D×N denotes the design matrix of the entire training set in the feature space. In this case, the model parameters can be directly described as a vector: θ := w ∈ R m . Given a Gaussian likelihood N (y; f (x, w), σ 2 ) and a Gaussian prior over the weights w ∼ N (0, α 2 I m ), then the posterior distribution over w after observing the dataset D is known to be Gaussian, yielding the following predictive mean and variance for a test point x * :

E[ f (x * )] = 1 σ 2 ϕ(x * ) ⊤ A -1 Φy, Var[ f (x * )] = ϕ(x * ) ⊤ A -1 ϕ(x * ), where A = 1 σ 2 ΦΦ ⊤ + 1 α 2 I. (4.1)
We observe that in order to make a prediction, one has to solve two D × D linear systems: A -1 Φy is solved only once, but A -1 ϕ(x * ) has to be solved for every new test point x * . So for n * test points, this translates to O(n * × D 2 ) complexity.

If the number of test points is large, then it is preferable to directly calculate the decomposition of the matrix A, so that it can be reused to solve the linear systems needed to calculate the predictive distribution for new test points. In many problems however, a large number of features might be required to obtain an accurate approximation of a GP. Later, we show that variational bootstrap can keep the computational cost down when both D and n * are large.

Variational Bootstrap for Neural Network Regression

In its original formulation, variational bootstrap was defined as an alternative to Bayesian inference for regression tasks. More specifically, the objective is to approximate the posterior distribution of a nonlinear model with a Gaussian likelihood and a Gaussian prior as follows:

p(D|w) = n ∏ i=1 N (y i ; f (x i , w), σ 2 ) and p(w) = N (0, α 2 I m ). (4.2)
For this model, variational bootstrap yields a set of samples that represents an empirical distribution q and approximates the posterior distribution of the parameters w. These samples are obtained by optimization of a set of particles. Each particle is defined as the maximum a posteriori (MAP) estimate for a regression task on a perturbed version of the joint log-likelihood, with its own set of training labels and mean parameters of the prior distribution. The training labels for each particle are obtained by a parametric bootstrap procedure. For each given label y i , we generate a perturbed label ỹi according to the likelihood function, that is Gaussian with variance σ 2 and mean y i :

ỹ(k) i ∼ N y i , σ 2 , k = 1...K (4.3)
Also, each particle is associated with a unique sample from the prior distribution w(k) . So the new prior for each model in the ensemble becomes as follows:

p(w (k) , w(k) ) ∼ N ( w(k) , α 2 I m ), where w(k) ∼ N (0, α 2 I m ), k = 1...K (4.4)
After resampling K perturbed sets of the labels and parameters of the prior distribution, we can obtain K samples from the approximate posterior by solving K optimization problems of the form:

argmin w (k) 1 2σ 2 n N ∑ i=1 ( ỹ(k) i -f (x i , w (k) )) 2 + 1 2α 2 ||w (k) -w(k) || 2 , k = 1...K (4.5)
The parameters of the particles are updated by a gradient descent algorithm. For the case of NNs under the additional assumption of linear or piecewise linear activation functions, it is shown in [START_REF] Milios | A Variational View on Bootstrap Ensembles as Bayesian Inference[END_REF] that each gradient step optimizes the joint log-likelihood of each model and moves the distribution of parameters q(w) closer to the true posterior p(w|D).

The main benefit of this procedure is that it allows to use modern gradient optimization techniques, like Adam, to approximate the posterior distribution of the parameters of a nonlinear model.

Variational Bootstrap for Random Fourier Features

We demonstrate here that variational bootstrap can induce computational advantages also for the RFF approximation of Gaussian processes. If the labels are perturbed according to the likelihood so that ε ∼ N (0, σ 2 ), and the regularization term is a sample from the prior so that w ∼ N (0, α 2 I m ), then the MAP solution is:

ŵ = 1 σ 2 A -1 Φ(y + ε) + 1 α 2 A -1 w, (4.6)
The MAP estimate ŵ is a Gaussian random vector, as it is a linear combination of two Gaussian random perturbations: ε and w. According to [START_REF] Milios | A Variational View on Bootstrap Ensembles as Bayesian Inference[END_REF], we can calculate the expectation and the covariance of ŵ to obtain the true posterior mean and covariance for the weights of the linear model:

E ε, w[ ŵ] = 1 σ 2 A -1 Φy =: w, E ε, w[( ŵ -w)( ŵ -w) ⊤ ] = A -1 (4.7)
Equation (4.6) yields samples from the true posterior distribution. From a computational perspective, we observe that it behaves differently from the predictive posterior in (4.1). Here, exactly two linear systems have to be solved for every sample, as opposed to every test point as in (4.1). Let K denote the number of samples; then the complexity becomes O(K × D 2 ). If the number of test points is larger than the number of posterior samples, this can induce significant computational gains, as we demonstrate in the experimental section.

Dirichlet Label Transformation

The goal of Bayesian classification is to estimate the distribution of class probabilities for an input data point. The Gaussian likelihood model we discussed in the previous section is not appropriate for a classification task; it is more reasonable to use a Multinomial likelihood instead. For a C-class classification problem, the class label y for the input point x is a sample from the Multinomial distribution y ∼ Cat(π). The authors of [START_REF] Milios | Dirichlet-based Gaussian Processes for Largescale Calibrated Classification[END_REF] propose to use a C-dimensional Dirichlet distribution to model the distribution of class probabilities π ∼ Dir(α), with parameters α = [α 1 , ..., α C ] T .

p (π 1 , . . . , π C ; α 1 , . . . , α C ) = 1 B(α) K ∏ i=1 π α i -1 i , B(α) = ∏ K i=1 Γ (α i ) Γ (α 0 ) (4.8)
Any label observations are treated as Dirichlet distributions: if an input point x belongs to a class k, then it corresponds to the following Dirichlet parameters:

α i =    1 + α ϵ , if i = k. α ϵ , if i ̸ = k. (4.9)
The term α ϵ > 0 represents a small quantity added to the Dirichlet parameters in order to guarantee a valid Dirichlet distribution. Then it is possible to represent samples from the C-dimensional Dirichlet distribution as samples from C Gamma distributions:

π i = z i ∑ C c=1 z c
, where z i ∼ Gamma(α i , 1), for i ∈ {1, ..., C}. The original paper proposes to approximate the Gamma distribution by a Lognormal( ỹi , σ 2 i ) distribution whose parameters are determined by moment matching. Because the logarithm of the log-normally distributed random variable has a Gaussian distribution N ( ỹi , σ 2 i ), it becomes possible to use a Gaussian likelihood and thus transform the classification problem into a regression problem: the transformed labels ỹi become the new targets for the inputs x i and σ 2 i becomes the variance parameter of the Gaussian likelihood.

σ 2 i = log(1/α i + 1) ŷi = log α i -σ 2 i /2 (4.10)
It is worth mentioning that the expression (4.10) produces different noise parameters σ 2 for each observation y. This means that after the transformation, we are dealing with heteroskedastic linear regression.

In order to train the whole classification model, it is required to solve C regression problems, one for each dimension of the Dirichlet distribution. To make predictions, one has to apply a softmax transformation to the outputs of the C regression models f = [ f 1 , ..., f C ] T , as follows:

E[π i |x] = exp( f i (x)) ∑ C c=1 exp( f c (x) p( f i (x)|X)df(x) (4.11)
During the inference stage this approach consideres the new type of transformation, which is different from the common sigmoid or softmax transformations. This choice is motivated by the opportunity of approximating the posterior distribution of the transformed labels ỹi with a Gaussian distribution.

Classification with Variational Bootstrap

One of the key components of variational bootstrap for regression is the data resampling via parametric bootstrap. In the case of classification, it is not straightforward to adjust this strategy to produce perturbed versions of the class labels in a way that reflects the nature of the Bernoulli (or the Multinomial) likelihood. For example, a class label y can take values in {0, 1}; if we locally fit a distribution Bern(p) to each y, the maximum-likelihood parameter will be the one-sample mean, i.e., p = 0 or p = 1. If we use this fitted model (i.e. Bernoulli with parameter 0 or 1) to resample new labels, this will deterministically produce either 0 or 1, depending on the original label.

Therefore, we have adopted a strategy that combines variational bootstrap with the Dirichlet labels transformation. Consider a dataset D = {x i , y i } N i=1 of input vectors x i and labels y i for a C-class classification problem. As a first step, we transform each label y i into a pair { ŷi , σ i }, ŷi ∈ R C , σ i ∈ R C by the means of the Dirichlet label transformation. The C-class classification problem is then transformed into a C-dimensional heteroskedastic regression problem with labels ŷi and observation noise variances σ 2 i ; each training sample and each dimension of the output has observation noise with its own variance.

The second step involves an application of the variational bootstrap method, where we generate K independent sets of regression labels, as well as K sets of prior parameters. Perturbations of each transformed label ŷi from the training set are sampled from the corresponding distribution with variance σ 2 i as follows:

ỹi ∼ N ( ŷi , Diag(σ 2 i )) where Diag(σ 2 i ) = I ⊙ (σ 2 i 1 ⊤ ) (4.12)
Algorithm 1 Variational bootstrap for MLP Draw sample w from N (0, α 2 I)

9:

Initialize w (k) ← w 10:

w (k) ← for each output dimension optimize (4.5) 11: end for Algorithm 2 Variational bootstrap for RFF 1: Input: Φ, y, α, α ϵ 2: Output: {w (1) , . . . , w (K) } ∼ p(w|Φ, y)

3: for i ← 1 to N do 4: ŷi , σ i ← DirichletTransform(y i , α ϵ ) 5: end for 6: for k ← 1 to K do 7: ỹ(k) 1 , . . . ỹ(k) N ∼ N ( ŷ1 , σ 2 1 ), . . . , N ( ŷN , σ 2 N ) 8:
Draw sample w from N (0, α 2 I)

9:
w (k) ← for each output dimension solve (4.6) 10: end for

The procedure for applying variational bootstrap to classification is summarized in Algorithm 1 for MLP models and in Algorithm 2 for the RFF Gaussian process approximation. In both cases, parameters of the prior distribution for each particle are sampled according to (4.4), while the labels are perturbed as in (4.12). Then in Algorithm 1, we optimize K models independently. For Algorithm 2, the main difference is that we obtain the solution of regression problem (4.5) analytically instead of using gradient optimization. After optimization of the K particles, it is possible to perform predictions using a Monte Carlo approximation of Equation (4.11). As a final remark, our choice to employ the Dirichlet label transformation and to treat the classification problem as regression implies that our methods enjoy any convergence guarantees that can be proven for the regression case.

Experiments: Toy dataset

We first demonstrate the application of variational bootstrap on a synthetic onedimensional binary classification problem. We considered a Gaussian process classifier based on a radial basis function (RBF) kernel, and we approximated it with a Bayesian logistic regression model on a set of random features. We used Random Fourier Features approximation of the RBF kernel proposed in [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF]. For our experiment, we considered 5, 000 random features. On the left panel of Fig. 4.1, we show the distribution over functions corresponding to the approximate posterior distribution over the parameters obtained with variational bootstrap with the Dirichlet label transformation; optimization was performed via the L-BFGS algorithm. In the middle panel of the figure, we report the distribution over functions obtained by the same approximation of the model, but where inference is carried out by Markov chain Monte Carlo (MCMC). On the right panel of the figure, we show the distribution of functions obtained with the Laplace approximation. For the toy dataset, we used the Metropolis-Hastings algorithm with 100 chains.

The Laplace approximation is described in Section 2.2.3.

For the prediction, we took the last sample from each chain. R-hat convergence diagnostic [START_REF] Stephen | General Methods for Monitoring Convergence of Iterative Simulations[END_REF] showed that it takes around 10 6 steps before convergence. The comparison shows a remarkable property of the proposed approach to accurately approximate the posterior over model parameters without the need for expensive or excessively long computations. In fact, the L-BFGS algorithm converged for variational bootstrap after 16 iterations only. We evaluated our method on several UCI classification problems outlined in Table 4.1. We applied variational bootstrap with the Dirichlet label transformation on two different models. First, we considered a two-hidden layer MLP model with a ReLU activation function and 512 neurons in each hidden layer. We used Adam [START_REF] Diederik | Adam: A Method for Stochastic Optimization[END_REF] to optimize the parameters of the model. This model is referred to as VBoot-MLP. The second model we considered is the RFF approximation [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF] of a GP, where 5, 000 random features were used to approximate an RBF kernel with fixed hyperparameters. In this case, we generated 50 samples as prescribed in The performance of variational bootstrap is compared against a regular BNN, whose posterior has been approximated by means of MCMC sampling, and in particular, SGHMC [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF]. The baseline BNN uses the same prior distribution over the parameters and a Bernoulli likelihood, while its output is given by a softmax activation function. VBoot-RFF is compared against the Dirichlet transformation with a Sparse GP approach proposed in [START_REF] Milios | Dirichlet-based Gaussian Processes for Largescale Calibrated Classification[END_REF]. The hyperparameters of this model were optimized as in the original paper. We evaluate performance using several metrics, including classification error and mean negative log-likelihood (MNLL). The results are outlined in Tables 4.2 and 4.3, respectively. The models considered are referred to as VBoot-MLP, VBoot-RFF, Sparse GP, and MCMC-MLP, respectively.

Experiments: UCI Datasets

Considering the MLP models, our experiments show that the proposed approach is competitive when compared against a principled framework such as MCMC. Most importantly, our approach has significantly different behavior in terms of convergence speed, as it relies on optimization rather than sampling. We demonstrate this property empirically by monitoring the progression of validation error for variational bootstrap and MCMC. In Fig. 4.2, we report classification error on held-out data over the first 200 training epochs. These results indicate that the variational bootstrap technique provides much faster convergence compared to SGHMC. On the other hand, only in some cases SGHMC eventually settles to significantly lower validation error. We used the Wilcoxon test [START_REF] Wilcoxon | Individual Comparisons by Ranking Methods[END_REF] to compare classification errors and MNLLs for VBoot-MLP and MCMC-MLP methods. The test did not show any statistically significant difference regarding the performance of these two methods within 0.05 significance level. The lack of statistical significance appears partly due to the small number of used data splits. Nevertheless, variational bootstrap is shown to achieve a good trade-off between accuracy and efficiency. This trade-off can be further exploited in practice, as variational bootstrap is trivially parallelizable.

Regarding the GP-based models, we also see that VBoot-RFF is highly competitive against traditional sparse GPs that rely on inducing points. In fact, according to the one-sided Wilcoxon signed-rank test, the VBoot-RFF results are slightly better than Sparse GP in a statistically significant way for some of the datasets (marked with "*" in Tables 4.2 and 4.3). Our variational bootstrap framework allowed us to use a large number of random features (i.e., 5, 000), which would not be possible for many of the datasets considered. This can be seen in Fig. 4.3, which shows the computation time for VBoot-RFF and Sparse GP models as a function of the number of random features and inducing points correspondingly. Of course, using more random features/inducing points results in a better approximation of the full GP model, but it also increases computational complexity. We see that VBoot-RFF has better scalability properties compared to sparse GPs. Regarding the computation times reported for Sparse GPs, we note that we have excluded the initial K-Means step needed to initialize the inducing locations, so in practice, sparse GPs are more expensive than what we report here. Interestingly, sparse GPs could not scale beyond 1, 000 or 2, 000 inducing points for some datasets due to memory errors. In this work, we proposed a novel way to carry out Bayesian inference for classification models based on Neural Networks and Gaussian processes. For NNs, this is important because, while they achieve state-of-the-art performance in many tasks, they lack a principled way to characterize uncertainty in predictions, so they represent a class of models for which a Bayesian treatment is highly desirable but mathematically and computationally challenging. For Gaussian processes, instead, while the Bayesian treatment is at the core of its formulation, classification tasks are difficult to handle because they require expensive approximations. Our work provides a practical and easily parallelizable way to tackle all these limitations. We are currently investigating parallel implementations of the proposed approach to considerably accelerate inference of large-scale problems by operating on clusters of computing machines. In addition, we are exploring the application of our approach to problems involving optical-based computing hardware, also known as Optical Processing Units [START_REF] Saade | Random projections through multiple optical scattering: Approximating Kernels at the speed of light[END_REF]. OPUs offer a fast and low-power way to approximate Gaussian processes through random features, and because they are capable of generating millions of these at the speed of light, we believe that our approach could be the key to exploit these computations effectively.

Another possible direction is the adaptation of our framework towards a more efficient marginal likelihood maximization for GPs, which is a standard practice to tune kernel hyperparameters [START_REF] Edward Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF]. In this work, we have treated GPs by means of fixed feature map approximations, which correspond to fixed hyperparameter values. Estimating marginal likelihoods through samples is an open research question, and it can be the subject of future work.

improving the dirichlet-transformation 4.4.1 Motivation

We aim to develop a method that provides well-calibrated predictions and accurate uncertainty estimates for Gaussian Processes Classification (GPC) by improving on the idea proposed in [START_REF] Milios | Dirichlet-based Gaussian Processes for Largescale Calibrated Classification[END_REF] and described in Sec. 4.3.4. The original method proposes to use a Log-normal distribution to approximate a Gamma distribution, which is used for sampling from a Dirichlet distribution. In fact, the Log-normal distribution poorly approximates Gamma distributions with a small value of the parameter α of the Gamma distribution. The experiments in the original paper [START_REF] Milios | Dirichlet-based Gaussian Processes for Largescale Calibrated Classification[END_REF] show that the best performance of the original method is achieved exactly for such cases. This means that an approximate posterior predictive distribution poorly approximates the true posterior predictive distribution when the classification error is the smallest. Also, this conclusion is indirectly supported by the visualization of the predictive posterior for 1D synthetic datasets, provided in [START_REF] Milios | Dirichlet-based Gaussian Processes for Largescale Calibrated Classification[END_REF], where we can see an underestimation of the variance of the predictive distribution.

The proposed work replaces the Log-normal approximation of the Dirichlet distribution used in the original work. The idea is to construct a transformation that allows to transform a Dirichlet distributed random variable to a normally distributed random variable. With this transformation, it will be possible to project the Dirichlet distributed class probability π to a latent representation in such a way that the latent representation is Gaussian. To obtain a class probability distribution for a test data point, it is necessary to solve a Gaussian process regression in a latent space with projected training labels. Then we have to apply the inverse of the initial transform to the prediction for the test data point representation in the latent space.

Cumulative Density Function Transformation

Considering a two-class classification problem, let's assume that there is a dataset D{(x i , y i )} of N input vectors x i and binary labels y i . In this case, the distribution of class probabilities for each input x is one-dimensional. Thus, it is possible to use a Beta distribution to model class probabilities instead of a Dirichlet distribution. For each training data point, there is a Beta distributed random variable π responsible for a class probability. If y = 0, the corresponding distribution of π will be Beta(α ϵ , 1 + α ϵ ), if y = 1, the corresponding distribution of π will be Beta(1 + α ϵ , α ϵ ). The hyperparameter α ϵ refers to the uncertainty in class labels.

The goal is to transform class probability π to a normally distributed latent random variable f centered at -µ for one class and µ for another class. We impose this symmetry in the latent space to reflect the symmetry that exists between Beta(α ϵ , 1 + α ϵ ) and Beta(1 + α ϵ , α ϵ ).

F Beta(α ϵ ,1+α ϵ ) (π) = 1 -F Beta(α 1+ϵ ,α ϵ ) (1 -π) (4.13)
After the transformation, the symmetric one should be as follows:

F N (-µ,σ 2 ) ( f ) = 1 -F N (µ,σ 2 ) (-f ) (4.14)
It is known that it is possible to transform a random variable X with a cumulative density function F X to a random variable Y with a cumulative density function F Y as follows:

Y = F -1 Y (F X (X)) (4.15)
We can apply this expression to transform a Beta distributed random variable into a normally distributed random variable.

g 1 (x) = F -1 N (-µ,σ 2 ) (F Beta(α ϵ ,1+α ϵ ) (x)) g 2 (x) = F -1 N (µ,σ 2 ) (F Beta(1+α ϵ ,α ϵ ) (x)) (4.16)
Also, the requirement of symmetry provides an intuition for another constraint. When we observe a latent prediction that lies in the middle between the centers of the Gaussians, we make a decision that the model is unable to assign a specific class label to this prediction. Thus, we have to force the inverse transformation to return 0.5 for this latent prediction. This will signify that the model assigns 0.5 probability for both classes.

g -1 1 (0) = 0.5 F -1 Beta(α ϵ ,1+α ϵ ) (F N (-µ,σ 2 ) (0)) = 0.5 (4.17)
The same derivations applied for the transformation g 2 will give the same result. Thus, we have obtained a connection between parameters µ and σ for the latent representations of the distributions of class probability.

µ = σ √ 2(erf -1 (2F Beta(α ϵ ,1+α ϵ ) (0.5) -1) (4.18)
This constraint gives the following property. Ratios between a probability mass of the whole distribution and the probability mass belonging to a region of overlap with the symmetric distribution are equal for the distributions of class probability and for their latent representations.

1 1 -F N (-µ,σ 2 ) (0) + F N (µ,σ 2 ) (0) = 1 1 -F Beta(α ϵ ,1+α ϵ ) + F Beta(1+α ϵ ,α ϵ ) (0.5) 1 2F N (-µ,σ 2 ) (0) = 1 2F Beta(α ϵ ,1+α ϵ ) (0.5) (4.19)
The latter equality is satisfied because of equality (4.17). Now there are two transformations instead of one. This fact poses the question: which one is necessary to use for the inverse operation? When we observe in the latent space a GP prediction with parameters -µ, σ 2 , it is better to use g 1 , and g 2 for the symmetric case. There is an intuition that it is better to use the transformation g -1 1 if the latent prediction is closer to the Gaussian, centered at -µ and the transformation g -1 2 , when the latent prediction is closer to the Gaussian centered at µ.

Both transformations g -1 1 , g -1 1 represent some activation functions that maps an interval [-∞, ∞] into [0, 1] (Fig. 4.4). There is an asymmetry of each function with respect to the point [0, 0.5]. This asymmetry imposes the following problem. If we apply both transformations to the latent predictions with a distribution N (-µ, σ 2 ) the goal is to obtain a distribution of class probabilities Beta(α ϵ , 1 + α ϵ ) for each transformation. It is possible to see that the transformation g -1 2 gives an incorrect distribution of class probabilities when we observe the distribution of the latent predictions, corresponding to y = 0, which is N (-µ, σ). The symmetric situation happens for y = 1 with g -1 1 transformation. It is possible to overcome the described difficulty by constructing a hybrid transformation that combines transformations g 1 and g 2

g(x) =    g 1 (x), if x < 0.5 g 2 (x), if x >= 0.5 (4.20)
The inverse transform is as follows:

g -1 (x) =    g -1 1 (x), if x < 0 g -1 2 (x), if x >= 0 (4.21)
This transformation is symmetric and thus produces almost correct distributions of class probabilities π for the extreme cases in the latent space when the GP predictions are distributed according to N (-µ, σ 2 ) or N (µ, σ 2 ). The quality of this approximation for non-extreme cases will be analyzed empirically in the next section.

Experiments

At the beginning, the method was tested on a toy classification dataset. Simple GP regression was used to solve the classification problem in a latent space. The hyperparameters of the kernel were optimized by log marginal likelihood maximization. Obtained predictions were compared with the predictions obtained with the log-normal approximation proposed in [START_REF] Milios | Dirichlet-based Gaussian Processes for Largescale Calibrated Classification[END_REF] (Fig. 4.6). The same noise hyperparameter α ϵ = 0.2 was used for both methods. As we see in Fig. 4.6, the proposed method assigns higher uncertainty estimates for the same input points. It is worth mentioning that for the toy dataset, the choice of µ hyperparameter for the latent space does not affect the distribution of class probabilities. Different values of parameter µ in a range [0.01, 50] were tested, and all of them gave exactly the same predictions in the space of class probabilities.

Also, the proposed method was tested on several UCI classification datasets and compared with the log-normal approximation in terms of classification error, mean negative log-likelihood (MNLL), and expected calibration error (ECE), which is outlined in Tables 4.4, 4.5, 4.6. In the tables, the proposed method is denoted as Hybrid Transformation, and the baseline is denoted as Log-normal approximation. The classification performances of both methods are very similar, but the proposed method requires solving only one regression problem in a latent space instead of two, as in the case of the baseline. The general goal of this work was to obtain an accurate posterior predictive distribution. The proposed method still uses an approximation, and it does not allow obtaining the true posterior distribution. But the most problematic part of this approach is multi-class classification because, in this case, we have to construct a C-dimensional activation function that maps logits to the [0, 1] interval. The The preliminary experiments conducted in this work for two class datasets show some advantages of the proposed technique compared to the Dirichlet label transformation in terms of negative log-likelihood and accuracy, but they are ambiguous in terms of Calibration Error. Because of these two arguments, this proposed approach is considered promising while requiring more investigation.

L O C A L G A U S S I A N P R O C E S S R E G R E S S I O N

introduction

Function estimation is a fundamental problem in Machine Learning. In supervised learning tasks applied to a data set composed of observed input data and labels, the goal of function estimation is to establish a mapping between these two groups of observed quantities. Function estimation can be approached in various ways, and we can broadly divide algorithms in two categories, as global and local. Examples of global algorithms are Neural Networks [START_REF] Radford | Bayesian Learning for Neural Networks[END_REF] and kernel machines [START_REF] Shawe | Kernel Methods for Pattern Analysis[END_REF], which impose a functional form yielding a global representation of the function. The functional form is parameterized by a set of parameters that are optimized or inferred based on all the available data. The estimated model can later be used to query the function at any input points of interest. In local algorithms such as K-Nearest Neighbors (KNN), instead, the target point is fixed, and the corresponding value of the function is estimated based on the closest data available.

Obviously, any global algorithm can be made local by training it only for the few training points located in the vicinity of the target test point. While it may seem that the idea of localizing global algorithms is not a very profound one, empirical evidence shows that localization could improve the performance of the best global models [START_REF] Bottou | Local learning algorithms[END_REF]. The idea of localization was therefore applied to global models such as SVMs [START_REF] Blanzieri | An adaptive SVM nearest neighbor classifier for remotely sensed imagery[END_REF][START_REF] Blanzieri | Nearest neighbor classification of remote sensing images with the maximal margin principle[END_REF]. In addition to performance gains, by operating on smaller sets of data points, these local approaches enjoy computational advantages, which are particularly attractive for kernel machines for which scalability with the number of data points is generally an issue [START_REF] Cheng | Localized support vector machine and its efficient algorithm[END_REF][START_REF] Segata | Fast and Scalable Local Kernel Machines[END_REF][START_REF] Segata | Local SVM approaches for fast and accurate classification of remote-sensing images[END_REF].

In this work, we develop novel ideas to implement a localization of Gaussian processes (GPs) in order to obtain performance gains, as well as computational ones. GPs are great candidates to benefit from computational speedups given that a naïve implementation requires expensive algebraic computations with the covariance matrix; denoting by n the number of input data, such operations cost O(n 3 ) operations and require storing O(n 2 ) elements, hindering the applicability of GPs to data sets beyond a few thousand data points [START_REF] Quinonero | A unifying view of sparse approximate Gaussian process regression[END_REF]. Another issue with GPs is how to choose a suitable kernel for the problem at hand so as to avoid problems of model misspecification. Both of these issues have been addressed in various ways, by proposing scalable approximations based on inducing points [START_REF] Hensman | Gaussian Processes for Big Data[END_REF] and random features [START_REF] Cutajar | Random Feature Expansions for Deep Gaussian Processes[END_REF][START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF], and by composing GPs to obtain a rich and flexible class of equivalent kernels [START_REF] Wilson | Deep Kernel Learning[END_REF].

In this work, we explore an alternative way to address scalability and kernel design issues by localizing GPs. In particular, we show how the localization operation leads to a particular form for the localized GP and what is the effect on the kernel of this model. Furthermore, the localization makes it apparent how to implement the model with considerable gains compared to other approaches to approximate GPs. We demonstrate such performance gains on regression tasks on standard UCI benchmarks [START_REF] Dua | UCI Machine Learning Repository[END_REF].

related work

Local learning algorithms were introduced by Bottou and Vapnik [START_REF] Bottou | Local learning algorithms[END_REF], with the main objective of estimating the optimal decision function for each single testing point. Examples of local learning algorithms include the well-known K-Nearest Neighbor regression [2] and local polynomial regression [START_REF] Fan | Local polynomial modelling and its applications[END_REF]. These methods provide simple means for solving regression problems for the cases where training data are nonstationary or their size is prohibitively large for building a global model. However, neither of these methods provides ways to quantify uncertainty in predictions, which is a highly desirable feature in cost-sensitive applications.

Gaussian Process Regression (GPR) [START_REF] Edward Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF] is a popular nonparametric regression method based on Bayesian principles, which provides uncertainty estimates for its predictions. Similarly to other kernel methods (e.g., SVMs and KRR), GPR is a global method, meaning that it takes into account the whole dataset at prediction time. Thus, GPR inherits the computational complexity of global kernel methods, which is prohibitive for large datasets. Among the large class of scalable approximations for GPR, successful ones are based on Random Fourier Features [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF] and on sparsification of the Gram matrix induced by the kernel [START_REF] Edward Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF].

Random feature approximation of the kernel proposed in [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF] is based on the Bochner theorem and allows representing the kernel function as a dot product of (possibly infinite) feature maps applied to the input data. In practice, infinite feature maps are replaced by a finite Monte Carlo approximation. The disadvantage of this approach is that it is necessary to construct a specific random feature mapping for each type of kernel. While random feature approximations are known for popular kernels such as RBF [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF], and polynomial [START_REF] Pennington | Spherical Random Features for Polynomial Kernels[END_REF], there is no straightforward application of this method to approximate arbitrary kernels.

The Gram matrix sparsification approach is based on the idea of introducing so-called inducing points in order to approximate the full Gram matrix. One of the most popular methods in this family is the Nyström approximation [START_REF] Edward Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF]. The main drawback of this approach is that a low number of inducing points might lead to a poor approximation of the original model, which affects predictive performance. An important advancement within this family of approaches which provides a scalable variational formulation, was proposed in [START_REF] Michalis | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF].

While providing good performance and scalability for large datasets, these approaches still require some design choices for the kernel. For stationary kernels, they assume that the same kernel is suitable for all the regions of input space, and if data are nonstationary, this may harm the predictive performance. The literature has a wide range of proposals to address kernel design by incorporating ideas from deep learning [START_REF] Cutajar | Random Feature Expansions for Deep Gaussian Processes[END_REF][START_REF] Wilson | Deep Kernel Learning[END_REF].

Recently, partitioning strategies have also gained some attention. The main idea is to divide the input space into regions where local estimators are defined [START_REF] Carratino | Park: Sound and efficient kernel ridge regression by feature space partitions[END_REF][START_REF] Meister | Optimal learning rates for localized SVMs[END_REF][START_REF] Müecke | Reducing training time by efficient localized kernel regression[END_REF][START_REF] Tandon | Kernel ridge regression via partitioning[END_REF]. In partition-based methods, the main challenge is to define an effective partitioning of the space.

There are several approaches that use the idea of local learning for training GP models. The method proposed in [START_REF] Meier | Local Gaussian regression[END_REF] and extended in [START_REF] Meier | Incremental local Gaussian regression[END_REF] mostly focuses on Bayesian parametric linear regression. The methods in these papers build an ensemble of local models centered at several fixed points, where each training point is weighted accordingly to the distance from the center of the model. Predictions are computed as a weighted sum of the local models. The authors claim that their approach extends to GPR, but in this case, each local model considers the full training set. This means that these methods use localization to address nonstationarity but poorly scale to large datasets. The method proposed in [START_REF] Snelson | Local and global sparse Gaussian process approximations[END_REF] proposes to build local GP models that use only subsets of the training data, but it lacks a mechanism that assigns importance weight for the training points for each local model according to the distance from the center of the model. That is why the model can make overconfident predictions for the points that lay far away from the centers of the local models. In [START_REF] Robert | Local Gaussian process approximation for large computer experiments[END_REF], in order to obtain fast approximate prediction at a target point, the Authors propose a forward step-wise variable selection procedure to find the optimal sub-design.

gaussian processes, kernel ridge regression, and localization

Gaussian Process Regression

We already discussed GP regression in Chapter 2, Section 2.2. The problem with this approach is that it is required solving a linear system involving a matrix of size n × n to obtain the parameters of a predictive distribution (2.28), (2.29). Direct methods to solve these operations require O n 3 operations and storing O n 2 . Iterative solvers, instead, can reduce these complexities by relying exclusively on matrix-vector products, which require O n 2 operations per iteration and do not need to store the Gram matrix [START_REF] Cutajar | Preconditioning Kernel Matrices[END_REF][START_REF] Filippone | Enabling scalable stochastic gradient-based inference for Gaussian processes by employing the Unbiased LInear System SolvEr (ULISSE)[END_REF]. However, a quadratic time complexity may still be prohibitive for large-scale problems.

There is rich literature on approaches that recover tractability by introducing approximations. One popular line of work introduces m so-called inducing points as a means to approximate the whole GP prior [START_REF] Quinonero | A unifying view of sparse approximate Gaussian process regression[END_REF]. This treatment of GPs was later extended within a scalable variational framework [START_REF] Krauth | Au-toGP: Exploring the Capabilities and Limitations of Gaussian Process Models[END_REF][START_REF] Michalis | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF], making the complexity cubic in the number of inducing points m. Another approach proposes ways to linearize GPs by obtaining an explicit set of features so as to obtain a close approximation to the original kernel-based model. Within this framework, a popular approach is based on random features [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF]. Denoting by Φ the n × D matrix obtained by applying a set of D random basis functions to the inputs in x 1 . . . , x n , these approximations are so that ΦΦ ⊤ approximates in an unbiased way the original kernel matrix K XX , that is E[ΦΦ ⊤ ] = K XX . For the Gaussian kernel, for example, a Fourier analysis shows that the basis functions that satisfy this property are trigonometric functions with random frequencies [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF]. This approach has been applied to GPs in [START_REF] Lázaro-Gredilla | Sparse Spectrum Gaussian Process Regression[END_REF] and later made scalable by operating on mini-batches in [START_REF] Cutajar | Random Feature Expansions for Deep Gaussian Processes[END_REF].

x 0 x 0 x 0 h h h (a) (b) (c)

Locally Smoothed Gaussian Process Regression

GP regression, as formulated above, is an example of a global learner; for a fixed training data D, it builds a posterior distribution over functions that can be used to calculate the predictive distribution for any test inputs. To construct a predictive distribution for any input point, it is necessary to use all the information available from the data D, resulting in the need to do algebraic operations with the matrix K XX + σ 2 I n . However, if we focus on the prediction problem locally, at a given target input x 0 , most of the information carried by the (potentially large) covariance matrix might be neglected with little loss of information. The main idea behind localized GPs is to down-weight the contribution of the data points far from x 0 , so that the structure of the covariance matrix is more adapted to the prediction task at a given point. To make this general idea work, we need to tackle two challenges. First, we need to specify what it means to be far or close to a given point; second, the change of the structure of the covariance matrix must give us a valid covariance matrix, i.e., the resulting covariance matrix should be symmetric and positive definite. Note that a simple truncation of the covariance function to obtain a compact-support covariance function may generally destroy positive definiteness [START_REF] Kaufman | Covariance Tapering for Likelihood-Based Estimation in Large Spatial Data Sets[END_REF].

We accomplish the localization of GPs in a straightforward manner as follows. We localize the target and the prior in the model (2.32) by multiplying them by the square root of the weighting function:

k h (x, x 0 ) := 1 h k ∥x -x 0 ∥ h , where k : X ⊂ R d → R is a non-negative, integrable function satisfying K(x)dx = 1 and ∥ • ∥ is Euclidean norm on R d .
Considering the square root of the weighting function will be convenient later when we discuss the link between local GPs and local Kernel Ridge Regression. Some classical examples of the weighting functions are given in Fig. 1. Because of the linearity of the weighting operation, the resulting model is another zero-mean Gaussian process f (x) with covariance function given by

K(x, x ′ ; x 0 ) = k 1 2 h (x, x 0 )K(x, x ′ )k 1 2 h (x ′ , x 0 ). (5.1)
In this formulation, we have localized the relationship between noisy targets and function realizations as with ỹi = k h (x i , x 0 )y i , and the prior is given by a zero-mean GP with the localized covariance kernel (5.1). The model (5.2) can be alternatively written as a model with heteroscedastic noise

ỹi = f (x i ) + ε i , ε i ∼ N (0, σ 2 ), (5.2)
-4 -2 0 2 4 -2 -1 0 1 2 (a) -4 -2 0 2 4 -2 -1 0 1 2 (b) -4 -2 0 2 4 -2 -1 0 1 2 (c) -4 -2 0 2 4 -2 -1 0 1 2 (d)
y i = f (x i ) + 1 k h (x i , x 0 ) ε i , ε i ∼ N (0, σ 2 ).
Making the noise parameter location-dependent can significantly improve the performance for problems where the assumption of a homoscedastic noise is not satisfied.

Proposition 1 Let I = {i : ∥x i -x 0 ∥ ≤ h}, X I = {x i : i ∈ I} and y I = {y i } i∈I ∈ R |I| .
Assume that (5.2) holds for the fixed target point x 0 . Then f (x 0 ) | y I is a Gaussian random variable with mean and variance given by

m(x 0 ) = K x 0 X I K X I X I + σ 2 W -1 x 0 -1 y I (5.3) K (x 0 , x 0 ) = K (x 0 , x 0 ) -K x 0 X I K X I X I + σ 2 W -1 x 0 -1 K X I x 0 (5.4)
where W x 0 is the diagonal matrix with main diagonal entries k h (x i , x 0 ), x i ∈ X I .

Proof: Let x 0 be any fixed target point. Then the observations y I ∈ R |I| and GP-function value at target point

f 0 = f (x 0 ) ∈ R are jointly Gaussian such that y f 0 ∼ N 0 I 0 , K X I X I + σ 2 W -1 x 0 K X I x 0 K x 0 X I K(x 0 , x 0 ) .
Then the proposition follows from the basic formula for conditional distributions of Gaussian random vectors (see, e.g., [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], Appendix A.2).

Compared to global GPs, in the local formulation, in order to compute the posterior mean and variance, we need to invert K X I X I + σ 2 W -1

x 0 . This might give a key advantage when dealing with large data sets, as the localization by the compactly supported kernel (local) could significantly sparsify the Gram matrix corresponding to K XX . Denoting by s 0 the number of inputs for which the localizing weights are nonzero for a test point x 0 , the complexity of performing such an inversion is O(s 3 0 ). Another interesting observation is that the kernel function K (x, x ′ ) is potentially more flexible than the original kernel function; this is due to the multiplication by the localizing weighting function, which may introduce some interesting nonstationarity even for kernel functions which are stationary, depending on the choice of the weighting function.

The calculation of the predictive distribution with Locally Smoothed Gaussian Process Regression (LSGPR) is described in Algorithm 3 The parameter selection in probabilistic models given by GPs is based on the marginal log-likelihood maximization, which, in our local formulation, can be defined as follows:

log p(y I |X I ) = - 1 2 y ⊤ I K X I X I + σ 2 W -1 x 0 -1 y I - 1 2 log K X I X I + σ 2 W -1 x 0 - n 2 log (2π) (5.5)
Unfortunately, gradient-based optimization cannot be used to find the optimal localization parameter h, as the marginal log-likelihood is not continuously differentiable w.r.t. this parameter when compactly supported local kernels are used.

The simplest way to resolve this problem is by using grid search for the localization parameter h, while kernel parameters can be optimized by gradient-based methods for any given h.

Algorithm 3 LSGPR

1: Input: X, y, σ 2 , h, x 0 2: Output: m(x 0 ), K(x 0 , x 0 ) 3:

I := {i : ||x i -x 0 || ≤ h} 4: X I := {x i : i ∈ I} 5: y I := {y i : i ∈ I} 6: W x 0 := diag({k h (x i , x 0 ) : i ∈ I}) 7: L := Cholesky K X I X I + σ 2 W -1 x 0 8: α := (L -1 ) ⊤ L -1 y 9: m(x 0 ) := K x 0 ,X I α 10: v := L -1 K x 0 ,X I 11: K(x 0 , x 0 ) := K(x 0 , x 0 ) -v ⊤ v

Local Kernel Ridge Regression

For every Gaussian process f (x) with covariance function K(x, x ′ ), there is a unique corresponding Hilbert space H K . This is commonly referred to as a reproducing kernel Hilbert space (RKHS) and constructed as a completion of the linear space of all functions:

x → k ∑ i=1 α i K (a i , x) , α 1 , . . . , α k ∈ R, a 1 , . . . , a k ∈ X , k ∈ N relative to the norm induced by the inner product k ∑ i=1 α i K (s i , •) , l ∑ j=1 β j K t j , • H K = k ∑ i=1 l ∑ j=1 α i β j K s i , t j .
It is well known that the posterior mean of Gaussian process regression can be alternatively derived by minimizing the regularized empirical risk over the RKHS [START_REF] George | A correspondence between Bayesian estimation on stochastic processes and smoothing by splines[END_REF]; see, e.g., [START_REF] Kanagawa | Gaussian processes and kernel methods: A review on connections and equivalences[END_REF] for a recent review. For local GPs, this corresponds to a weighted least square minimization over the RKHS with the weights given by k h (x, x 0 ), that is

m(x) = arg min f ∈H k n ∑ i=1 (y i -f (x i )) 2 k h (x i , x 0 ) + σ 2 n ∥ f ∥ 2 H k . ( 5.6)
Note that in the local formulation, for a given point x 0 one has to estimate both the parameters of the reproducing kernel and the width of the local kernel h. Here are two examples of well-known classical local methods, which are the solution to the empirical risk minimization problem (5.6).

k-nearest neighbors This model corresponds to the noise-free case (σ = 0) with a positive constant reproducing kernel and a rectangular local kernel whose width is adjusted to contain exactly k data points. The solution of the minimization problem (5.6) is the mean of the outputs corresponding to the k closest to x 0 input points.

local polynomial regression If we use the polynomial kernel K(x, x ′ ) = (1 + xx ′ ) k for the space H K , and use any smooth local kernel (i.e. exponential), then in the noise-free case the solution of the minimization problem (5.6) is so called local polynomial regression [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF]. In this special case, when the degree of the polynomial is 0, we have Nadaraya-Watson regression, which is the minimizer of the local squared loss over the constant function.

experiments

Toy dataset

In order to illustrate the behavior of the proposed Locally Smoothed Gaussian Process (LSGPs), we start from a toy dataset generated from the Doppler function (Fig. 5.3).

y(x) = x(1 -x) sin 2.1π x + 0.05 + ε, 0 ≤ x ≤ 1, ε ∼ N (0, 0.1) (5.7)
For this experiment, we used the RBF kernel, and the Epanechnikov localizing kernel

k(x) = 3 4 (1 -|x| 2 )I(|x| ≤ 1). (5.8)
We tuned the lengthscale parameter of the RBF kernel by optimizing the marginal log-likelihood of the model. We used the L-BFGS algorithm for gradient optimization [START_REF] Pytlak | Conjugate gradient algorithms in nonconvex optimization[END_REF]. We chose the value of the parameter h of the localizing kernel that gave Because our model is more flexible than standard GPR, we also added to the comparison Deep GP models based on random features expansion. The size and dimensionality of these problems are outlined in Table 5.1. Since the Euclidean norm used in the local kernels depends on the units in each coordinate, all the datasets except the Protein were scaled within the [0, 1] range. We used a standardization procedure for the Protein dataset because the baseline model worked much better with this type of preprocessing. For the Deep GP model, we used the hyperparameters and the data preprocessing described in the original paper.

The LSGPR method requires creating a new local model with its own set of hyperparameters for each input point where the prediction has to be made. During the optimization, the kernel parameters of each model were constrained to be equal among all local models. Considering the hyperparameter h, we found that it is hard to find values of h that perform well across all regions of the input space. Thus, for each input point of interest, we chose values of h that ensured that the localizer considers at least m neighboring training points. In this experiment, we used 3-fold cross-validation to choose the noise variance σ 2 , the lengthscale of the GP kernel, and the parameter m of the localizing kernel. We report the results on the held-out test set.

In this experiment, we also used the Hilbert localizing kernel [START_REF] Devroye | The Hilbert kernel regression estimate[END_REF][START_REF] Shepard | A two-dimensional interpolation function for irregularlyspaced data[END_REF] k(x) = ∥x∥ -1 I(∥x∥ ≤ 1), (

which showed good performance for most of the datasets. In Table 5.2, locally smoothed Gaussian Process Regression based on Hilbert kernel is referred to as LSGPR Hilbert, while the same model based on Epanechnikov kernel is referred to as LSGPR Epanechnikov. GP regression, Deep GP regression, and KNN regression are referred to as GP, DeepGP and KNN, respectively.

The results indicate that LSGP offers competitive performance with respect to GPR and Deep GP baselines. The results also clearly show that LSGP offers superior performance to KNN, suggesting that the localization alone is not enough to obtain good performance and that this works well in combination with the GP model. To make a comparison with the baselines, we used the one-sided Wilcoxon test [START_REF] Wilcoxon | Individual Comparisons by Ranking Methods[END_REF]. For each method, we measured its performance on 10 data splits, and we used exactly the same splits for testing the performance of each method, so we matched samples of the MSEs. Then we used the test to compare methods in pairs, where the alternative hypothesis was that the MSE of any given method is smaller than the MSE of competitors. We used confidence level α = 0.05. In Table 5.2, the results that are statistically better than the competitors are marked in bold.

discussion of the method

In this work, we developed a novel framework to localize Gaussian processes (GPs). We focused, in particular, on Gaussian Process Regression (GPR), and we derived the GP model after applying the localization operation through the down-weighting of contributions from input points that are far away from a given test point. The form of the localized GP maintains positive definiteness of the covariance, and it allows for considerable speedups compared to standard global GPR due to the sparsification effect of the Gram matrix. The proposed method may suffer from the curse of dimensionality. In this work we chose the localization hyperparameter to ensure, that the local model considers at least m neighboring training points. This is a reasonable choice for low-dimensional problems, but it may not be sufficient for high-dimensional problems.

The proposed method requires cross-validation to tune the scale parameter of the localizing kernel, while other GP-based techniques use a less expensive marginal log-likelihood (MLL) gradient optimization to tune these types of parameters. We found MLL gradient optimization problematic because of the discontinuity of the local kernel with respect to the scale parameter, which in turn makes the MLL function discontinuous with respect to this parameter. It would be interesting to investigate ways to extend the idea of localization for GPR to other tasks, such as classification.

I M P O S I N G F U N C T I O N A L P R I O R S O N B AY E S I A N N E U R A L N E T W O R K S

introduction

Artificial Neural Networks (NN) currently represent a general class of successful models for various machine learning tasks, including computer vision, natural language processing, and many others. Bayesian Neural Networks (BNN) combine the representation power of NNs with Bayesian inference, making them an attractive choice in applications where predictive performance and accurate uncertainty quantification in parameter estimates and predictions are both important. However, BNNs are typically difficult to work with due to the intractability of the posterior over model parameters, which requires one to resort to approximations. Furthermore, recent works point to the need to choose sensible priors over model parameters in order to be able to obtain good performance [START_REF] Vincent Fortuin | Priors in bayesian deep learning: A review[END_REF][START_REF] Ba-Hien Tran | All You Need is a Good Functional Prior for Bayesian Deep Learning[END_REF]. In BNNs, the prior over the weights and the network architecture determine a distribution over the outputs of such BNNs [START_REF] Sun | FUNC-TIONAL VARIATIONAL BAYESIAN NEURAL NETWORKS[END_REF], and we refer to this induced prior as a functional prior.

The functional prior should encode any prior information on the conditional distribution of the labels given the inputs. However, it is unclear how to encode this type of information when having to specify a prior distribution over the weights.

In this paper, we propose a novel framework to enforce such meaningful functional priors. In particular, we rely on scalable Markov chain Monte Carlo (MCMC) sampling from an approximation to the posterior distribution over BNN weights, and we specify the prior over the weights implicitly through a prior over the induced functional prior. Our approach is different from the literature on Implicit Process Priors (IPPs) [START_REF] Ma | Variational Implicit Processes[END_REF], where the goal is to obtain an approximate framework to handle the functional prior implicitly induced by choice of a prior distribution over the weights. In our work, we operate in the opposite direction by imposing a functional prior, which implicitly determines a prior over the weights; we do not know such a prior over the weights in closed form, but we implicitly determine it through the specification of the induced functional prior.

Stochastic Processes are natural mathematical objects suitable to define distributions over functions [START_REF] Kallenberg | Foundations of modern probability[END_REF], and GPs represent popular examples that are routinely used in numerous machine learning tasks. This type of stochastic process is well investigated and has strong theoretical foundations [START_REF] Edward Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF]. There are theoretical guarantees for the generalization error of GP regression, and this method has a strong connection with non-Bayesian Kernel Ridge Regression (KRRs) [START_REF] Kanagawa | Gaussian processes and kernel methods: A review on connections and equivalences[END_REF]. Also, it was shown in [START_REF] Radford | Bayesian Learning for Neural Networks[END_REF] that in the infinite width limit, shallow BNNs are equivalent to GPs. We propose to use GPs to impose functional priors over BNNs because GPs provide a flexible set of instruments for encoding different types of beliefs about functions, such as periodicity or smoothness, through the specification of kernels. However, our approach is not restricted to GPs, and it can handle any functional priors that can be written down in the closed form.

a review of types of prior in bnn

A popular way of choosing prior distributions for BNNs is to employ a Gaussian distribution over the weight of the model [START_REF] Graves | Practical Variational Inference for Neural Networks[END_REF][START_REF] Radford | Bayesian Learning for Neural Networks[END_REF]. While this choice of prior does not render inference tractable, it does offer some practical advantages, for instance, when employing Variational Inference (VI) [START_REF] Graves | Practical Variational Inference for Neural Networks[END_REF]. In this case, it is possible to calculate part of the VI objective without the need to resort to Monte Carlo sampling, and it requires a linear complexity in the number of parameters when making a mean-field assumption. This computational benefit comes with a cost of poor approximation of the true posterior due to the lack of flexibility of the family of the approximate distributions.

Even when adopting more advanced and generally more accurate inference techniques, such as Hamiltonian Monte Carlo (HMC) [START_REF] Radford | MCMC using Hamiltonian dynamics[END_REF], and Stochastic-Gradient HMC (SGHMC) [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF], the Gaussian assumption on the prior over model parameters is still very common. By studying the entropy of the predictive distribution, it was shown by [START_REF] Vincent Fortuin | Bayesian Neural Network Priors Revisited[END_REF] that Gaussian priors are problematic in terms of model performance and the ability to detect Out-of-Domain (OOD) input examples. This work also shows how Gaussian priors over the weights could be responsible for the cold posterior effect described by [START_REF] Wenzel | How Good is the Bayes Posterior in Deep Neural Networks Really?[END_REF]; this effect is characterized by the necessity of applying temperature scaling to the prior density term in Bayes theorem in order to obtain good performance.

The poor performance associated with the choice of Gaussian priors gave rise to alternatives involving flexible distributions, such as a mixture of Gaussians [START_REF] Blundell | Weight Uncertainty in Neural Network[END_REF], Student's t-distribution [START_REF] Vincent Fortuin | Bayesian Neural Network Priors Revisited[END_REF], hierarchical Gaussian distribution [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF], and many others [START_REF] Vincent Fortuin | Priors in bayesian deep learning: A review[END_REF]. However, all these types of prior distributions imposed over the weights of the model share the same issue of preventing from understanding the effect that they impose on the outputs of the model.

Another way of encoding the desired output prior constraints into the model is to use an auxiliary probabilistic model which learns the desired distribution of parameters to fulfill the required output constraints and to use this distribution as a prior over the weights of the original model [4].

Instead of focusing on the prior over the weights of NNs, an alternative is to study their effect on the distribution on the NNs output, which we refer to as functional priors. [START_REF] Sun | FUNC-TIONAL VARIATIONAL BAYESIAN NEURAL NETWORKS[END_REF] proposes to use VI to tackle the problem of finding a Bayesian posterior in the space of functions for a functional prior, defined by a stochastic process. They introduce a variational objective, which includes the evaluation of the KL divergence between the approximate predictive posterior and the functional prior using a finite set of evaluations of the function. It was shown that the supremum of the KL divergence over all possible sets of input points is equal to the true KL divergence in functional space. In this setting, the model is then trained in an adversarial manner, which means that the optimization procedure simultaneously minimizes the optimization objective with respect to the parameters of the model and maximizes the KL term with respect to the input data points, which makes the optimization process unstable. Also, the opti-mization objective requires evaluating the gradient of the approximate posterior density, which is performed with the Stein gradient estimator [START_REF] Shi | A Spectral Approach to Gradient Estimation for Implicit Distributions[END_REF], and this, in turn, requires a careful choice of a kernel function. The work in [START_REF] Ma | Variational Implicit Processes[END_REF] focuses on representing the functional prior as a BNN and uses GPs to obtain an approximate posterior over functions. The problem with this approach is that GPs may yield a poor approximation quality for the true functional posterior. The authors in [START_REF] Sun | FUNC-TIONAL VARIATIONAL BAYESIAN NEURAL NETWORKS[END_REF] and [START_REF] Ma | Variational Implicit Processes[END_REF] use VI to find an approximate posterior distribution, which means that the optimization objective contains a functional KL divergence term. However, in [START_REF] Burt | Understanding Variational Inference in Function-Space[END_REF][START_REF] Tim | Rethinking Function-Space Variational Inference in Bayesian Neural Networks[END_REF], it is claimed that the KL divergence between the functional approximate posterior and the GP process functional prior is problematic as it may go to infinity. On the other hand, they acknowledge that it does not mean that parametric models cannot approximate GPs well.

The authors of [START_REF] Ba-Hien Tran | All You Need is a Good Functional Prior for Bayesian Deep Learning[END_REF] propose to impose GP functional priors so as to constrain the parametric prior over the weights of BNNs. They propose several options to define such a prior, including Gaussian, Hierarchical Gaussian, and Normalizing Flows [START_REF] Rezende | Variational Inference with Normalizing Flows[END_REF]. They propose to optimize the parameters of the prior over the weights so as to minimize the Wasserstein distance between the outputs of the untrained BNN and samples from the GP prior. Then, the posterior over the weights is characterized by means of MCMC.

In our work, we aim to avoid the computation of the KL divergence or any other distance metric in functional space. Instead, we propose to enforce the choice of a functional prior directly when carrying out approximate inference of BNN weights.

imposing functional priors on bnns

Consider a supervised learning task with a dataset D{(x i , y i )} i=1...n of n input vectors X = {x i } i=1...,n and corresponding labels y = {y i } i=1...,n , and imag- ine employing a NN-based model with parameters w to establish a parametric mapping between inputs and labels. We denote the input/output mapping by f w (x), and for convenience we also define f ⊤ = [ f w (x 1 ), . . . , f w (x N )] and f * ⊤ = [ f w (x 1 ), . . . , f w (x N ), f w ( x1 ), . . . , f w ( xM )] as the evaluation of the function f w (x) at the inputs X and an augmented set of inputs X * = [X, X], respectively. The set X * has cardinality N * = N + Ñ, and the Ñ inputs in X are drawn from a given p(x). Note that the sets X and X * can be disjoint, but in order to keep the notation uncluttered, we assume X ⊂ X * A Bayesian treatment NNs requires specifying a prior distribution p(w) over the parameters and a likelihood function for the labels given the inputs, that is p(y|X, w). For this BNN, it is possible to write down an expression for the posterior distribution over model parameters as: p(w|y, X) = p(y|X, w)p(w) p(y|X, w)p(w)dw (6.1)

Carrying out inference in BNNs is extremely difficult for at least two reasons. One main difficulty stems from the complex way in which parameters affect the likelihood function, and this requires approximation techniques to characterize the posterior over model parameters; popular approaches involve MCMC and variational approximations. A second and more subtle challenge is how to specify priors for BNNs because it is difficult to establish what is the effect of prior parameters on the distribution over the functions that BNNs can represent. In this work, we propose a novel way to address the challenge of choosing sensible priors for BNNs by working with implicit priors over the weights induced by the choice of functional priors while we follow the recent trend to employ MCMC techniques to address the intractability of the inference process. We begin by focusing on the distribution over the functions represented by BNNs. In particular, we consider the distribution of f * , which is the distribution of f w (x) evaluated at the set of input points X * , and we impose a prior over this set of variables which encourages functions to behave in a sensible way a priori. Later we will study in particular Gaussian process priors, but any functional prior can be incorporated as long as it can be expressed in closed form.

We now rewrite the likelihood function in terms of f rather than w:

p(y|X, w) → p(y|f). (6.2)
The main idea behind our work is to now define a prior over f instead of w, and to perform inference over w. With this change of variables, we should account for the change of measure through a Jacobian term. However, such a change of variables involves groups of variables of different dimensions in general, and even when this is not the case, computing this term would be computationally costly. For this reason, we are going to ignore the Jacobian accepting to settle for an approximate posterior over w. With this choice, we rewrite Bayes theorem as:

log p(f * |y, X * ) = log p(y|f) + log p(f * |X * ) + const. (6.3)
Note that in this equation, we introduced the functional prior:

p(f * |X * ) = p(f * |X * , w)p(w)dw, (6.4) where p(f * |X * , w) is a Dirac's delta placed at the evaluation of f w (x) at the inputs X * due to the deterministic way in which inputs are mapped into outputs in NNs. Again, we stress that while we focus on the distribution of functions represented by BNNs, we actually use the objective in eq. 6.3 to perform MCMC sampling in the space of the weights w. Note that we carry out inference over w through MCMC, but given that we are working with an approximation to the posterior over w, we could alternatively employ other fast approximate inference techniques such as Variational Inference. Here, we focus on MCMC so as to isolate the effect of the way we impose functional priors compared with alternatives that try to characterize the exact posterior over w [START_REF] Ba-Hien Tran | All You Need is a Good Functional Prior for Bayesian Deep Learning[END_REF].

bayesian interpretation. From a Bayesian point of view, imposing a prior over function by specifying a prior over f * induces an implicit prior over the weights through eq. 6.4. In other words, the prior over f * is, in practice, a prior over a deterministic transformation of w, and this is implemented by the NN. It is interesting to note that in the literature eq. 6.4 is usually interpreted in the opposite way; that is, one uses eq. 6.4 starting from a prior over the weights p(w) to define a functional prior in an implicit way [START_REF] Ma | Variational Implicit Processes[END_REF]. The likelihood function establishes what is the likelihood of the labels y, and it is conditioned on f or equivalently on w and X. Therefore, the expression in eq. 6.3 can be seen as an expression for the (approximate) posterior over the weights w (due to the lack of a Jacobian term), where the prior is assumed over a transformation of such weights. In this paper, we take this view to carry out Bayesian inference over w using MCMC techniques. We also note that our approach has some close similarity with the Product of Expert approach proposed in [START_REF] Wenk | Fast Gaussian process based gradient matching for parameter identification in systems of nonlinear ODEs[END_REF] for inference of parameters of Ordinary Differential Equations using Gaussian Processes.

regularization interpretation. While we proceed with a Bayesian treatment of w, it is useful to interpret eq. 6.3 as a regularized objective in the following way. The first term log p(y|f) is the negative loss, which can be equivalently seen as a function of w and X * , so this provides a constraint on w because the objective promotes values of f which are compatible with the labels y, and f depends on w and X * . The second term is a regularization term, which penalizes functions deviating from a behavior established by the functional prior. Because f * is a function of w and X * , this translates into a regularization term for w.

gaussian process priors

The proposed formulation focusing on functional representations has the advantage of putting the emphasis on the functions that BNNs can represent, and for which it is possible to assume sensible priors. Here we specify how to operate in the case of Gaussian processes (GPs), which yield a prior term in eq. 6.3 as:

log p(f * |X * ) = - 1 2 f * ⊤ C -1 f * + const, (6.5)
where the covariance matrix is C = (K X * X * + σ 2 n ), and K X * X * contains the evaluation of the kernel function κ among all the inputs in X * . For simplicity, we assumed a zero-mean GP, but other mean functions can be easily included.

In the next subsections, we elaborate on how to use this GP prior in practice by proposing a way to operate with mini-batches for scalability purposes, by discussing hyper-parameter optimization, and by discussing the properties of the proposed approach when N * goes to infinity.

Mini-batching

In this work, we aim to employ advanced MCMC sampling methods based on stochastic gradients, and in particular Stochastic Gradient Hamiltonian Monte Carlo (SG-HMC) [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF] to sample from the weights w of BNNs. In order to do so, we need to formulate our MCMC objective in a way that is suitable for mini-batching. However, extending the previous formulation to operate with minibatches without care would produce a biased estimation of the quadratic term

f * ⊤ C -1 f * ̸ = E[f * ⊤ b C -1 b f * b ],
where f b and C b are computed over a mini-batch X b . The main difficulty of full batch training is the necessity of solving linear systems with the matrix C, which has O(N * 3 ) complexity in the number of inputs in X * . The literature on GPs offers many cues on how to circumvent this problem. In particular, there exist formulations of GPs based on inducing points [START_REF] Hensman | Gaussian Processes for Big Data[END_REF] and random features [START_REF] Lázaro-Gredilla | Sparse Spectrum Gaussian Process Regression[END_REF][START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF] which operate on mini-batches [START_REF] Cutajar | Random Feature Expansions for Deep Gaussian Processes[END_REF]. In this work, we focus on approximations based on random features, but inducing points formulations is also possible.

The Random Features (RF) expansion of the kernel κ(•, •) allows us to obtain a finite-dimensional representation for an explicit feature map that approximates its value. Using this expansion, we can express the Gram matrix as a dot product of feature maps computed over the data K ≈ ΦΦ ⊤ . We can use this property and the Woodbury identity to rewrite the quadratic term as follows:

f * ⊤ C -1 f * = f * ⊤ (ΦΦ ⊤ + σ 2 f I) -1 f * = 1 σ 2 f f * ⊤ f * - 1 σ 2 f f * ⊤ Φ(Φ ⊤ Φ + σ 2 f I) -1 Φ ⊤ f * . (6.6)
In this case, instead of inverting the matrix of size N * × N * , we invert the matrix of size D × D, where D is the dimensionality of the RF vector. But this approach has two drawbacks. First, it is extremely unstable when σ 2 f → 0, because after application of the Woodbury identity the term 1

σ 2 f f * ⊤ f * → ∞.
Second, this approach still does not allow mini-batch training.

We can reformulate our MCMC objective by replacing the nonparametric term pertaining to the GP with a parametric one based on RFs. For the set f * , we can factorize its prior probability as:

p(f * |X * ) = p(f * |β, X * )p(β)dβ, (6.7)
where β are the parameters of RF approximation of the GP, that is p(β) ∼ N (0, I) and p(f * |β, X * ) ∼ N (Φβ, σ 2 f I). In this case it is easy to verify that p(f * ) = N (0, ΦΦ ⊤ + σ 2 f I) and according to the property of the RF approximation, the covariance matrix coincides with the prior term of the objective in eq. 6.6. Instead of sampling directly from the unnormalized posterior p(f * |X * , y) marginalized over β, we can sample from the joint density p(f * , β|X, y) and discard samples over β:

p(f * , β|X * , y) ∝ p(y|f)p(f * |β, X * )p(β). (6.8)
Again, when we refer to the fact that we sample f * , in practice, we sample w. This RF-based approach avoids the necessity of inverting the matrix (ΦΦ ⊤ + σ 2 f I) during the computation of the objective.

Resuming, the expression for the unnormalized log-posterior in eq. 6.3, where the GP regularization is approximated using RFs, is as follows: log p(f * , β|X * , y) = log p(y|f)-

1 2σ 2 f ∥f * -Φβ)∥ 2 - ||β|| 2 2 + const. (6.9)
It is straightforward to verify that this MCMC objective can be written as a sum of terms involving individual input points, and it is therefore, amenable to minibatching. It is also easy to verify that one can proceed with a Gibbs sampling scheme whereby f * (that is w) is sampled from the conditional p(f * |β, X * , y) using SG-HMC and β is sampled directly from p(β|f * , X * , y), which has a Gaussian form.

Hyper-parameter optimization

The choice of a GP prior opens to the need to specify its kernel parameters. In the absence of any way to determine such hyper-parameters, we propose to optimize them by marginal log-likelihood (MLL) optimization, which is a popular way to proceed with GP models. In our case, the random feature approximation lends itself to a scalable solution, avoiding the need to invert large matrices. In particular, the marginal likelihood is:

p(y|X) = N (y|0, ΦΦ ⊤ ). (6.10)
Again, using Woodbury matrix identities, it is possible to rewrite the marginal likelihood so that the cost of computing it is cubic in the number of random features instead of cubic in the number of input points.

Classification

While for regression, it is natural to specify functional priors through GPs and to obtain a tractable framework to scale these through random features, for other likelihoods, things may become more involved. For instance, in classification problems, we may wish to specify functional priors such that the distribution over classes is uniform a priori. Alternatively, following an empirical Bayes approach, we could optimize the GP prior hyper-parameters so as to maximize the marginal likelihood. In this case, the random feature approximation of GPs leads to so-called Generalized Linear Models (GLMs) and this requires approximations to be able to compute the marginal likelihood. For classification tasks, there exist solutions to bypass the need to work directly with Bernoulli or Multinoulli likelihoods p(y|w, X).

Here we follow the idea proposed by [START_REF] Milios | Dirichlet-based Gaussian Processes for Largescale Calibrated Classification[END_REF], in which labels are transformed so that classification models can be replaced by regression models with heteroscedastic observation noise. In particular, for each one-hot encoded label y we can obtain real-valued vectors ỹ, σ 2 n (see [START_REF] Milios | Dirichlet-based Gaussian Processes for Largescale Calibrated Classification[END_REF] for details): .11) In this case, we can use a Gaussian likelihood, and thus we can obtain a closedform solution for the marginal likelihood of the model using the fact that the Gaussian likelihood and Gaussian prior are conjugate.

ỹi = log(α i ) - σ2 i 2 ; σ2 i = log 1 α i + 1 . ( 6

Connections with Kernel Ridge Regression

As already pointed out, in the proposed approach, we are free to choose the set X * . In this section, we try to get insights as to what is the behavior of the Gaussian process functional prior as a function of the cardinality of X * , that is N * .

In particular, we aim to show that the term associated with the functional prior in eq. 6.3, which is the quadratic form f * ⊤ C -1 f * , does not grow indefinitely with N * .

To do so, we resort to some connections with Kernel Ridge Regression. According to [START_REF] Veit Wild | Connections and Equivalences between the Nyström Method and Sparse Variational Gaussian Processes[END_REF] we can express

f * ⊤ C -1 f * = min g∈H 1 σ 2 f N * ∑ i=1 ( f w (x * i ) -g(x * i )) 2 + ||g|| 2 H .
(6.12)

In the case when f w belongs to the RKHS induced by the kernel κ of the GP prior, we can easily see that the quadratic term on the LHS of (6.12) is bounded by ∥ f ∥ 2 H . Indeed, if ĝ is a unique minimum of the RHS of eq. 6.12, we can write 1

σ 2 f N * ∑ i=1 ( f w (x * i ) -ĝ(x * i )) 2 + || ĝ|| 2 H ≤ 1 σ 2 f N * ∑ i=1 ( f w (x * i ) -f w (x * i )) 2 + || f w || 2 H = || f w || 2 H . (6.13) This means that if || f w || 2 H is bounded, the quadratic term f * ⊤ C -1 f * is bounded as well.
As an illustration of this fact, we picked a 1D function that has a following form:

f (•) = M ∑ i α i κ(•, x i ). (6.14)
We used M = 20, x i ∼ Unifrom(-1, 1), α i ∼ Unifrom(-1, 1). For this function, we can compute the Hilbert norm analytically. Then, for this function, we evaluated the value of the quadratic term f * ⊤ C -1 f * for different amounts of points x i (Fig. 6.1). We can also reuse a result from [START_REF] Veit Wild | Connections and Equivalences between the Nyström Method and Sparse Variational Gaussian Processes[END_REF], which allows us to express the projection of f w ∈ H on a subset of H, H ′ ⊂ H. And H ′ is a linear span of functions

ϕ i (•) = k(•, x i ), i ∈ {1, . . . , N * }. P( f w ) = arg min g∈H ′ || f w -g|| 2 H = k X * (•)K -1 X * X * f * , (6.15)
and we can note that

∥P( f w )∥ 2 H ′ = f * ⊤ K -1 X * X * f * . (6.16)
By increasing the number of basis functions in H ′ to infinity we get

min g∈H ′ || f w -g|| 2 H = || f w -P( f )|| 2 H N→∞ ---→ 0. ( 6
.17)

Using the expression for P( f ), we obtain:

|| f w -P( f )|| 2 H = || f w || 2 H + f * ⊤ K -1 X * X * f * -2⟨ f w , k X (•)K -1 X * X * f * ⟩ = || f w || 2 H -f * ⊤ K -1 X * X * f * . (6.18)
In the last step, we used the reproducing property of the RKHS H. Substituting this result into (6.17) we obtain:

f ⊤ K -1 X * X * f N→∞ ---→ || f w || 2 H . (6.19)
In cases where the function f w is represented by NNs with smooth and bounded activation functions, like tanh or sigmoid, at the layer before output, the output function of the NN f w is smooth and bounded for any bounded w. This is because a linear combination of smooth and bounded functions with bounded coefficients is smooth and bounded too. This means that if we are using kernels that have universal approximation properties, like RBF, exponential or binomial kernels [START_REF] Steinwart | Support vector machines[END_REF], each function generated by such NNs can be modeled by functions belonging to the RKHS imposed by this kernel.

With a finite amount of points in X * , the overall expression for the unnormalized log-posterior becomes:

log p(f|y) = log p(y|f) - 1 2 || ĝ|| 2 H + const, ĝ = arg min g∈H ′ || f w -g|| 2 H . (6.20)
This means that, in fact, we are using an approximation of the correct unnormalized log-posterior, but still, this objective guarantees that for the points in X * , we will obtain evaluations of functions sampled from the true posterior.

6.5 experimental evaluation of mcmc for gp priors

Toy regression dataset

We begin by testing our approach on a synthetic 1D dataset. We employ a twohidden layer NN with tanh activation function and with 256 neurons at each layer.

For the functional GP prior, we use an RBF kernel with a length-scale parameter l = 1 and output variance σ 2 out = 1. In Fig. 6.2, we report functions sampled from the predictive posterior of the BNN with such a GP prior (GP in the figure) along with the same GP prior approximated with 100 random features with and without mini-batching (GP RFF and GP RFF mini-batch in the figure). As a regularization set X, we used a grid of 200 equally spaced points. These points were used as a test set. For comparison, we also include the approach proposed in [START_REF] Ba-Hien Tran | All You Need is a Good Functional Prior for Bayesian Deep Learning[END_REF] where the prior over the weights is optimized so that the functional prior of the BNN minimizes the Wasserstein distance to the GP prior (WDGPi-G in the figure).

UCI regression datasets

To evaluate the effectiveness of our approach on more challenging problems, we consider a suite of UCI datasets [START_REF] Dua | UCI Machine Learning Repository[END_REF]. For these tests, we use a two-hidden layer NN with tanh activation function and with 100 neurons at each layer for all data sets, except for the Protein dataset, for which we used 200 neurons. We imposed a GP prior with an RBF kernel. Input vectors and labels were initially standardized.

For all experiments, we used the extended dataset X * , which contained 90% of real trained data and 10% of uniformly sampled vectors from the input domain.

The uniform distribution has support bounded at each dimension by minimum and maximum values of the training samples over this dimension. For all the datasets, we used full-batch training, while for Kin8nm, Power, and Protein, we used mini-batch training with a batch-size size of 512. As a baseline, we consider the method proposed in [START_REF] Ba-Hien Tran | All You Need is a Good Functional Prior for Bayesian Deep Learning[END_REF], and we choose the best-performing version from their paper, which is called GPi-G. We denote this as WDGPi-G here to emphasize the fact that this approach optimizes the Wasserstein distance of the BNN functional prior to the GP prior to determine the prior over the BNN weights. This method uses a Gaussian distribution as a parametric form for the prior over the weights and a Hierarchical GP with a LogNormal distribution over the length-scale and output variance of the GP kernel. The comparison of our method and WDGPi-G in terms of RMSE is shown in Table 6.1. We also report comparisons with deep ensembles [START_REF] Balaji Lakshminarayanan | Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles[END_REF], where each model in the ensemble had the same architecture as the NN used in the proposed method.

From these results, we see that our method outperforms WDGPi-G for most of the datasets. It is particularly interesting to note that WDGPi-G uses a form of Hierarchical Gaussian Process as a functional prior, whereas we consider a simple GP in our case. Hierarchical GPs represent a richer functional prior than standard GPs used in our method, but we still obtain competitive performance.

Toy classification dataset

We now illustrate how the proposed approach behaves on a 2D toy example involving the so-called banana dataset. Again we used a two-hidden layer NN with tanh activation function and with 256 neurons at each layer. We transform the labels using the approach from [START_REF] Milios | Dirichlet-based Gaussian Processes for Largescale Calibrated Classification[END_REF] to obtain labels which allow us to employ a Gaussian likelihood, as explained in Sec. 3. We then use an RBF kernel σ out = 5 with different values for the lengthscale and show how these choices affect the final results. As a regularization set X, we used a grid of 40×40 equally spaced points. Also, these points were used as a test set. For comparison, we report results with the WDGPi-G approach from [START_REF] Ba-Hien Tran | All You Need is a Good Functional Prior for Bayesian Deep Learning[END_REF]. As we can see on the plot, WDGPi-G method fails to incorporate the GP prior for a smaller lengthscale (l=0.1) and the prediction function behaves more smoothly than it should.

UCI classification datasets

In this section, we test the effectiveness of our approach on various UCI classification datasets. We use a two-hidden layer NN with tanh activation function. For the EEG, HTRU2, Letter, and Magic, we use 100 neurons in each hidden layer, while we use 200 neurons for Miniboo, Drive, and Mocap. Because all the datasets are rather large, in all cases, we use the Random Feature approximation of the functional GP prior with D = 1000 random features and mini-batches of size 512. GP hyper-parameters are optimized using the marginal likelihood after the label transformation presented in Sec. 3. We found that using the same transformation for the classification task with the BNN itself gave slightly better results than using classification likelihoods, and therefore we report these results in the table. We attribute this to the fact that the GP hyper-parameters are optimized with the transformed labels.

As baselines, we chose the same methods as in the regression experiments. The results in terms of classification accuracy and mean negative test log-likelihood are presented in Tables 6.3 and 6.4. The results show that even in the classification setting, the performance of the proposed approach is competitive with the state-ofthe-art. While we perform optimization of GP hyper-parameters, our approach allows us to get rid of the Wasserstein optimization phase, which is used in WDGPi-G, without a significant loss in classification performance.

limitations of the approach

While we consider our approach quite elegant in encoding prior information in the form of functional priors, we believe that it is important to point out some limitations compared to other works.

First, we point out that the form of posterior distribution we are targeting is approximate. Then, the main limitation is that the functional prior density needs to have a closed form. Even though the class of functional priors which have this property is large, this might be too restrictive in applications where it is possible to sample from such priors, but no closed form is available. Prior works which perform a preliminary optimization of the prior over the weights (e.g., [START_REF] Ba-Hien Tran | All You Need is a Good Functional Prior for Bayesian Deep Learning[END_REF]) can operate on samples from functional priors without the need to express these in closed form.

Another limitation related to the choice of a GP functional prior is how to set hyper-parameters. In this work, we resort to marginal likelihood optimization, but it is possible that this choice induces overfitting. One way around this would be to include hyper-parameters in the set of variables to be sampled in SG-HMC to obtain samples from their posterior at the expenses of having to deal with a more costly MCMC sampling. Having said that, there are situations where functional priors are easy to elicit and express without the need to carry out hyperparameter optimization.

discussion of the mcmc inference for functional priors

In this paper, we proposed a novel way to incorporate prior knowledge in Bayesian NNs (BNNs) in the form of functional priors. In our view, such functional priors implicitly determine priors over BNN weights, and the proposed formulation yields an approximate posterior over the weights from which it is possible to sample through MCMC or any other approximate inference techniques. In this paper, we studied the scenario where functional priors are expressed in the form of Gaussian processes (GPs), but our formulation can handle any functional prior, which can be expressed in closed form. We then discussed how to scale our approach to handle large data sets by operating on mini-batches, despite the complications stemming from the use of GP priors.

We can also point out that if we look at the proposed method from a different perspective, we find out that it produces samples from the true posterior distribution p(w|X * , y) of another probabilistic model, which is depicted on Fig. deterministically depends on w, so

p(f 1 |w) = δ(f 1 -f (X * , w)), (6.21)
where f (X * , w) is an output of the BNN. Let's derive a joint distribution of variables w, β, f 1 , f 2 , y:

p(w, β, f 1 , f 2 , y|σ n , σ f ) = p(w)p(f 1 |w)p(y|f 1 , σ n )p(β)p(f 2 |β, σ f )p(f 1 , f 2 ) = p(w)p(f 1 |w)p(y|f 1 , σ n )p(β)p(f 2 |β, σ f )δ(f 1 -f 2 ).
In the third line we used the fact that f 1 and f 2 are deterministically connected, thus p(f 1 , f 2 ) = δ(f 1f 2 ). Using the property of the Dirac delta function, we can note that:

p(f 2 |β, σ f )δ(f 1 -f 2 ) = p(f 1 |β, σ f ).
Also, we can substitute the expression for the p(f 1 |w) and obtain

p(w, β, f 1 , y|σ n , σ f ) = p(w)δ(f 1 -f (X * , w))p(y|f 1 , σ n )p(β)p(f 1 |β, σ f ).
Then, using the same property of the Dirac delta function p(w, β, y|σ n , σ f ) = p(w)p(y| f (X * , w), σ n )p(β)p( f (x, w)|β, σ f )

Here we note that this expression is exactly the sampling objective that was used in the proposed method, with an assumption that the prior distribution over the weights of the BNN p(w) is uniform. Thus, we can conclude that the proposed method actually produces samples from the true posterior distribution but from a different type of model with certain assumptions. We tested our proposal on regression and classification tasks and compared it with state-of-the-art approaches to carry out inference and prior optimization for BNNs. Our results demonstrate that the proposed approach is competitive in terms of performance and quantification of uncertainty while being easy to implement.

We are currently investigating ways to handle GP priors with priors over hyperparameters for increased flexibility and alternative ways to specify functional priors. Furthermore, we are investigating applications of BNNs for image classification tasks for which BNN architectures use convolutional layers.

variational bootstrap for bnns with functional priors

In this section, we develop another type of VI for BNNs with a GP prior over the outputs of the model. As it was mentioned before [START_REF] Burt | Understanding Variational Inference in Function-Space[END_REF][START_REF] Tim | Rethinking Function-Space Variational Inference in Bayesian Neural Networks[END_REF], the concept of explicit utilization of KL divergence as an objective for functional VI, like in [START_REF] Sun | FUNC-TIONAL VARIATIONAL BAYESIAN NEURAL NETWORKS[END_REF], [START_REF] Wang | Function Space Particle Optimization for Bayesian Neural Networks[END_REF], is questionable. One of the possible ways to avoid this problem is to abandon the explicit KL computation and switch to methods that optimize the divergence implicitly. Recently, a family of Particle Optimization Variational Inference (POVI) methods that exploit exactly this trick [START_REF] Francesco | Repulsive Deep Ensembles are Bayesian[END_REF], [START_REF] Liu | Stein Variational Gradient Descent: A General Purpose Bayesian Inference Algorithm[END_REF] was proposed.

POVI allows one to obtain samples from the posterior of the BNN without modification of the model. To achieve this property, methods from this family use particles which are different instances of the model. The particles are optimized independently, with special constraints for each particle. These constraints provide convergence of a set of particles to an approximate posterior distribution while preserving their diversity.

The performance of POVI methods depends on a choice of a kernel, which is independent of a type of prior in general. This work aims to establish some connection between a kernel used to construct a prior distribution of functions and a kernel used to solve a POVI problem.

The proposed work takes the idea of Variational Bootstrap (VB) from [START_REF] Milios | A Variational View on Bootstrap Ensembles as Bayesian Inference[END_REF]. We already described this idea in Section 4.3.2. In the proposed work, the diversity of the particles is introduced by perturbations of the training labels and different samples from the functional prior distribution. Mahalanobis regularization, proposed in Section 6.4, ensures the closeness of the set of particles to samples from the functional prior distribution. The set of perturbed labels and the samples from the functional prior are unique for each particle.

Linear model case

If there is a dataset D = {x i , y i } N i=1 , a likelihood of the model p(D| f ) and a functional prior GP, our goal is to get samples from the predictive posterior distribution of the model with a parametric bootstrap procedure.

Let's assume, that the functional prior is GP {0, κ} with some kernel κ,a likelihood of the model is Gaussian p(y|f) = N (y|f, σ 2 n ), where f is a vector of latent function values and σ 2 n is an observational noise. To get samples from a functional posterior, we can apply a functional version of the parametric bootstrap. Each particle represents a kernel ridge regression (KRR) problem with perturbed data labels ỹ ∼ N (y, σ 2 ) and regularization || f -f || 2 H , where f ∼ GP (0, κ).

argmin f 1 σ 2 n N ∑ i=1 ( ỹi -f (x i )) 2 + || f -f || 2 H (6.22)
The solution for each particle is as follows:

f (x * ) = k ⊤ (K + I) -1 ( ỹ -f) + f (x * ) (6.23)
where f is a vector of evaluations of the function f on the training set. It is easy to see that So, the mean and the variance of the predictions of the particles are equivalent to the mean and the variance of the GP prediction. We can observe that So, we can claim that g T K -1 g is the Hilbert norm of the smoothest function that passes through all the points f -f. So, this quadratic form gives a lower bound on || f -f || 2 H . For the KRR, we can reformulate the optimization objective as:

E f ∼p(0,κ) [ f (x * )] = k ⊤ (K + σ 2 n I) -1 y (6.
argmin f 1 σ 2 n N ∑ i=1 ( ỹi -f (x i )) 2 + (f -f) T K -1 (f -f) (6.27)
This optimization problem will have exactly the same solution as the initial formulation of the KRR. We can note that the regularization term is equivalent to a squared Mahalanobis distance. To avoid numerical problems, we can relax the equality constraint for (6.26)

h ′ = argmin g || g|| 2 H s.t. || g(x) -g(x)|| 2 < c, ∀x ∈ {x} N i=1
The optimization objective for this problem is as follows: Table 6.5: RMSE for GP regression (GPR), the proposed method (FBVI) and WDGPi-G

h ′ = argmin g 1 λ N ∑ i=1 ( g(x i ) -g(x i )) 2 + ||

Emprical evaluation

In the beginning, the Functional Bootstrap Variational Inference (FBVI) was tested on a toy regression dataset. A fully connected network was used for the experiments on synthetic data. The network had two hidden layers with 64 neurons in each layer and used tanh activation function. An equally spaced grid of input points was used for the Mahalanobis regularization term. Then the proposed method was tested on several UCI regression datasets. A network with the same architecture was used. It had 100 neurons in each layer, and it used the ReLU activation function. We performed the experiments only for the small UCI datasets, where full-batch training is available because the proposed method in its current state is not adopted for mini-batch training.

Discussion of the FVBI method

The proposed idea gave good-looking predictions on toy data (Fig. 6.5, that are much more similar to a GP than the ones produced by WDGPi-G method (Fig. 6.

2).

But the FBVI method failed to outperform the WDGPi-G in the preliminary experiments on the UCI data, and it demonstrated worse results than MCMC for GP functional priors proposed earlier in this chapter (Section 6.3). Besides, this method does not require a Wasserstein distance minimization phase like WDGPi-G. And as we see in Chapter 4, Fig. We proposed a unified approach for integrating OPUs into deep learning models that allows addressing two issues posed by this device. On the one hand, the proposed approach allows for dealing with the restrictions that OPUs put on the input data type. The ability to work exclusively with binary input data made OPU applicable only to a very narrow set of problems, which was a significant obstacle to the commercialization and widespread use of the technology. The proposed approach allows binarizing any input data of any type using an appropriate neural network in an optimal fashion regarding the loss of the specific machine learning problem. On the other hand, the proposed method allows integrating OPUs into deep learning models so that the whole pipeline can be trained in an End-to-End fashion. It allows for solving problems that could not be solved with the Kernel Ridge Regression method implemented by the device. As we emphasized in Chapter 1, machine learning became so successful due to feature extraction procedures trainable with backpropagation. And nowadays, it is the only way to achieve high performance for many machine learning tasks. That is why it was crucial to resolve the problem of non-differentiability of the OPU. While this work proposes solutions for some problems of the OPU, there are still challenges that have to be addressed. One of the most problematic things with optical hardware is the high latency of the data exchange between the OPU and other devices, such as CPU, GPU, and RAM. This drawback significantly alleviates the benefits of the OPU. The approach proposed in this work necessitates a considerable exchange of data between the device and the GPU. By conducting the data exchange in a batch manner, the proposed method prioritizes high bandwidth over low latency. However, the high latency of the OPU is still a challenge because it results in a slower training process compared to models that only use the GPU. Thus, one of the possible directions for future research could be decreasing the number of data exchanges between OPUs and other hardware. Another issue that requires further consideration is the possibility of stacking several OPUs in one pipeline. As it was said, the main benefit of optical hardware is the ability to extract features in nanoseconds, so if OPUs work on top of deep learning pipelines, the primary time complexity will come from GPUs. Moving as many computations as possible to OPUs will significantly decrease the inference time for the whole pipeline.

Efficient and accurate inference for Bayesian classification (Chapter 4)

Neural networks are widely used for achieving state-of-the-art performance in various tasks. However, they lack a systematic way to characterize uncertainty in predictions. This makes a Bayesian treatment of this type of models highly desirable but also mathematically and computationally challenging. On the other hand, GPs are founded on Bayesian principles, but classification tasks using GPs require computationally expensive approximations. Inspired by the idea of variational bootstrap, which allows performing fast and accurate Bayesian inference for regression models, like GPs, we proposed a novel method for conducting Bayesian inference in classification models such as neural networks and Gaussian Processes. Our method provides a practical solution to the abovementioned limitations by offering an easily parallelizable approach for Bayesian inference.

One of the possible directions for further research on the proposed method is various parallel implementations of our approach to significantly speed up inference for large-scale problems through the use of clusters of computing machines. Another possible improvement for the proposed method is to find a better way to approximate a Dirichlet likelihood used in the approach. Currently, the approximation is not accurate for the cases when the level of confidence in training labels is high. This may lead to the degradation of the model's performance in terms of classification accuracy and uncertainty estimation. We made some steps to improve the quality of approximation in Section 4.4, but this problem requires further investigation. And the third possible research direction for the proposed method is combining it with the methods presented in Chapter 3 for training Bayesian models on OPUs. The proposed method is preferable for the OPU setting because it provides a faster inference procedure compared to other Bayesian inference techniques.

Making Gaussian Processes scalable with localization (Chapter 5)

In Chapter 5, we introduced a new approach for scaling up GPs. We focused on regression problems and created a novel GP-based model by applying the localization technique. We achieved the localization effect by reducing the contributions of training points that are distant from a target test point. The training points outside the radius, defined by the localization hyperparameter, do not contribute to the prediction. The resulting GP maintains the positive definiteness property for a covariance matrix and allows for a significant performance improvement compared to a vanilla GP model that uses the full training set for each test point, thanks to the effect of sparsification of the Gram matrix. As a result, the model uses different subsets of training data to perform predictions for test points located in different regions of the input space.

Currently, the main drawback of the Local Gaussian Process Regression is the necessity to choose the localization hyperparameter with cross-validation. It is possible to leverage the marginal log-likelihood optimization approach to select the appropriate values for this hyperparameter. In a current setting applying this approach to the proposed method is problematic because the MLL function is not smooth with respect to the localizing hyperparameter. Addressing this issue could be one of the directions for future research.

Another direction for future research is extending this approach to random feature approximation of the localized GP. Some steps in this direction were already made in [START_REF] Wacker | Local Random Feature Approximations of the Gaussian Kernel[END_REF], but this work addresses a very special combination of GP kernel and localizing kernel. A more general way of obtaining such approximations for localizing GPs is an open question.

Bayesian neural networks with functional priors (Chapter 6)

In this chapter, we proposed a novel approach for incorporating functional prior knowledge into Bayesian neural networks using GP priors. Our approach connects the prior distribution over functions with a prior distribution over BNN weights. Using such type of priors allows for encoding into the deep model some functional properties, like smoothness, seasonality, or trend, which adds an additional mechanism to control the behavior of the deep model, which is often considered as a black-box from an interpretability point of view.

Our method allows for estimating the approximate posterior distribution of the weights using MCMC. Also, we showed the possibility of using a variational bootstrap approach for the same purpose. In the proposed work, we used GPs as a functional prior, but our formulation allows for extending the method to any type of functional distribution that has an analytical expression for the probability density. For the case of GP prior, we managed to make our approach scalable to large datasets with mini-batched training by using random feature approximation to address the scalability issues typical for GPs.

As it was mentioned in Chapter 6, the proposed method performs approximate posterior inference, and the main challenge for future research is making this procedure exact. The most promising direction, in this case, is POVI methods that, under some conditions, can provide an exact solution [START_REF] Francesco | Repulsive Deep Ensembles are Bayesian[END_REF].

Another possible extension of the proposed approach is speeding up the inference by variational bootstrap. We made some steps in this direction in Section 6.8. In this case, we are unable to obtain an exact posterior distribution over the weights, but our experiments in Chapter 4 showed that the variational bootstrap provides a satisfactory approximate posterior in terms of the performance of the model and its uncertainty estimation.

concluding remarks

In conclusion, this work explored some intersections between deep learning, Bayesian methods, and linear models, with a particular focus on Gaussian Processes. In four separate projects, we demonstrated how these approaches can be combined to enhance their flexibility, improve inference speed and scalability, and incorporate desirable functional properties into a model.

Deep learning showed outstanding results in complex machine learning problems, but its utilization in risk-sensitive applications is arguable due to the lack of uncertainty quantification and control over the properties of its predictions. Bayesian methods allow addressing both of these issues by forcing it to produce a predictive distribution instead of single predictions and by introducing prior knowledge into the model. The first property allows reasoning about the level of confidence of the predictions, and the second property allows forcing the model to respect some given constraints, which can be essential for some specific task.

As we explored in this thesis, Bayesian linear models can be modified to stay relevant in the age of big data on their own and also to work together with deep models, compensating for some of their drawbacks. The latter gives hope that the extensive theoretical framework behind these classical models will play some role in understanding the theoretical properties of deep models.

25 3 . 6 Figure 3 . 1 24 Figure 3 . 2 25 Figure 3 .3 26 Figure 3 . 4 29 Figure 3 . 5 31 Figure 3 . 6 31 Figure 3 . 7 34 Figure 4 . 1 43 Figure 4 . 2 46 Figure 4 . 3 46 Figure 4 . 4

 3631243225326342935313631373441434246434644 Figure 3.1 Initialization of the Method of Auxiliary Coordinates (MAC) procedure 24 Figure 3.2 The state of the MAC algorithm after one iteration 25 Figure 3.3 26 Figure 3.4 Mean Squared Error (MSE) for regression (top) and negative error on classification (bottom) datasets comparison. 29 Figure 3.5 Error comparison on classification (bottom) datasets for experiments on the real hardware. 31 Figure 3.6 Convergence of the training procedure on regression problem: boston dataset (left) and classification problem: mocap dataset (right). 31 Figure 3.7 Training loss with and without variance reduction 34 Figure 4.1 Comparison between the predictions with the regression weights, obtained by the variational bootstrap (left), MCMC (middle) and Laplace approximation (right) 43 Figure 4.2 Convergence of the classification error on validation data for variational bootstrap and SGHMC methods 46 Figure 4.3 Computation time complexity for VBoot-RFF and Sparse GP as a function of the number of random features and inducing points correspondingly. 46 Figure 4.4 The plot of the activation functions g -1 1 , g -1 2 for µ = 12.74, σ = 10 48 Figure 4.5 The histograms of predictions in the latent space (right) and corresponding transformed values in the space of class probabilities (left). 49 Figure 4.6 Toy dataset predictions obtained with the proposed transformation (top) and with the log-normal approximation (bottom) 50 Figure 5.1 Examples of local kernels: (a) Rectangular kernel, (b) Epanechnikov kernel, (c) Gaussian kernel. 55 Figure 5.2 (a) Samples from a global GP prior with the exponential kernel. (b) Samples from a GP prior with exponential kernel localized by rectangular smoother centered at x 0 = 0. (c) Samples from a GP prior with exponential kernel localized by Epanechnikov smoother centered at x 0 = 0. (d) Samples from a GP prior with exponential kernel localized by Gaussian smoother centered at x 0 = 0. 56 Figure 5.3 Illustration of the predictive distribution of GPs (left) and LSGPs (right) applied to data sampled from the Doppler function for 400 training points (top), 200 training points (middle), and 100 training points. 59

Figure 4 . 6 50 Figure 5 . 1 56 Figure 5 . 3 59 ixFigure 6 . 1

 46505156535961 Figure 6.1

Figure 6 . 2 71 Figure 6 . 3 73 Figure 6 . 4 75 Figure 6 . 5 FBVI predictions of 20 particles 79 L

 6271637364756579 Sampled predictions of BNNs where the GP functional prior is imposed implicitly (our work) and by means of the optimization of the Wasserstein distance with the functional BNN prior (WDGPi-G). Sampled predictions of the neural network with GP prior using Mahalanobis regularization and WDGPi-G methods Graphical model representation of an alternative view on the proposed method. Double line connections are used to indicate a deterministic relationship. I S T O F TA B L E S

Figure 3 . 1 :

 31 Figure 3.1: Initialization of the MAC procedure

Figure 3 . 2 :

 32 Figure 3.2: The state of the MAC algorithm after one iteration

Figure 3 . 3

 33 Figure 3.3

Figure 3 . 4 :

 34 Figure 3.4: Mean Squared Error (MSE) for regression (top) and negative error on classification (bottom) datasets comparison.

Figure 3 . 5 :Figure 3 . 6 :

 3536 Figure 3.5: Error comparison on classification (bottom) datasets for experiments on the real hardware.

Figure 3 . 7 :

 37 Figure 3.7: Training loss with and without variance reduction

Figure 4 . 1 :

 41 Figure 4.1: Comparison between the predictions with the regression weights, obtained by the variational bootstrap (left), MCMC (middle) and Laplace approximation (right)

Figure 4 . 2 :Figure 4 . 3 :

 4243 Figure 4.2: Convergence of the classification error on validation data for variational bootstrap and SGHMC methods

Figure 4 . 4 :

 44 Figure 4.4: The plot of the activation functions g -1 1 , g -1 2 for µ = 12.74, σ = 10

Figure 4 . 5 :

 45 Figure 4.5: The histograms of predictions in the latent space (right) and corresponding transformed values in the space of class probabilities (left).

Figure 4 . 6 :

 46 Figure 4.6: Toy dataset predictions obtained with the proposed transformation (top) and with the log-normal approximation (bottom)

Figure 5 . 1 :

 51 Figure 5.1: Examples of local kernels: (a) Rectangular kernel, (b) Epanechnikov kernel, (c) Gaussian kernel.

Figure 5 . 2 :

 52 Figure 5.2: (a) Samples from a global GP prior with the exponential kernel. (b) Samples from a GP prior with exponential kernel localized by rectangular smoother centered at x 0 = 0. (c) Samples from a GP prior with exponential kernel localized by Epanechnikov smoother centered at x 0 = 0. (d) Samples from a GP prior with exponential kernel localized by Gaussian smoother centered at x 0 = 0.

  f * C 1 f * Value of the Hilbert norm

Figure 6 . 1 :

 61 Figure 6.1: The value of the quadratic regularization term with respect to the number of points used for evaluation of the function.

Figure 6 . 2 :

 62 Figure 6.2: Sampled predictions of BNNs where the GP functional prior is imposed implicitly (our work) and by means of the optimization of the Wasserstein distance with the functional BNN prior (WDGPi-G).

Figure 6 . 3 :

 63 Figure 6.3: Sampled predictions of the neural network with GP prior using Mahalanobis regularization and WDGPi-G methods

6

 6

. 4 .Figure 6 . 4 :

 464 Figure 6.4: Graphical model representation of an alternative view on the proposed method. Double line connections are used to indicate a deterministic relationship.

  [START_REF] Cho | Kernel Methods for Deep Learning[END_REF] andvar( f (x * )) = κ(x * , x * )k ⊤ (K + σ 2 n I) -1 k (6.25)

6. 8 . 2

 82 Approximation of the Hilbert space normLet's consider functions f , f ∈ H and denote f -f as g. A vector g is an evaluation of the function g at the training points {x} N i=1 . Let's consider a functionh h = N ∑ i=1 α i κ(•, x i )with coefficients α α = g T K -1

||h|| 2 H

 2 = g T K -1 g and that h = argmin g || g|| 2 H s.t. g(x) = g(x) ∀x ∈ {x} N i=1(6.26)

g|| 2 H( 6 . 28 ) 1 λ I) - 1 K(K + 1 λ I) - 1 g(Figure 6 . 5 :

 2628111165 Figure 6.5: FBVI predictions of 20 particles

4 . 2 ,

 42 Variational Bootstrap allows obtaining faster convergence than MCMC. Also, this method could be extended to mini-batch training by means of Random Features similarly to the idea in Section 6.4.1.C O N C L U S I O N S7.1 summary of contributions and open questionsOptical Processing Hardware (Chapter 3)

Table 3 . 1

 31 LeNet-based binary encoder architecture 32

Table 3 . 2

 32 Image datasets used in the experiments 32

Table 3 . 3

 33 Classification error obtained by the model with the LeNet encoder 32

Table 3 . 4

 34 Architecture of the resnet-based binary encoder 34

Table 3 . 5

 35 Residual block structure. The number of filters is specified in Tab. 3.4 34

Table 3 . 6

 36 Classification error for models with the resnet encoder 35

Table 4 .

 4

	1	UCI datasets used for evaluation.	43
	Table 4.2	Classification error for variational bootstrap with the Dirich-
		let transformation and Markov-chain Monte Carlo approaches.	45
	Table 4.3	MNLL for variational bootstrap with the Dirichlet transfor-
		mation and Markov-chain Monte Carlo approaches.	45
	Table 4.4	Classification error of the proposed method and log-normal
		approximation	50
	Table 4.5	MNLL of the proposed method and log-normal approxi-
		mation	51
	Table 4.6	ECE of the proposed methods and log-normal approxima-
		tion	51
	Table 5.1	UCI datasets used for evaluation.	60
	Table 5.2	Comparison in terms of test set MSE between LSGPR a
		standard GPR.	61
	Table 6.1	Average RMSE for UCI regression datasets	72
	Table 6.2	MNLL for UCI regression datasets	72

Table 6 . 3

 63

  Bayesian model comparison and hyperparameter optimizationLet's now consider Bayes Theorem (2.24) from a different point of view. We can use this theorem to perform model comparison and selection. Given a set of possible models M l |l = 1, . . . , L and a dataset D, our goal is to define which model describes the data in the best way. The set of the models is discrete so that we can write a posterior probability of each model given the data D:

	2.2.2 p(M l |D) =	p(D|M l )p(M l ) ∑ L l=1 p(D|M l )p(M l )	.	(2.36)
	If it is necessary to choose only one model, we can find a maximum a posteriori
	model M MAP . And given the uniform prior over the models, the best model
	M MAP is given by
			We can sub-
	stitute the test point x * to these expressions and use the property of the sum of
	Gaussian random variables to obtain the expressions (2.28) and (2.29).

Table 3 . 1 :

 31 LeNet-based binary encoder architecture

		Layer	Dimensions
		Conv2D	5×5, 6 filters
		MaxPooling	2×2
		Conv2D	5×5, 16 filters
		MaxPooling	2×2
		Linear	576×512
		Linear	512 × d binary
		Binarization layer	d binary
		Table 3.2: Image datasets used in the experiments
	Dataset	Train/test size	C	D
	MNIST	60000/10000	10	1 × 28 × 28
	F-MNIST	60000/10000	10	1 × 28 × 28
	CIFAR10	50000/10000	10	3 × 32 × 32

Table 3 . 3 :

 33 Classification error obtained by the model with the LeNet encoder

	Dataset	AE	DFA	End-to-End SE
	MNIST	0.06±0.02	0.31±0.03	0.01±0.00
	F-MNIST	0.20±0.01	0.47±0.01	0.09±0.00
	CIFAR10	0.55±0.01	0.81±0.02	0.32±0.01

Table 3 . 4 :

 34 Architecture of the resnet-based binary encoder

	Layer	Dimensions
	Conv2D	3 × 3, 64 filters
	Conv2D	4 × 4, 64 filters
	Residual Block	3 × 3, 64 filters
	Conv2D	4 × 4, 64 filters
	Residual Block	3 × 3, 64 filters
	Conv2D	4 × 4, 128 filters
	Linear	4096 × d binary
	Binarization layer	d binary

Table 3 . 5 :

 35 Residual block structure. The number of filters is specified in Tab. 3.4

	Layer	Dimensions
	Conv2d	3 × 3
	Conv2d	3 × 3

image reconstruction task. It seems that the binary latent projection of image data, which is suitable for image reconstruction, is unsuitable for image classification.

Table 3 . 6 :

 36 Classification error for models with the resnet encoder

	Dataset	LBAE	End-to-End SE
	MNIST	0.15±0.01	0.01 ± 0.00
	F-MNIST	0.24±0.02	0.06±0.01
	CIFAR10	0.63±0.02	0.17±0.01

Table 4 . 1 :

 41 UCI datasets used for evaluation.

	Dataset	Classes	Training instances	Test instances	Dimensionality
	Magic	2	14020	5000	10
	HTRU2	2	12898	5000	8
	MiniBoo	2	120064	10000	50
	Drive	11	48509	10000	54
	Letter	26	15000	5000	16
	Mocap	5	68095	10000	37

Table 4 . 2 :

 42 Classification error for variational bootstrap with the Dirichlet transformation and Markov-chain Monte Carlo approaches.

	Dataset	VBoot-MLP	MCMC-MLP	VBoot-RFF	Sparse GP
	Magic	0.12±0.01	0.12±0.01	0.13±0.01	0.13±0.01
	HTRU2	0.02±0.01	0.02±0.00	0.02±0.00	0.02±0.01
	MiniBoo	0.11±0.01	0.09±0.01	0.08±0.01	0.08±0.01
	Drive	0.01±0.00	0.001±0.000	0.01±0.00 *	0.02±0.01
	Letter	0.05±0.01	0.03±0.01	0.05±0.01 *	0.08±0.01
	Mocap	0.005±0.000	0.007±0.000	0.02±0.00 *	0.03±0.01

Table 4 .

 4 3: MNLL for variational bootstrap with the Dirichlet transformation and Markovchain Monte Carlo approaches.

	Dataset	VBoot-MLP	MCMC-MLP	VBoot-RFF	Sparse GP
	Magic	0.31±0.00	0.29±0.00	0.33±0.01 *	0.35±0.01
	HTRU2	0.08±0.00	0.07±0.00	0.07±0.01	0.07±0.01
	MiniBoo	0.25±0.01	0.22±0.00	0.20±0.01 *	0.21±0.01
	Drive	0.06±0.00	0.01±0.00	0.08±0.03	0.08±0.01
	Letter	0.31±0.00	0.24±0.00	0.28±0.01	0.25±0.01 *
	Mocap	0.02±0.00	0.03±0.00	0.09±0.01 *	0.13±0.01
	4.3.8 Discussion of the approach

Table 4 . 4 :

 44 Classification error of the proposed method and log-normal approximation

	Dataset	Hybrid Transformation	Log-normal approximation
	EEG	0.18±0.01	0.21±0.03
	Magic	0.133±0.002	0.139±0.003
	HTRU2	0.0224±0.0006	0.022±0.001
	MiniBoo	0.083±0.003	0.085±0.003
	4.4.4 Discussion of the approach

Table 4 . 5 :

 45 MNLL of the proposed method and log-normal approximation

	Dataset	Hybrid Transformation	Log-normal approximation
	EEG	0.40±0.03	0.43±0.04
	Magic	0.336±0.006	0.351±0.005
	HTRU2	0.074±0.005	0.081±0.005
	MiniBoo	0.226±0.008	0.205±0.006

Table 4 . 6 :

 46 ECE of the proposed methods and log-normal approximation

	Dataset	Hybrid Transformation	Log-normal approximation
	EEG	0.027±0.005	0.022±0.007
	Magic	0.0263±0.0009	0.022±0.004
	HTRU2	0.043±0.001	0.039±0.001
	MiniBoo	0.0274±0.0003	0.0286±0.0009

empirical method of construction of the activation function proposed in this work is not applicable in a C-dimensional space.

Table 5 . 1 :

 51 UCI datasets used for evaluation.

	Dataset	Training instances	Dimensionality
	Yacht	308	6
	Boston	506	13
	Concrete	1030	8
	Kin8nm	8192	8
	Powerplant	9568	4
	Protein	45730	9

Table 5 . 2 :

 52 Comparison in terms of test set MSE between LSGPR a standard GPR.

	Dataset	LSGPR	LSGPR	GP	DeepGP	KNN
		Hilbert	Epanechnikov
	Yacht	0.63±0.12	2.02±0.58	1.09±0.05	0.93±0.13	57.80±16.65
	Boston	14.78±0.88 15.30±1.28	17.94±0.71 7.92±0.14	23.30±2.58
	Concrete	34.79±1.07 40.43±3.16	37.81±0.61 130.94±3.93 94.23±7.89
	Kin8nm	0.01±0.000 0.01±0.00	0.01±0.000 0.06±0.00	0.01±0.00
	Powerplant 14.65±0.34 14.40±0.67	16.85±0.69 14.57±0.15	15.35±0.49
	Protein	36.85±2.15 12.50±0.26	17.03±0.57 16.94 ±0.16 19.89±4.26

Table 6 . 1 :

 61 Average RMSE for UCI regression datasets

	Dataset	Functional	WDGPi-G	Deep
		MCMC		Ensembles
	Boston	2.73±0.02	2.83±0.92	3.69±1.15
	Concrete	4.06±0.12	4.80±0.41	3.07±0.26
	Energy	0.48±0.18	0.34±0.07	1.37±0.32
	Kin8nm	0.04±0.00	0.06±0.00	0.06±0.00
	Power	3.24±0.06	3.72±0.18	3.86±0.21
	Protein	3.61 ±0.04	3.65±0.02	4.45±0.02
	Wine	0.60±0.01	0.60±0.04	0.62±0.02

Table 6 .

 6

		2: MNLL for UCI regression datasets
	Dataset	Functional	WDGPi-G	Deep
		MCMC		Ensembles
	Boston	2.45±0.01	2.48±0.12	3.19±1.12
	Concrete	2.74±0.16	3.03±0.05	3.07±0.26
	Energy	0.80±0.05	0.35±0.15	2.07±0.98
	Kin8nm	-1.46±0.11	-1.23±0.01	-1.32±0.08
	Power	2.73±0.08	2.74±0.04	2.74±0.05
	Protein	2.73±0.01	2.75±0.00	2.80±0.01
	Wine	0.76±0.04	0.92±0.06	1.08±0.20

Table 6 . 3 :

 63 Average classification accuracy for UCI classification datasets

	Dataset	Functional	WDGPi-G	Deep
		MCMC		Ensembles
	EEG	92.51±1.82	94.13±1.96	89.04 ± 5.01
	HTRU2	98.10±0.26	98.03±0.24	98.03 ± 0.20
	Magic	88.16±0.33	88.37±0.29	87.90 ± 0.24
	Miniboo	92.54±0.21	92.74±0.39	91.49 ± 0.19
	Letter	98.22±0.18	96.90±0.29	96.38 ± 0.30
	Drive	99.45±0.09	99.69±0.04	99.33 ± 0.05
	Mocap	99.10±0.12	99.24±0.10	99.10 ± 0.08

Table 6 . 4 :

 64 Average test NLL for UCI classification datasets

	Dataset	Functional	WDGPi-G	Deep
		MCMC		Ensembles
	EEG	0.33±0.04	0.18±0.04	0.24 ± 0.10
	HTRU2	0.06±0.002	0.06±0.00	0.07 ± 0.01
	Magic	0.31±0.00	0.29±0.00	0.30 ± 0.01
	Miniboo	0.18±0.01	0.18±0.00	0.20 ± 0.01
	Letter	0.09±0.01	0.17±0.00	0.15 ± 0.01
	Drive	0.08±0.01	0.03±0.00	0.05 ± 0.01
	Mocap	0.19±0.00	0.03±0.00	0.04 ± 0.00

colophon This document was typeset using the typographical look-and-feel classicthesis developed by André Miede and Ivo Pletikosić. The style was inspired by Robert Bringhurst's seminal book on typography "The Elements of Typographic Style". Final Version as of April 25, 2023 (classicthesis v4.6).

the best Mean Squared Error (MSE) on a validation set. The LSGP model used on average 7 training points to make a prediction. We compared the predictions of LSGP with the predictions of standard GP regression.

As we can see from Fig. (5.3), the Gaussian Process with RBF kernel is unable to make reasonable predictions in the region where the target function contains highfrequency components. While some nonstationary covariance functions might be appropriate for this example, the combination of a standard stationary covariance function with the localization approach offers substantial modeling improvements.

UCI datasets

We evaluated the performance of the LSGP method on several problems from the UCI datasets collection and compared it against standard GPR, Deep GPs approximated with random features [START_REF] Cutajar | Random Feature Expansions for Deep Gaussian Processes[END_REF], and k-nearest neighbors (KNN) regression. In particular, we aim to compare the predictive performance offered by the localization against the baseline of exact GPR, and to verify that any performance gains are not just due to localization, meaning that we expect to outperform KNN.