Gaussian Process for Time Series Analysis

Speaker 1 [00:04]

Thank you very much. So I would actually like to thank Vincent because he did a great introduction for the topic that I'm going to talk about, although I think we're going to play good cop, bad cop because I want to go a step further and try to understand the Gaussian process in terms of the Bayesian approach, so to say. So that's the idea of the talk, and at the very end I want to focus on how you can use this in a concrete use case. So I'm going to start with a Bayesian linear regression. Then I will go to the definition of a Gaussian process as it is and try to understand and make sense of it. And at the very end, I will, let's say, give a little bit of detail on the kernels and how you can use them to do forecasting. So I guess we already have seen what the normal distribution is. So essentially, whenever you have a random vector, it says it has a normal distribution if the density looks like this. So essentially you have a mean and a covariant functions which, let's say, specify the location and the shape, the support of the multinormal distribution. So here we have two examples. On the left-hand side we have seen how this nice ellipse shows how a normal multinormal distribution looks like. It can happen that each of the components is normal And this weird support actually is an example of something that is not multivariate. So, you can see that the devil is in the details, so I want to, let's say, point out some of these things. So, let us start with the basic linear regression. So, what I do have here is around 100 data points, which I generated with random noise. And, let's say, the model has intercept 1 and a slope 3. And I just want to know what's the best fit for that. So what I want to do is get the intercept and the slope. But actually, I can do more or aim for more just by trying to find the distribution of the parameters in the Bayesian point of view. And how do I do that? Well, I do that through the Bayes theorem, which says that if I have this set of data and I model it with random noise, then we can find the parameter distribution as posterior, which is proportional to the likelihood, times the prior. That's the Bayesian game. So for the prior of the parameters, I'm going to select a multinormal distribution, which in this case has dimension 2, so the intercept and the slope. and before, let's say, seeing the data, I'm going to sample from the prior distribution to see which are the possibilities that I have. So the game is as follows. I have this specific matrix 2, 1, 1, 2, which I decided to take because that's my prior knowledge and I just sample from that. So I took one sample from the multinormal distribution that gives me two values, the intercept and the slope, and I draw a line. And I do this many times by sampling from this multinormal distribution, and what I get are all of these lines. So these are all the possibilities that I can actually get from the specification of this prior. But now, this is just from my prior knowledge. I actually want to see how the data might affect that, and this is going to be encoded in what's called the likelihood function, which is essentially the probability of the data given the input x and b. And it turns out that this is a multinormal distribution again. So in Bayesian analysis, getting the posterior distribution is sometimes not easy. And instead of doing the math, you actually sample from it. So you can run a Monte Carlo simulation to try to estimate the posterior distribution of these parameters. So what we see here is that the estimation with the sampling approach of beta 0 is around 1. This is OK-ish, and also we have the standard deviation saying, OK, how certain we are of getting those values. And we have something similar for the slope. But actually, the multinormal distribution, the Gaussian distribution, has such nice properties, which it's really strange, but a lot of the results can be obtained analytically. so instead of sampling if you do a little bit of math just a little bit then you can show that the posterior distribution of this parameter is again a Gaussian and I guess maybe what it's important is that the mean is somehow proportional to this product and maybe you can see that this is essentially capturing the OLS solution of the linear regression so essentially it's hidden in the mean And now that I have this analytical solution, well, actually, I could do this by sampling just from this data. I can just play the same game. I sample from the normal distribution, take two values. That's going to give me an intercept and a slope. And I play the same game, and I plot this in green. And we see how all of these lines actually fit the data much, much better. So this is how we usually work the linear regression problem in the Bayesian setting. So if I want to do inference, if I want to do a prediction, then what I need to do is weight the probability of, let's say, this new input point times the probability of the parameters and take the sum, or take the sum of that. And again, it turns out that just because of the form of the multinormal distribution, this is again a multinormal distribution. And if I want to do a prediction, I will sample from that, and I was unlucky enough that actually the prediction point went quite well to the mean, but if it's here within the credible interval, it's fine, because let's say that's the uncertainty that I have in my model. And this is something that says that if you want to do something that is a polynomial regression or are more nonlinear terms in an additive way, you can do that and you can do the same thing. So just you need to replace, for example, if you want to take a polynomial, instead of having x, you will map this to x and the second component x squared and the third component x cubed. And you can actually play the same game and you'll get exactly the same results. and independently of the math you see that at the very end you get some explicit expressions which seem quite mysterious but the only important thing and this is the main message is that everything in this problem can be obtained as the product of phi which is the embedding or x sigma, the covariant matrix and another phi these are the building blocks so you can glue this thing together to get the predictive distribution of your linear model. So essentially, if you have that, you have the solution of your problem. So why not define this as a thing by itself? So we define the covariance function or kernel in this specific case as a function that takes two points and it just applies the phi to x or you can just think of it as a vector as x in the linear model multiplied by the covariance matrix and then multiplied by phi of the second point. So this is a real number 3, 4, whatever and this is the meaning of the kernel. And it's important because it specifies the solution of your linear model. and this little phi here it's what is called the kernel trick where I can just embed this let's say polynomials the x into a polynomial space but at the very end you don't care about it itself, you just care about the kernel right? so yeah, so when I also first read the definition of regression process I was super confused and the reason I'm giving this talk is to try to see how much I can understand this And I did want to match the definition that I found in Wikipedia with some experiments that I'm going to show you in a second. So, I guess the moral is that whenever I give you or someone gives you a kernel function, this actually implies a distribution of a function. So, at the moment, this makes not a lot of sense, but we're trying to make sense out of it. So, if you go to Wikipedia, this is the notion of Gaussian process. So a Gaussian process is a collection of random variables, any of which a finite number of them have a multinormal distribution. So what? And it says that you can have the mean and the covariance is what specifies the Gaussian process properties, which is the usual mean and the covariance in the usual case using the expectation value. But I want to just mention that we already have seen a Gaussian process and the Gaussian process is just the linear regression with this prior so let us see why if I take a prior which is a normal with mean 0 and sigma as the covariant matrix and I take this product with this phi or even the vector this is a linear combination of normal distribution which by definition since this thing is multinormal is it defines a joint multinormal probability. So this is a Gaussian process in this definition, and it has mean zero and the covariance given by Wikipedia's definition actually coincides with the natural brick which we define via the simple linear regression. So this is an example of a Gaussian process and now I want to show why people often say that the Gaussian process is the infinite dimensional analog of the multinormal distribution in this function spaces. But instead of doing the math, I'm going to, let's say, just explain you how to sample from a function. So to sample from a function is really hard because these things are infinite dimensional, and you don't know how to do that. But you can specify a function just by saying what's the image of any grid of points. so in order to sample for functions what I do is first I set a grid of points which I call x star and then I do the following thing I define this matrix capital K which has each kernel combination for each point so in this example I have 80 points x star and this is a matrix of dimension K times K so whenever I sample from this I will get 80 points and I draw them here and this would be one of these iterations so it turns out that when I sample from this in this way the functions that I get are linear functions but that's not clear why because I'm just sampling from this multinomial distribution with this kernel So what is important is whenever I give you a different kernel or the kernel changes, this actually is going to give you different properties of the sampling function that you will get. Yeah, and we haven't seen the data, so if we have, let's say, the same linear regression problem, the next thing that we want to use Gaussian process for is to try to find a fit. But this is, I guess, the key slide on which, instead of sampling from parameters, so here we don't have, like, intercept or slope. We're just sampling from this multinormal distribution, and it happens in this case that we get a straight line. But given the prior, if I just construct the same, let's say, Gaussian process with mean zero and the kernel matrix, well, essentially the same, but just including the training data, I can now say how I would find the best fit well I will sample from this one and then I will take out all of the lines which don't fit my training data so like if they're really far away, take them out if they're good, then just take them but we don't need to do that because we're actually what we want to do is to condition on the training set so instead of just let's say staying with this pure Gaussian process we can condition on the training data and the input data. And this thing is, again, a multinormal distribution. And what is important is not the formula, but it just depends on the kernel. So the kernel is capturing absolutely everything here. And if I do the same game where I sample, again, I'm not sampling on parameters, but I'm sampling on functions, this is what we get. So, if you, let's say, work out these formulas, these are exactly the same expressions as the Bayesian linear regression. They are exactly the same. The only thing that changed was the point of view on which I forget about the parameters and I start sampling in this Bayesian way in a space of functions. And that, for me, is when I began understanding what a Gaussian process was actually doing. But then you have a new universe open up for you because, as we've seen in the previous talk, the selection of the kernel gives you different types of curves for your Gaussian process. So, the first one is, so there's a zoo of kernel functions, and these are ones which have somehow an interpretation. So the first one is just the one that we have seen for the linear model, which is called the dot product kernel. And actually, you can take powers of it to generate high-order polynomials, as we've seen before. Then we have the square exponential, which is just taking the exponential of the differences. And actually, one can prove the following. If you take your linear regression with these polynomial terms, like x, x squared, x cubed, and you take all of the polynomials, the kernel that you will get from this linear regression in the limit is this one. So the square exponential can be thought as the natural kernel coming from a linear regression with infinite polynomials. And that's why all of these lines drawn from this Gaussian process are really smooth. Because a natural question is like how smooth these lines are going to be because we're seeing the interaction between the points depend that all of these kernels depend essentially on the difference. It doesn't need to be like this. But if you want for whatever application control the smoothness, then you want to go for something which is maybe a rational quadratic kernel. And this has another intuition or like a way of, let's say, thinking about it. So it actually, it's a weighted sum of square exponentials, right? And as we have seen, you can also model periodic signals by introducing periodic functions like the sine function. So let us try to see an example on which this can be applied to the nonlinear setting. So it's not important about the function that I'm taking. It's a sum of sine functions. But if you are giving this, trying to fit a linear model by yourself, trying to get the features, like is it cubic or polynomial order 4? Well, actually, this is where Gaussian processes can come really handy. So what I'm doing here is just generating random noise out of this nonlinear curve. And I want to see how I could fit the points. So it's essentially the same game. I have my prior, and the kernel that I'm going to use is the square exponential. So there's a little bit of code you can unfortunately not see, but I will share this, and I promise you can actually run these notebooks. And what is important is that I'm running a sample from the distribution. Here I have 80 points where I want to test. And then I do the selection, and I plot the images. And you see that in this Gaussian process, we get a lot of these complex and curly curves. And the fact that there's some curvature and complexity is just because this is a different kernel. And different kernels give different natures on the curve space. So these are all your possibilities. So now let's see how, let's say, the joint distribution looks like. And as we have seen, the kernel just depends on the differences. And that's why, if I plot a heat map of the covariance of the joint distribution, this is how it looks. But actually, we're interested in this little square on which is the new data coming. And the only thing I need to do is to condition. And this is how the new covariance matrix would look like. This is not really interesting. What is interesting is this, that whenever I condition in the joint distribution with this kernel, I get a really, really nice fit. And there's no parameters here. This is a non-parametric way, but it is a nice, let's say, method to, let's say, do a density estimation or model these non-linear signals. and maybe I didn't comment on that but all of these kernels actually depend on parameters and these parameters actually are also important in your fit you cannot just take the one that you like the most because this will have an impact on the distribution on functions so if I take these parameters let's say in a different way I'm going to get different fits so L for this Gaussian for this square exponential somehow the scale. So if I take it like really big or really low, then this is going to, let's say, define how complex or how flexible my Gaussian process is, even if it's really, really smooth. So there are various ways of selecting these hyperparameters. I don't want to go into details. You can't do it by maximizing the marginal likelihood or via a cross-validation approach. But this is just point out that this is something that is going to be a could be learned and a the kernel space is a quite rich so a the basic kernels that i show which have some intuition can actually be combined to just get a huge zoo of kernels which actually this is where the thinking happened because no one is going to give you a kernel coming from heaven you need to think about the problem and to see how you model your signal is not just doing model.fit. This is, in my opinion, where in the Gaussian process methods, the thinking happens, selecting the kernel. For example, how complex do you want it to be, how smooth, and so on and so forth. And there's a rich theory about how these kernels can be constructed and, let's say, how different spectral properties would give different kernels for different types of applications. Anyway, I just want to say that it's a really rich space, and actually sums, products, convolutions, direct sums, direct products, you can get new kernels from old kernels. So just to finalize, I want to just show how you could do this in practice, and there are many ways of doing this. yesterday we have seen talks with Stan, also we can do it yourself, I'm going to show you how to do it with Scikit-Learn unfortunately the resolution is not the best, but I will share this anyway so what is important is that here I have a periodic signal with noise and the natural candidate for that is just an exponential sine squared function and you can do this with Scikit-Learn just by getting the kernels and specifying the parameters. And in scikit-learn, you can define the bounds so that they can optimize for the best parameters by maximizing the marginal likelihood. If you add a linear trend, and here the red line specifies where I had the training and the test set, one way of modeling that is with exponential kernel or a radial basis function. And here you take a really long scale so that this, let's say, captures the long-term movement. So just by adding these two, you can model the signal quite well. Yeah, so I can add yet another seasonal component so that it mixes with the other one in such a way that the periods are different. And how do I model this new signal? Well, I just add another periodic kernel. And if I want to, depending on how I see the signal and where, like, domain knowledge and knowing where your data is generated from, this is where you will get insight of how to select the kernel better. And I don't think this is an easy task, but in some of these applications, it has been proven to be quite successful. And it gives you, even though it's not a parametric model, you understand what's happening if you really think carefully about how to select the kernel. Yeah, so there are some computational challenges. So when people talk about Gaussian process, there's always this question that in order to do inference, you're going to use essentially all of the data because you need to have this huge covariance matrix and then you take some inverse of certain matrices. So yeah, that's an important topic in Gaussian processes, but nevertheless, there are ways of doing this in a clever way. So try not to invert the matrix, but use some matrix factorization algorithm or do this by approximation. So that, by itself, is a quite interesting topic. About some reference, I can just recommend the book of Gelman and others on Bayesian data analysis. They have a chapter devoted on Gaussian process, and they go really fast on the definition of Gaussian process, but they spend most of the chapter studying how kernels should be thought. I have written a little bit of experiments about how to show this so that you can play around with it and see how Gaussian process, let's say, works from scratch. So from the basic linear regression towards doing this time series forecasting. And most of the basis of my talk, actually, was trying to understand Chapter 2 of this book on Gaussian process for machine learning of Rasmussen and Williams. So this book, I mean, this chapter has, like, I don't know, 10 pages or even less, but it took me a decent amount of time to run these experiments and convince myself that what they claim makes sense, so to say. And this is why I decided to give this talk. Thank you very much. If you have any questions, just let me know. Thank you. Okay, any questions? Thank you. You only showed examples for one-dimensional time series. Can you use Gaussian processes to also fit two-dimensional or three-dimensional time series? Yeah, actually, when you see the definition on Wikipedia of Gaussian process, they don't have a time. They can define Gaussian process in any index set. So in principle, if you have an assignment of give me a two-dimensional or n-dimensional thing and you can give me a normal distribution, then you can use it in this context. So, yeah. So, Gaussian processes are defined either for, let's say, more dimensional indices space. Thank you for the talk. I have one question concerning the data that you were fitting. So, mostly what I have seen that you were fitting some data which can be fit also with standard Fourier transform or the discrete value of retransform, how good they are fitting the data compared to the Gaussian process? I haven't experimented with that. I did this for myself, so I always take the easiest example first. What I do know is that you can relate a Gaussian processor with splines, so you could also try to do, okay, because essentially I could also fit this with the splines. And for the spline, maybe you need to add the degree of smoothness or how many nodes. In this book, they claim that Gaussian processes could be a better option there because you don't need to select this. But in case it depends on the application, so I don't have a concrete answer of why using Gaussian processes against these ones. Okay, yeah. Hi. So you have showed that the kernels have a kind of algebra, so you can add them and stuff. And then in applications, you figure out the best combination. So is there also a kind of basis within these kernels? I mean, that you just can stick with the linear kernels or the polynomial kernels, and then you know that it will converge eventually? Or do you have to... Well, I guess, let's say, once you talk about function spaces, this becomes like an infinite dimensional space, right? So I guess you could try to use some Hilbert space decomposition, but I don't think this is going to help you in practice because, yeah, then you need to specify, okay, how do I project these spaces into my concrete application? What I've seen is that people will go into the data and try to see how it was generated, get domain knowledge, and from that, for example, if you see that the amplitude should be increased with the trend, you will multiply your square sine kernel with a rated basis kernel. And you believe that smoothness is going to be something important for you, then they can control it by adding a rational quadratic. So, yeah, I guess it's more on the data knowledge and domain expertise rather than a general framework. Okay, last question. Thank you. Just a short one here. So I just discussed with the other author, Vincent, on the applicability of GPs outside academia. What is your opinion on GPs and their use in industry? Yeah, so actually I've used them. I think they've come really handy for modeling, at least in the applications that I've done, seasonality. Like strange seasonality, again, you can use splines. And, yeah, I guess you, in some sense, also get, like, the credible interval for splines. What I like about, for example, using GARCH processes is that you can encode a lot of information in the kernel with a lot of intuition, and you get the credible intervals for it. So if I need to remove a strange seasonality, which doesn't look like a sine function, I would just throw a random process and model it like this. but I do know a lot of people using Gaussian processing in practice and I guess for them the challenge is about the computations but I do know they are still quite active in industry Okay, thank you very much Juan very interesting talk and another round of applause