A conceptual and practical introduction to Hilbert Space Gaussian Process (HSGP) approximation methods

Speaker 1 [00:05]

All right, thank you for being here. I'm very happy to be able to talk about this topic at PyData. So what I want to talk about today, it's a relatively new method to fit Gaussian processes kind of at scale, so to say. So it's a method that I've been using quite a lot in applications, and I know other PyMC developers are using as well, but we believe this hasn't kind of reached a wider audience because maybe it's too hard or it seems kind of still kind of new to try it out so my intention with this talk is to give you an overview of how and why this works and at the beginning might be a little bit confusing so to say but that's fine because the whole point is that we are preparing a lot of material to cover various aspects of this topic so this material would be about more about why does this work and kind of a glimpse on an example and we'll be preparing other type of documentation regarding kind of practicalities so yeah please bear with me and I hope you can find it as useful as we are doing it so the most important slide is this one because on this url which I can share on discord later on there's a fully and I would say lengthy notebook on which I go through all the details meaning code and also kind of the math so you don't need to kind of take notes I also share these slides you can also see the slides already if you if you go there and so just sit relax and take notes on what you don't understand so that you can then go to the notebook and and then also ask questions but yeah I prepare kind of really nice material that hopefully brings this into something tangible. I'm also kind of comparing, well, not comparing, but kind of providing various implementations. So in PyMC, which is a quite friendly, there's a quite friendly API for fitting such models, and also non-Pyro, and we do this a little bit by hand. All right. So I'm going to briefly talk about kind of the motivation and some applications, just go to the straight point and say we are actually using this in industry. I'm going to talk about Gaussian processes and I'm not going to spend a lot of time here so kind of this is the sacrifice of the trade off. I want to of course give you an example and provide kind of the problem the main kind of component of this talk would be about the approximation and why it works and we are not going to really go up to the last epsilon but really try to kind of hint the ingredients and the reason is because the APIs kind of that wrap this function has certain parameters and if you know what these parameters actually mean it might help you debug this type of of models and yeah i'll share some references if i i'm gonna kind of take it easy and if i don't cover the whole the material that's fine because you have everything online and you can always reach to me so about applications i mean when i took when people talk about gaussian processes they usually talk about the toy example that i'm going to present today but in reality it can it can serve many purposes so the first kind of case study that i where i saw gaussian processes being a fundamental part of a of an interesting real model was this blog post by alexandora and he was trying to model the popularity of presidents in france so you see kind of the the y-axis is between zero and one this is a share of popularity and he built a very nice model where he was using kind of these gaussian processes as kind of latent variables for this share kind of a variable and also doing contrafactors so for example okay what if for example the unemployment rates would drop to five percent what would that actually mean in terms of popularity so i sent i have the link down there please check it out it's super nice and kind of this was the first kind of study which actually hooked me into trying to work on Gaussian processes for real applications. Most of the times when I work with Gaussian processes, it's actually to encode time-varying coefficients in linear regression. So in the usual linear regression, you have, let's say, y as a function of x, and you have a regression coefficient, they call me beta, which is fixed. But in many applications, if you allow this actually vary smoothly over time, then you you can grasp dynamics which are relevant for real case scenarios. So this is taken from a little case study that I wrote about this classical bikes data set where you want to model the number of bikes rented in a city as a function of many covariates on temperature. And if you fit a linear model, you have like a coefficient. Let's say how many bikes do you expect to be rented if you increase the temperature degree by one. But if you actually do like a time variant coefficient, you see that after 25 degrees or so, kind of the effect is actually lower than what you would get from a linear model. And actually, it's very interesting because I also compare it against having this type of model, this type of data fit into kind of an XGBoost, and if you do kind of this partial dependency plot, you will get the same phenomena. But I just find having control on this kind of time-varying coefficients, it's very nice. where am I using it and we are using it a lot well we work on trying to kind of capture marketing efficiency over time but it's called kind of media mix models and you would expect in general kind of the the efficiency of certain channel it's a the spend you put a to to have an effect on conversion so to say and you're interested in the cost per acquisition and actually I believe that the kind of base models where you have kind of regression coefficients which are not changing over time are a good baseline, but if you want to kind of try to understand the dynamics over time, then this comes very handy. So I link also a blog post by the folks from PyMC Labs about this. And you can also vary the coefficients, but you can also vary the intercept. And by mixing these two, you can actually get models which fit the data pretty well. So I hope this conveys the message that what I'm going to talk today, even though it might look a little bit weird, it actually can be very helpful. So let me go into the section of Gaussian processes. So, oh no. Okay, let me try classic. One second. Otherwise, it's going to be very weird. Okay, so I hope you can compile LaTeX. This is the first time it happens, I apologize, and it's going to look very odd, so apologies for the recording as well. Okay, let's do that. Thank you very much. You can still get the raw data because I provide that as well. So typical use case for a Gaussian process is that you have kind of the blue data, which is your training data, which looks weird. I know kind of where this is coming from, this is just a combination of sine and cosine, but actually in real practice I don't know how this is generated, and I want to generate predictions kind of out of sample data, let's say on orange dots. And this is what we want to do. We can do this in different ways. We can use splines, we can use regressions, but we're going to see how this works for Gaussian processes. So the key ingredient, and I think this is why I like Gaussian processes a lot, and this is a very opinionated way of thinking, is that you can pass a lot of information through what is called the kernel. And you, as a modeler, are telling the model, let's say, kind of a measure of similarity between points. because if you see this plot, I mean, I show in black the real data generation process which is very smooth but you could also tell the model like I'm expecting the Gaussian process to be not that smooth or to be periodic or to have certain properties so this, you as a modeler, pass as a kernel so this is not kind of a black magic that you're doing so a kernel is a way of encoding similarity between points and the most used kind of kernel It's called the square exponential, and this kernel takes two elements and gives you a number. In this case, I'm just taking kind of this exponential of the difference and multiplying by kind of certain parameters. And these parameters are important. The A parameter is the amplitude, and it's going to control, let's say, how this shifts on the noise, so to say. and the l is the length scale which kind of it's really actually the most important parameter because it's capturing how let's say how i'm measuring the closeness so to say so if l is very big as compared to the range of x then if i change x prime a little bit it's going to make very little difference whereas if l is small as compared to the range of these values a little small change in x prime will actually change the value of the kernel quite a lot and yeah there are many other type of kernels and if you have what is nice is that these gaussian processes can work as lego breaks because you can actually have many kernels you can add them together and multiplying them and kind of it's a bit it requires some practice but if you you can practice this intuition to actually capture the kind of what you expect as structural components of your data within these kernels. You, of course, don't know what A and L are, and you would like to infer them. This particular kernel has a property that let's say the dependency on x and x prime is through the difference, which means that, and this is important because it's going to come later, this is what is called a stationary kernel. So it really depends about the relative difference between these two. So if you plot this as a matrix, you can now evaluate this, let's say, with respect to the whole training data set, and you have a matrix that is, the dimension is kind of the size of the training data set, and you see that kind of the kernel is concentrated along the diagonal, meaning that kind of points that are close together are going to be more similar. This is obvious, but this is what the kernel is looking. But notice something important, is that this matrix actually scales with this kind of, it's as big as the training data set, which means if the training set is very big, this matrix is going to be quite big. So what's the whole point about a Gaussian process model? Well, you are going to assume that you can model this data as a multivariate normal. So this is kind of a generalization of a normal distribution where you have certain mean, which is not very important, but you care actually about the covariance of this matrix, which is kind of encoding the similarity between these two points. So I denote x and y, kind of the training data set, and by x star, the test set. I would like to know what y star is. I don't know that. But kind of the way the Gaussian process works is that you are going to assume that the observations y and kind of the latent Gaussian process f, because kind of you have this line that fits the data, but also you have the noise. So I'm assuming that the observations and this latent kind of fit on the test set come from a multinormal distribution. And this matrix kind of bold K is actually splitted in between kind of the components. So it's evaluating the kernel matrix on kind of here the training set, here the training and the test, here the other way around, and here the test set. And this little sigma here is just to account for noise. and that's kind of your your your parametrization and what you do and you can actually do it by hand is that this fit is just take this multinormal a normal multinormal distribution and condition on the training data which is x and y and on the points of the test data where you want to generate predictions so it turns out that if you do this game you can all you will also get a multiviral distribution, multinormal distribution, with certain mean and certain covariance. It's not important. What is important is that there are explicit formulas that you can compute by hand, but there's an interesting thing, is that if you want to do it this way, you have to take the inverse of this kernel function, of this kernel matrix. So, that means this is going to require a lot of computation, and you can do kind of Cholesky the composition if you want but there is no easy way of getting out of the fact that doing this for large kind of data sets is actually quite expensive so the the takeaway of this is that Gaussian processors are nice but there's a price to pay and it's through this kind of matrix inversion formula but nevertheless we can do it for this specific example and in the Bayesian kind of paradigm you want to understand also the the the values of of L and A and we can impose prior distributions and this is a huge topic by itself about how to set priors on GPs on Gaussian processes but the intuition says that I will expect the length scale prior to be less than the amplitude because kind of if you see the plot before kind of the range between 0 and 1 and I have a lot of data points so I want this distance to be kind of sensible respect to this distance and the amplitude is between minus one and two I don't want to go into details but this is kind of an ingredient to the Bayesian inference framework that you want to do and you can do this very easily in PyMC so I also do it in non-pyro but kind of what you tell PyMC if you are not familiar with the PyMC syntax just think about like input output So you pass the data, you pass the priors of your kernel parameters and also about the noise. And there's a very convenient kind of mean function, which we're going to assume is zero, and a covariance function on which you pass the kernel, and you pass the amplitude here, and your length scale, and we just fit that. and you are going to assume, okay, the likelihood, it's a normal likelihood. And this F that I have here is the latent Gaussian process. So 50 light of codes, and you can fit this. Yeah, this is kind of how the model looks. The input are the priors. The output is the data. And within the Bayesian kind of framework, before kind of looking into the data, you can run what is called these prior predictions. I just mentioned that this is important if you're doing Bayesian modeling. to assess whether the priors are good and here I'm just kind of generating kind of feasible curves from the model without conditioning on the data so this is on the first step and we see that we're not being very restrictive and it's plausible and then you feed your model and what you get out of this feed is distribution of the kernel amplitude and length scale and the noise and at the very end And PyMC gives you a very nice method on which you pass the test set, and then you condition. And this is doing this huge matrix thingy. So you actually don't need to know much about how this is happening. But if you fit this, this is going to take maybe four minutes, even though we have kind of 80 or 100 data points. So it's not great if you want to scale it to like 1,000 data points. But you get this nice fit on which the orange plot is kind of the likelihood and the pink one is the Latin-Gaussian process. So within the range of definition, we capture this quite well. But out of the set of definition, actually the values that are kind of, the predictions are kind of converting back to the mean. And also you get a lot of uncertainty because the Gaussian process says, like, I haven't seen that much data outside, so I want to be conservative and just go to the prior again. So this is just a little bit of caution. Let's say extrapolation with Gaussian processes is something that has to be done carefully. It can be done, but kind of depends on the data and the kernel. And then I can do this for time-varying coefficients. I can do all of these predictions on the presidential kind of popularity, and that's fine. But the key thing is that this doesn't scale well. And what we want to kind of present here is an alternative that has proven to be very successful. So this, I think, was published in 2020. There have been a lot of approaches, but this actually works quite well. So again, the problem is that if I want to do something with classical GPs, then I need to invert this matrix, and this is expensive. So the approach of the summary it's the following we want to find a way of approximating this kernel matrix with something which has a smaller rank because if we do so and if the approximation is nice then inverting something smaller it's easier so that's what we want to do and here kind of i need to also have a disclaimer i i like kind of the ideas behind it so because it connects various fields so this is going to sound a bit weird, but we'll try to make this more tangible. The whole point is that, and this is why we're talking about Hilbert spaces, is that there's a way of interpreting the kernel as kind of linking that to a certain differential operator. And differential operators on kind of functions form also kind of linear spaces, which are called Hilbert spaces, and you have a lot of spectral theory, like eigenvalues and eigenfunctions. So we're going to borrow a lot of tooling from that world to try to solve this problem. It's important not to get lost on the details but to focus on the objective. But also, what I like about this is that it brings a lot of this spectral theory into the game. And at the very end of the game, we're going to essentially reduce fitting a Gaussian process to a linear regression. So hopefully that gets you excited. Well, not as excited as getting matjacks to work, but close. So this is our diagram. I'm going to be very emphatic about the strategy and kind of try to go fast through the details. So we are starting here from the kernel. And what we're going to do is to associate what is called a spectral density, which is just a nice representation of a kernel. and then we're going to try to do kind of a Taylor expansion to see what we get. We get nothing unless we do something. And this other something is actually bringing what is called the Laplace operator. And the Laplace operator is being studied in physics, for example, quite a lot, where it appears on the wave equation, on the heat equation, so it's quite well studied, so we have a lot of tooling. and we're going to link this through the Fourier transform and at the very end we're going to arrive to an approximation formula. So this is kind of what you need to get in mind when you think about these APIs and filling the data. I want to give you a glimpse on why this works just because I think having a reference is useful but at the end of the day you can simply use the approximation formula. So what are we going to do? Going back to the diagram. So we're going to give a detour on eigenvalues and eigenvectors because I think this is kind of the key component of this and kind of it's often kind of overlooked and I think it's a very strong theory. Then as I mentioned we're going to associate to each kernel what is called a spectral density and we're going to run a polynomial expansion on this thing and I haven't done anything unless I bring the Laplace operator and I'm going to show or explain that if you take the Fourier transform of the Laplacian you can actually plug it in into this approximation but you need to set certain boundary conditions to get all of this spectral theory working and then we will do that in a second and then you put this together just by kind of eyeballing the formula and you get to a final linear regression and this linear regression is what we have down there so it's not really important about kind of what's happening the thing is i'm going to be able to fit a very simple linear model because i have like a like a usual beta coefficients here i'm going to have this what are the eigen functions of the laplacian which is something that is easy to compute so they're in the books they're a function i don't care about that and what is most important is that they don't depend on the hyperparameters of the kernel which means like a and l so i just can store them somewhere and reuse them as i wish and the only thing i need to do is to compute the spectral density which is a function again that i know there's no nothing i need to do there on the eigenvalues of the laplacian which is something that i know so the story is that i do know how to do this very easily with a linear regression So what is important for the detail is to count the factor and understand why here's a root square and so on. But at the very end, it's a linear regression. So thinking about applications, what in PyMC run kind of the simple model in four minutes with 100 data points, it runs in milliseconds with the approximation. And there's already kind of a nice API to do this. So I guess the main component here is about eigenvalues and eigenvectors, so you have a square matrix, remember that an eigenvector and eigenvalue pair satisfy this property, meaning if I take a matrix, multiply it with a vector, I get a scalar lambda times v. I want v not to be zero, so that it's not trivial, and I say that the spectrum of this matrix is the state of eigenvalues and eigenvectors. And probably you've done this quite a lot in school, but it's still very, very useful. So, if you want to see an example, I have these two by two metrics, and one can show that the eigenvalues are three and minus one, with eigenvectors one, one, and minus one, one, and one can add this kind of funky factor is so that kind of the inner product between these two is zero and the norm is one. And there's kind of an easy way of doing this. I also provide resources if you are a little bit rusty, but it's kind of the gist. But kind of eigenvalues and eigenvectors by themselves might not seem that interesting. What I think is one of the most interesting theorems or resulting linear algebra it's what is called the spectral theorem and it just simply says that if you do that procedure and you use these eigenvectors as a change of basis so trying to change the coordinates your operator or your matrix which in this if it's symmetric is always going to become diagonal which means you can just work with diagonal matrices which are just much nicer So this is a very interesting result, and it can generalize to a much wider kind of scope. So I probably don't want to go into details, but I provide this example, so you can use JAX or NumPy, and you can actually test that what I'm telling you is true, which means if I give you a matrix, minus 1, 2, 2, 1, you can compute the eigenvalues, and the JAX gives you this in a normalized way which is 1 over the square root of minus 2 of 2 and if you do the change of basis kind of game you actually get a diagonal matrix which is super nice and one interesting thing here is that you can actually write your operator as a sum of the eigenvalues times kind of taking the inner product with V and multiplying, which if you apply this to a vector it will give you a scalar and then multiply it with the vector again. And this is kind of, it's very, very useful for what we're going to do because instead of working with matrices we can just work with these sums. And these sums are actually quite convenient. So this is what is called a spectral projection because it's indeed a projection. If you do this twice for these operators, you'll get the same thing. And in the notebook, I show you some examples where I see if I kind of do this, this actually works for all vectors. But the most important thing about this is that if you have such a representation of an operator A, if you take a function of this matrix, let's say take a square or take an exponential, This is actually as simple as taking the function on the eigenvalues here. So remember that multiplying a matrix is like you need to do this game with the rows and columns, but here I'm telling you just write the matrix as this kind of spectral decomposition, so to say, and if you want to apply any function, then the only thing you need to do is apply the function to the eigenvector. So easy examples, if you have the identity function, you recover the operator a if you put f z equals to one then this is just the going to be the identity matrix and you can take the square just by taking kind of squares of the eigenvalues and you take and take exponentials of matrices just by taking the exponentials of the eigenvalues so this is what we're going to do and this is why the spectral theorem is so useful is because it allows you to take any type of functions on matrices just by evaluating this so you can do this by pen and paper you don't need to mess up with any matrix and it's just very convenient so this was kind of a small detour but what we want to do is actually go to our kind of kernel so one thing that I that one can show is that if I have a stationary kernel which means it depends on on x and x prime as a difference then you can always write it kind of in a as an integral representation this is the ugliest slide that i have so apologies for that but think about it as a fourier transform it just happens to to work and there are very kind of clever mathematicians that actually compute this so given a kernel you can go to rasselman's book and just look for the spectral density and they say okay i have a formula and they look like this so they depend of course on the parameters of your kernel so this is for the square exponential and this This depends on actually a variable omega, and this goes as a kind of Gauss-Bell shape in omega. But what is important is that you can just write these formulas in NumPy, and I provide such a code. So it's very simple. So again, think about this as a Fourier transform, and where you don't need to think because these quantities are given. And here I show how the spectral densities for the co-exponential look. have this feature that they always decay at infinity. So what we can do, actually, is assume that there's going to be a Taylor expansion of these spectral densities. And you can just assume this with certain parameters, A, which you don't know. And you say, OK, if I were able to read this from somewhere, then I could maybe truncate this series and do something clever. But at this point, this is where other side of the story comes in and it's the Laplacian. So the nice thing and I guess the key idea of the approach is that the expansion they do on the square of omega. I'm going to give you a reason for that that the paper doesn't give because they go really fast and it has to do with the Laplace operator. So in another world, you have the Laplacian which physicists have studied quite a lot and this is this operator that takes the second derivative of a function and it has a lot of properties that people like but one property that we're going to use is that if you take the Fourier transform of this operator it gives you a polynomial of order 2 so the whole idea is, ok, I have this in my expansion this actually can be interpreted as the Fourier transform of the Laplacian so let's try to plug this in by brute force and see what happens so you can actually do that so you just can consider the operator on which you replace the omega square by the Laplacian and see what can happen because that's another way of interpreting this omega square but here you haven't gotten anything because if you have studied kind of the differential operators you know it's important to define a domain of definition or boundary conditions. So one very common boundary condition is to assume that the functions in the domain of this operator vanish in a boundary of certain domain. So you can put your data into certain box and assume that all of the functions are going to vanish here And if you consider the Hilbert space, so that's why the name of the approach, of all of the functions that are square integrable and satisfy the boundary condition that they vanish on the boundary, then this is actually a symmetric operator or self-adjoint. And what you gain with this is that you have the spectral theory. You can do exactly what I said for matrices, meaning the spectral theorem, taking functional calculus and taking powers of the operator just by taking the eigenfunctions of the Laplacian. So you get a lot of kind of spectral theory, and we are very happy because we know how to actually take powers of this polynomial. That's the motivation of kind of the bottom formula. We want to do this in a clever way, and we know if we have a spectral decomposition, we can do it. just very easily by taking powers of the eigenvalues. So this operator on the set of functions that vanish on the boundary actually has a discrete spectrum which grows to infinity and I have eigenfunctions which actually I can compute. And why I can compute? Just because you know how to solve this differential equation which is just take which function if you take the second derivative would give you minus that function and vanish at the boundary. Well, this is the sine functions. So I guess if you have taken any physics course or basic math or kind of differential equations, you know that this is kind of the first exercise that you can do. And so we all know how to do that. You can wrap this into a little function and just move forward. So if you want to see the plots, these are just sine and cosines. sorry, just signs moving around the boundary and usually take the box, let's say this minus LL to be sufficiently big so that you have the training data here because that's going to be where you want to do kind of extrapolation. And kind of the final step is actually noticing that if I want to take powers of the Laplacian, the only thing I need to do is take powers of kind of the spectral decomposition and this is what i've written in red if you can see so this is just the last step it's about kind of identifying the Taylor expansion of the spectral density and the and computing this in terms of the of the eigenfunctions of the Laplacian but notice that you don't need to do anything because the if you want to take powers of this huge thing you just need to take powers of the little eigenvalue. This is what you're winning. You don't need to compute anything if you have such spectral representation. And through these identifications, you realize that in order to kind of identify this coefficient, kind of the parameters are actually evaluating the spectral density on the eigenvalues of the Laplacian. This is just by eyeballing it. So if you look into the papers, this takes kind of less than one page but they don't explain the why so kind of my my take-home message is that it everything has to do with the fourier transform this is kind of the link and at the very end you get this final representation which kind of if you translate it into the world of the gaussian processes it just says that you can truncate you just compute the eigenvalues to the Laplacian compute the eigenfunctions and compute these things that you can actually do very easily because you can pre-compute it, they don't depend on the Gaussian process, this blue one, and this one has formulas you can get in the book and we get to the final linear regression so at the very end we get to this very easy formula which again you don't need to know it or the derivation but I think it's helpful when you're working with these models and for one dimensional case I do this in the notebook I do this explicitly and actually show that what I do in non-Pyro by hand actually coincides to what PyMC is doing behind the scene so again, take home it's this, I can share a link I will share the link for the papers but it's just another way of seeing how all of this spectral theory can help us solve a very concrete problem so back to the example these are now priors that I can get but instead of having a GP I have this linear regression so this looks quite reasonable everything I told you about this weird theory can be written in less than 30 line of codes because these are just kind of sine, cosine, exponential and it's provided there So this is something that is just kind of a linear regression. So you can just use these helper functions for your applications. You don't need to actually rethink it yourself. I want to skip this because I want to show you that the results that you get is the traces on the parameters, which are the kernel amplitude, the length scale, and the noise. These are actually the same, like almost the same as what you get with the Gaussian processes, kind of vanilla implementation but you also have all of these beta coefficients which is kind of the number of terms in the expansion which are just kind of linear terms and this fits very very fast and this is the result that we get which is very similar to what we got before it's pretty much the same PyMC actually gives you a very nice way of doing it and the only thing you need to change is GP to SGP, or like HSGP. And it provides kind of a very simple API, so going from vanilla Gaussian process to Hilbert spread approximations, you just need to change one line of code. But now you have certain control of the number of functions to use. In this case, I'm using 20. And the L, which is kind of the size of the box. So hopefully, if you use PyMC, you already know that this L has to do with the Laplacian where I'm doing the approximation and this M has to do with the eigenvalues, eigenfunctions that I take. And I can use the conditional method. I have to do the same thing and I get exactly the same results. So you can take this as for granted. PyMC gives you a way of actually getting the spectral density yourself and the eigenvalues of the Laplacian, eigenfunction of the Laplacian if you want to do something funky it's also supported so please check it out again i do both implementations in the notebook so a kind of i can do out of sample predictions in a similar way yeah and that's what i wanted to talk about and about references as i mentioned kind of we have the classical papers here on the regression process but this little notebook is actually an open PR where we are working on having even more practical examples on how you can use hierarchical models on these Hilbert space approximations and do this on, for example, two dimensions, which Hilbert spaces are also quite useful. Yeah, so I'm sharing this with you, and yeah, thank you very much.

Speaker 2 [39:46]

Thank you, Juan, for your introduction to HSGP. We have some questions on Slido. The first one is, how does HSGP compare to other GP approximations, especially to stochastic GPs fitted with variational inference?

Speaker 1 [40:04]

So, I don't know the details, I haven't studied that much those, but I do see kind of a difference from what I heard, is that here you can use full MCMC, whereas in variational inference you need to kind of specify kind of a family of the posterior distribution. So, I guess what we gain here is that you don't need to impose that, and it could work well if you have kind of bimodal distributions, but I don't have any concrete answer for that.

Speaker 2 [40:40]

Are your choice of boundary conditions causing bias, especially for large Lange scales?

Speaker 1 [40:48]

So kind of what it's recommended, this is kind of the theory, but what it's recommended in practice is that you center the data. So if you have your data like this, you actually move it in such a way that the mean is zero and you take L big enough. And this kind of, the choice of L, it's a little bit tricky because if it's too small, it won't be able to extrapolate that well, but if it's too big, then these things are going to vary quite slow. So, there's a paper here, which is this one, which I strongly encourage you to read, which is Practical Hebrew Spell Approximation, and this is where the authors go into the empirical studies about how to select this, and this is more kind of trial and error, and kind of rule of thumb tricks.

Speaker 2 [41:38]

So, short HSGP is much faster than vanilla GP? Any drawbacks of HSGP?

Speaker 1 [41:45]

So it's definitely much faster, and I don't think you lose a lot. Of course, there are things that don't work. For example, if you go to very high dimensions, then it doesn't work that well. So for dimension 1 and 2, which is the cases where I've used this a lot, it works pretty well, but it also has the limitation of the input dimension.

Speaker 2 [42:09]

For which kernels does HSGP work?

Speaker 1 [42:12]

So it works for most of what you can think, so square exponential, all of the Martin families, which are the ones which are non-smooth, and you have all the spectral densities, so you can just use that from the books. A periodic, it's a little bit weird because it doesn't work exactly like this, but the authors of this paper actually make it work. So I say in practice, all of them that you will use in a typical use case. So for the periodic ones, you can. Yeah, no, but for example, the linear, then it's faster to fit. So yeah, I think the linear won't, but kind of the ones that you usually use, yeah. But yeah, I agree that some of them won't.

Speaker 2 [42:59]

What again is the speed up of HSGP in your example? Sorry? What again is the speed up of HSGP in your example?

Speaker 1 [43:06]

in your examples? So I didn't benchmark quite thoroughly, but at least in PIMC, it went from four to five minutes to milliseconds. And in practice, I usually never work with vanilla GPs.

Speaker 2 [43:18]

If the main advantage of HSGP is speed, can we apply it to the deep Gaussian process?

Speaker 1 [43:25]

I'm not an expert on that, but I would if you are using vanilla caution processes I wouldn't imagine why you shouldn't or you couldn't do this, but I haven't tried myself

Speaker 2 [43:35]

Last question. What about high-dimensional data?

Speaker 1 [43:38]

Yeah, that's what I mentioned, that dimension one or two works fine, higher dimension gets more complex because now, and this is why it's important to know about the theory, you'll have the eigenfunction of the Laplacian on all of the dimensions, and then you need products of them. So this will increase quite a lot with the dimensions. But in most cases that I've worked, which is dimension one and two for geospatial data, it works quite well.

Speaker 2 [44:03]

Thank you again for answering the questions. All right. Thank you for your talk.

Speaker 1 [44:07]

Thank you.