I am familiar with basic regression theory and numerical optimization; I have a qualitative understanding of VC bound. Given this background:
What is the number of parameters in a GP model?
A parametric regression/classification model has some number of parameters, and having estimated them I only need to carry them around. If a model has more parameters, it is less parsimonious and less desirable than a model with fewer parameters that achieves the same error.
If I add regularization/bias/prior to a model, it reduces the dimension of parameter-space.
In a GP model I need to hold on to my entire training set. But clearly the number of parameters I'm fitting is much smaller. So how do I quantify that?
To be clear, I'm not talking about the hyperparameters of the covariance function.
I am not exactly sure what you are asking here, but I think that is more a terminology disconnect than anything. I haven't even heard someone say Vapnik-Chervonenkis dimension in a long time, so we may be coming at this from two different standpoints. But I'll try to answer as best I can.
When you say parametric regression model, do you mean something like s(x) = a + bx where the parameters here are a and b? If that's the case then I might say that a GP model has N parameters, where there are N pieces of data. You could argue (and I think Rasmussen does) that GPs are nonparametric, since there are N parameters and N pieces of data, but I don't usually use that term about GPs.
As you point out, you do not need to carry around the data after you have fit something like a s(x) = a + bx model ... you only need a and b to make predictions. This is not true for GPs, where you will always need the {x_1, ..., x_N} points in order to make predictions. The benefit to using these kernels as your basis s(x) = c_1K(x, x_1) + ... + c_NK(x, x_N) instead of s(x) = a + bx + cy + dx^2 + ... is that there is a guaranteed unique interpolant to your data. It would depend on how you define "error" (and this is a constant issue in any sort of real-world setting) but if you are defining error by how well your model fits the data then a GP model would perfectly fit any data. I'm not saying that is really what you want (I talk in http://blog.sigopt.com/post/132959177928/sigopt-fundamentals... about why that is a bad idea at times), but I am suggesting that there are certain properties that GPs satisfy very nicely.
I would need to understand why you are saying that the number of parameters you are fitting is much smaller. I don't immediately know how you can guarantee that just by looking at a bunch of x and y values. If there is a specific reason you think your data came from a very simple function, then that would make sense, but if you can't make that assumption I don't know why it would immediately be preferable to use one model over another.
Given that, if you do think that your data came from a relatively low parameter source (to put it in tangible terms, suppose you have 20000 pieces of data and you think you should only need 50 components in your s function) then you can create a 50 term model using reproducing kernels instead of polynomials (or whatever else). This is often called an RBF network (or kernel network). You can find a great reference online (http://www.eng.buffalo.edu/~psingla/Teaching/MAE502/ReadingP...) or a decent reference in chapter 18 of my book (kernel-based approximation methods in Matlab). If you choose to use this strategy you would be creating a model where you would only have to carry around 50 parameters, instead of the 20000 that you would have to use in the GP setting. The benefit over polynomials would be flexibility and locality, the benefit over GPs would be cost (50<<20000 so predictions are cheaper). There are also Support Vector Machines/Regression which do this stuff in a sparse framework (less that 20000 basis functions) using reproducing kernels.
Basically, I would totally agree with you that simpler is better. A the end of the day, I think of GPs as a useful tool for identifying local behavior, but I try not to use them globally (although in this paper I looked at the global limit, http://epubs.siam.org/doi/abs/10.1137/110824784). Probably the best strategy, if you believe that your data comes from a low complexity function and should only ne...
With N parameters for N data points I have far, far too many degrees of freedom and have grossly overfit my training set. The only way to guarantee good generalization outside one's training set is to have a lot of data compared to the number of degrees of freedom in the model.
To avoid this, we use bias/prior/regularization which reduces the effective number of degrees of freedom. For example, fitting a quadratic polynomial to some data, we have 3 degrees of freedom including the constant term. If I add a regularization term \lambda*||a|| where a is the vector of coefficients, I can alleviate overfitting. We can determine an optimal value for \lambda using validation or cross-validation. For \lambda=0, we have 3 parameters. As \lambda->infinity, we force the answer closer and closer to the constant 0 polynomial ie. no degrees of freedom.
In other words, regularization avoids overfitting by reducing the dimension of the model.
I am not worried about the amount of data I have to carry around; rather, it was a way of thinking about the dimension of parameter-space.
Don't big alarm bells go off in your head when you hear N parameters for N data points?
> Don't big alarm bells go off in your head when you hear N parameters for N data points?
Gaussian process models typically also have hyperparameters that control things like the smoothness of the learned surface. My hand-wavey sense of how GPs work is that you have a prior distribution over functions (think curves) and a hyperparameter controlling the smoothness/length-scale of those functions. You observe data, then find the posterior distribution over functions that make your data most likely - i.e. of all these possible functions, which ones have the best chance of actually generating the data you're observing? The shape of that posterior distribution depends on the set of N points you observe, and it will be different if you have different data. In that sense, then, the data are parameters, since the posterior depends exactly on the set of data.
The analogy to traditional regression comes by having the hyperparameters (controlling the covariance function IIRC, i.e. the smoothness of the posterior curves), which allows the resulting posterior distribution to explain the data perfectly (overfitting, sharp curves) or not (generalizability, smooth curves). See for example slides 23-24 of this slide deck:
Well, this is the difference in phrasing between the traditional statistical modeling setting and the terminology used in numerical analysis; perhaps this is why GPs are sometimes referred to as a nonparametric model. Also, I think you're oversimplifying the role of regularization, but that's more of a semantic than salient point.
You are correct that you should not use "N parameters for N data points", but the role that covariance kernels play is not the same as the role played by (as you suggest) polynomials. Covariance kernels have some sort of localized behavior which limits the degree to which they interact with the data. Not to limit the discussion to things I'm comfortable with, but I wrote some of this up at a very high level here http://blog.sigopt.com/post/131498069303/sigopt-fundamentals.... As a result, I think one could argue that the compactness of their domain (practical compactness, regardless of their potentially infinite support) produces the actual dimension of the parameter space. As such, different locations in the domain may have different dimensionality, though most of the time the swift decay of kernels away from their center will lead to a much lower "effective dimensionality" than N.
Maybe a good way to think of this is the following: Why are you using a quadratic polynomial to fit your (presumably, if you have only 3 parameters) 1D data? Probably just because you feel like it - which is as good a reason as any, short of some rigorous statistical testing. I use kernels because I like their theoretical properties, not because "the data demands it" or some more rigorous reason. But suppose you actually thought of the biggest polynomial you could use: an N-1 degree polynomial:
That's the biggest (unique) polynomial you could use to fit your data. You feel, probably correctly, that N terms is too many, so you zero out some of these terms by manually setting most of the coefficients equal to 0 and solving (somehow) for the rest:
This is one, very common and very useful way, to fit data, and it is the first thing I teach when I introduce people to approximation theory. But think about the main drawback behind this method: this is the same function everywhere that you want to make approximations. You have no ability to do anything different in different regions of your domain. You could have a problem for which this is true, but I think a lot of problems have some degree of spatial variability/locality, and this is where kernel methods shine. Kernel methods propose that you fit a model that looks like
Now, does this require you to invert an NxN matrix - absolutely. But because there is some locality in these K functions, different kernels will end up being 0 in different locations. Suppose that x_4 through x_N are very far from x_1; your interpolation in that setting looks like:
s(x) = a_1K(x,x_1) + a_2K(x,x_2) + a_3K(x,x_3) + a_4 0 + ... + a_N 0 , for x close to x_1
Notice I use the term "close" here; the decision of what close means is the hyperparameter solution that this blog post was trying to allude to, and that this other post http://blog.sigopt.com/post/134931028143/sigopt-in-depth-pro... addresses more thoroughly. These summation values are probably not identically 0 (although it could be, if you like the Wendland or Wu kernels) but it may be close enough to think of it as 0. Suppose, instead, that you are over near x_N and that only x_4 is close to x_N:
> Notice that both the polynomial and kernel methods have less than N terms which are nonzero, but [...]
starts to answer my worry a little bit. It really is not the case that there are N true degrees of freedom, we just nominally started out with that many and then gave up certain ones of them.
What I would love to have is a quantitative analysis of the "number of free parameters" in my model after I have fit it, so that I can compare it to (say) a 50-parameter polynomial model that is equally good at reproducing the training set.
I am also very interested in knowing what to do when the dimension of the input space can vary.
For example, suppose that I'm interested in learning 4-body gravitational motion using 3-body training data ie. predict total energy of a system given mutual distances.
Notwithstanding the fact that this is trivial to compute directly, how do I set this up as a GP? Are there general strategies for this? What do I look up on google scholar, or what application field of ML most deals with this?
> I am also very interested in knowing what to do when the dimension of the input space can vary.
It is frequently the case that one's data lies on a (curvy) low dimensional manifold embedded within a higher dimensional space. Because of this, performing an orthogonal projection to a lower dimension is often a bad idea. Consider the famous "Swiss roll" dataset: http://www.convexoptimization.com/dattorro/manifold_learning... Although the data lies on a 2D manifold, there exists no affine subspace to project the data onto; one must instead perform a non-linear mapping.
Problems like this are tackled using a variety of methods that attempt to "unroll" the manifold into its own vector space. The rest of the problem can then be solved using the tools of linear algebra. See http://repository.upenn.edu/cgi/viewcontent.cgi?article=1001... for example.
I have a question that is only tangentially related to this so feel free to ignore.
Say I flip 50 fair coins, P(heads)=P(tails)=0.5, 10 times each a day for 100 days. For each coin we have 100 records like HTHTHHTTHT. Now I count the number of tails (5 in the example) for each day for each coin. So each coin has a record with length 100 like [5, 8, 5, ...,6]. Now I can calculate the variance of these values for each coin and end up with 50 variances. What do we expect the distribution of these variances to look like?
I need a bit of clarification here. In the list of length 100, each entry in that list is a binomial random variable X~Binom(10, .5). The variance of each of these binomials is Var(X) = 10(.5)(.5) = 2.5.
But I don't think that's the variance you're talking about. Do you mean the sample variance of the list of length 100? I'm sure that we can figure that out too (or at least what it should look like) but I want to confirm that's the variance you are interested in. Off the top of my head, I think this sample variance should have a roughly Chi-Square distribution, since Binom(10, .5) is close to normal (sort of). In the Wikipedia article, it explains how Chi-Square can be expressed a sum of normals squared (https://en.wikipedia.org/wiki/Chi-squared_distribution). I could maybe figure out the actual distribution if you really need it.
Then you ask what should the distribution look like? Is this answered by the answer above, or did you mean something related to the 50 samples in particular?
Thanks, I'm wondering how to calculate the distribution of sample variances. Yes, each variance is calculated from the list of 100 counts of #tails. It's just something that has bothered me for awhile.
I am a frequentist by training and got a little confused by Bayesian and engineering terminology while flipping through the posts. Bear with me.
Did I understand well that kernel interpolation is what we'd call a kind of non-noisy kernel regression? If it's the case and dimension d of the regressors is large, multivariate non-parametric estimation will have very slow convergence, won't it?
I also get confused by the notation, so you are not alone.
I think it is fair to say that kernel interpolation is a non-noisy kernel regression. Of course, this would also depend on your choice of terminology ...
I use the word regression to mean any sort of fitting of specific basis functions to data - the fit doesn't have to be perfect. Interpolation, at least in my mind, means that the fit does have to be perfect.
So when I say kernel interpolation I think it is fine to think of it as regression that perfectly passes through any data points.
Technically, I like to think of the fact that the data has noise as independent of the strategy that you use to fit it. So you could perfectly interpolate noisy data (a bad idea) and you could do a regression on noise-less data (also bad, but not as bad). But yeah, I think it's probably safe to say "non-noisy kernel regression".
As far as your real question, chapter 9 of my book talks about that, and the reason I mention that here is that it is a complicated question. Just to confirm, when you say dimension you mean physical dimension, not the number of pieces of data, right?
The classical convergence bounds all deal with the smoothness of the data (or rather the function that generated the data) and the smoothness of the reproducing kernel. The relevance of the smoothness is augmented by the dimension, so the dimension can be relevant. Arguably, this is why statisticians like Michael Stein would say that Gaussians are not appropriate for low dimensions (too much smooothness) but are necessary for high dimensions (because without the smoothness the convergence is too slow).
Some recent results (http://epubs.siam.org/doi/abs/10.1137/10080138X) talk about dimension independent bounds, but with some caveats, and almost entirely from the numerical analysis standpoint (which can be tough to read).
Sorry to sort of fork the answer into pieces there. Anything else I can help with?
Thanks for the answer and the SIAM reference, I did manage to pick up something interesting from there, I think. Smoothness, dimensionality and also adaptivity are vast topics indeed. Good luck with the future work!
I tried to do some digging to find easy to access (both through the web, but also not terribly complicated) references on error bounds for kernels interpolation. I don't think there is one ... although the internet will surely correct me if I'm wrong.
The original source that I often return to (even before I look in my book) is Wendland's 2005 book "Scattered Data Approximation". But that book is really tough to read, even for me. Lots of tough math. Fasshauer's 2007 book "Meshfree Approximation Methods in Matlab" is much easier to read, but references Wendland's book for most of the heavy lifting.
There is a newer branch of study on convergence dealing with what are called "sampling inequalities". In many ways, the math behind these is even worse as it usually requires a bunch of polynomial theory, but the results are more easily accessible. In particular, Christian Rieger and Barbara Wolmuth have outstanding content on this which is readily available on the web. In particular, I was able to find Christian's PhD thesis (https://www.deutsche-digitale-bibliothek.de/binary/JOWXLCGSV...) which provides a good journey from the start to actual results.
Sorry I can't provide a cleaner statement about this. Maybe the most basic point I can make about convergence here would be to reference the early part of Christian's thesis where he alludes to the theorem that says the quality of an interpolant sort of looks like:
error of interpolant = O(h^{k-d})
That is a gross simplification, but gives the gist of the result. h is the "fill distance" (or grid width for structured data) k is the smoothness of the kernel and d is the dimension of the data. This indicates why smooth kernels (large k) are needed to get the convergence for high d.
Excellent post. Would you have a recommendation on a resource for learning about reproducing kernel Hilbert spaces? In particular, I'm interested in learning more about how they are applied to data analysis. I've been reading about them in a very theoretical sense, but very few math texts provide "real-world" examples of how they are used.
You are correct that there are few texts that discuss real-world examples of reproducing kernel Hilbert spaces without a significant amount of overhead to get there. In reality, I'm not sure there is an easy way to talk about reproducing kernel Hilbert spaces without all that functional analysis, which is a tough topic. I would argue that's probably the greatest thing holding back kernels - their high barrier to entry.
I'll spend one line here to mention my book (kernel-based approximation methods in Matlab) and also say that it is pretty heavy on the theory; the second half of the book deals with reproducing kernels as used in a variety of fields, but getting there requires a lot of very tough math.
Lemme take a moment to list places (some more real-world than others) where I know that reproducing kernels pop up - and anyone else reading this please feel free to add to this list:
1) Approximation theory - My bread and butter and the most important role for kernels as far as I'm concerned. When you have scattered data (especially in higher dimensions) and you want to interpolate and predict unobserved values, reproducing kernels give you several very nice guarantees. Through this, they also pop up in computer vision or surface reconstruction from point cloud data.
2) Spatial statistics - They actually play an identical role in kriging as they do in approximation theory, albeit for totally different reasons. People who are doing remote sensing or topological mapping may use reproducing kernels (or software that uses these kernels.
3) Machine learning - Both RBF networks and support vector machines use reproducing kernels, although in two different ways. The mechanism by which they are used in support vector machines (feature maps) is difficult for me to understand which is why I don't talk about them much. They are still a very important field that uses kernels.
4) Differential equations - The same theory that makes reproducing kernels work for approximation theory also works for numerically solving differential equations. You could lump stuff regarding Green's functions in here too, although those may not be reproducing kernels.
5) Numerical integration - Especially for integration in higher dimensions, basically all the theory of convergence for quasi Monte Carlo methods relies on reproducing kernel Hilbert spaces.
Maybe the best thing I can point you to is the same place where I first learned about this stuff: math.iit.edu/~fass/590. That's an old (recently updated) class website from a class I took at the Illinois Institute of Technology when I was an undergrad there. Dr. Fasshauer does an outstanding job putting notes and slides for his class together.
If you're feeling adventurous, I would recommend our book, but a good starting point would be his class notes, which I still turn back to often and are written to introduce graduate students (who know nothing about kernels) to the topic. I'll work on putting together a list of good application-oriented sources as well.
Thanks for your interest! Did you have a specific example/topic of kernels from data analysis that piqued your interest? I could look into more targeted content.
That may be. And I don't think kernels will ever be a hot topic. But part of that is that they are a very old topic (you can see Gauss referring to them in slide 22 of http://math.iit.edu/~fass/590/notes/Notes590_Ch1.pdf). I think that kinda prevents funding agencies from pushing too hard for research into them, since there is a more likely quick return on investment into a new topic. And I don't begrudge them that - funding should be given to people with new ideas which have the potential to revolutionize things. I'd also mention that there are people trying to bridge the gap between how GPs work and the benefits of neural nets (http://arxiv.org/abs/1502.05700).
Probably another thing holding them back is their high barrier to entry - you can't do much with the theory of kernels unless you know statistics and functional analysis (and of course numerical analysis and linear algebra,) which is one of the reasons I had difficulty finding research students. To draw a parallel to numerical analysis (where I am much more comfortable than machine learning,) this is why finite element methods are much more popular than, say, boundary element methods for PDEs, despite the fact that boundary element methods are better in many circumstances. But when implementing things, especially in industry, simplicity and robustness count for something and those are not strong points of kernels.
But kernels haven't gone anywhere yet, and they've been fundamental to analysis since at least David Hilbert, so they are not going anywhere any time soon. Deep learning ... we'll see. Maybe it has the legs to stick around or maybe it'll be swept aside by the next hot thing. 8 years ago all the research funding was going into uncertainty quantification, 5 years ago compressed sensing was the hot idea (thanks Terry Tao,) and now deep learning is going to light the world on fire. All of those are, and should be, hotter topics than kernels, but we'll see what we're talking about in 5 years. It won't be kernels (unless something crazy happens and I become the leader of the US, in which case I have a great idea for a new reality show) but I don't know that it will be deep learning either.
26 comments
[ 2.8 ms ] story [ 68.0 ms ] threadWhat is the number of parameters in a GP model?
A parametric regression/classification model has some number of parameters, and having estimated them I only need to carry them around. If a model has more parameters, it is less parsimonious and less desirable than a model with fewer parameters that achieves the same error.
If I add regularization/bias/prior to a model, it reduces the dimension of parameter-space.
In a GP model I need to hold on to my entire training set. But clearly the number of parameters I'm fitting is much smaller. So how do I quantify that?
To be clear, I'm not talking about the hyperparameters of the covariance function.
When you say parametric regression model, do you mean something like s(x) = a + bx where the parameters here are a and b? If that's the case then I might say that a GP model has N parameters, where there are N pieces of data. You could argue (and I think Rasmussen does) that GPs are nonparametric, since there are N parameters and N pieces of data, but I don't usually use that term about GPs.
As you point out, you do not need to carry around the data after you have fit something like a s(x) = a + bx model ... you only need a and b to make predictions. This is not true for GPs, where you will always need the {x_1, ..., x_N} points in order to make predictions. The benefit to using these kernels as your basis s(x) = c_1K(x, x_1) + ... + c_NK(x, x_N) instead of s(x) = a + bx + cy + dx^2 + ... is that there is a guaranteed unique interpolant to your data. It would depend on how you define "error" (and this is a constant issue in any sort of real-world setting) but if you are defining error by how well your model fits the data then a GP model would perfectly fit any data. I'm not saying that is really what you want (I talk in http://blog.sigopt.com/post/132959177928/sigopt-fundamentals... about why that is a bad idea at times), but I am suggesting that there are certain properties that GPs satisfy very nicely.
I would need to understand why you are saying that the number of parameters you are fitting is much smaller. I don't immediately know how you can guarantee that just by looking at a bunch of x and y values. If there is a specific reason you think your data came from a very simple function, then that would make sense, but if you can't make that assumption I don't know why it would immediately be preferable to use one model over another.
Given that, if you do think that your data came from a relatively low parameter source (to put it in tangible terms, suppose you have 20000 pieces of data and you think you should only need 50 components in your s function) then you can create a 50 term model using reproducing kernels instead of polynomials (or whatever else). This is often called an RBF network (or kernel network). You can find a great reference online (http://www.eng.buffalo.edu/~psingla/Teaching/MAE502/ReadingP...) or a decent reference in chapter 18 of my book (kernel-based approximation methods in Matlab). If you choose to use this strategy you would be creating a model where you would only have to carry around 50 parameters, instead of the 20000 that you would have to use in the GP setting. The benefit over polynomials would be flexibility and locality, the benefit over GPs would be cost (50<<20000 so predictions are cheaper). There are also Support Vector Machines/Regression which do this stuff in a sparse framework (less that 20000 basis functions) using reproducing kernels.
Basically, I would totally agree with you that simpler is better. A the end of the day, I think of GPs as a useful tool for identifying local behavior, but I try not to use them globally (although in this paper I looked at the global limit, http://epubs.siam.org/doi/abs/10.1137/110824784). Probably the best strategy, if you believe that your data comes from a low complexity function and should only ne...
To avoid this, we use bias/prior/regularization which reduces the effective number of degrees of freedom. For example, fitting a quadratic polynomial to some data, we have 3 degrees of freedom including the constant term. If I add a regularization term \lambda*||a|| where a is the vector of coefficients, I can alleviate overfitting. We can determine an optimal value for \lambda using validation or cross-validation. For \lambda=0, we have 3 parameters. As \lambda->infinity, we force the answer closer and closer to the constant 0 polynomial ie. no degrees of freedom.
In other words, regularization avoids overfitting by reducing the dimension of the model.
I am not worried about the amount of data I have to carry around; rather, it was a way of thinking about the dimension of parameter-space.
Don't big alarm bells go off in your head when you hear N parameters for N data points?
Gaussian process models typically also have hyperparameters that control things like the smoothness of the learned surface. My hand-wavey sense of how GPs work is that you have a prior distribution over functions (think curves) and a hyperparameter controlling the smoothness/length-scale of those functions. You observe data, then find the posterior distribution over functions that make your data most likely - i.e. of all these possible functions, which ones have the best chance of actually generating the data you're observing? The shape of that posterior distribution depends on the set of N points you observe, and it will be different if you have different data. In that sense, then, the data are parameters, since the posterior depends exactly on the set of data.
The analogy to traditional regression comes by having the hyperparameters (controlling the covariance function IIRC, i.e. the smoothness of the posterior curves), which allows the resulting posterior distribution to explain the data perfectly (overfitting, sharp curves) or not (generalizability, smooth curves). See for example slides 23-24 of this slide deck:
https://www.cs.toronto.edu/~hinton/csc2515/notes/gp_slides_f...
Rasmussen has some useful video tutorials of all this. It's been a few years since I've watched them but if you're interested I'd recommend them.
You are correct that you should not use "N parameters for N data points", but the role that covariance kernels play is not the same as the role played by (as you suggest) polynomials. Covariance kernels have some sort of localized behavior which limits the degree to which they interact with the data. Not to limit the discussion to things I'm comfortable with, but I wrote some of this up at a very high level here http://blog.sigopt.com/post/131498069303/sigopt-fundamentals.... As a result, I think one could argue that the compactness of their domain (practical compactness, regardless of their potentially infinite support) produces the actual dimension of the parameter space. As such, different locations in the domain may have different dimensionality, though most of the time the swift decay of kernels away from their center will lead to a much lower "effective dimensionality" than N.
Maybe a good way to think of this is the following: Why are you using a quadratic polynomial to fit your (presumably, if you have only 3 parameters) 1D data? Probably just because you feel like it - which is as good a reason as any, short of some rigorous statistical testing. I use kernels because I like their theoretical properties, not because "the data demands it" or some more rigorous reason. But suppose you actually thought of the biggest polynomial you could use: an N-1 degree polynomial:
s(x) = a_1 + a_2x + a_3x^2 + a_4x^3 + ... + a_Nx^{N-1}
That's the biggest (unique) polynomial you could use to fit your data. You feel, probably correctly, that N terms is too many, so you zero out some of these terms by manually setting most of the coefficients equal to 0 and solving (somehow) for the rest:
s(x) = a_1 + a_2x + a_3x^2 + 0x^3 + ... + 0x^{N-1}
This is one, very common and very useful way, to fit data, and it is the first thing I teach when I introduce people to approximation theory. But think about the main drawback behind this method: this is the same function everywhere that you want to make approximations. You have no ability to do anything different in different regions of your domain. You could have a problem for which this is true, but I think a lot of problems have some degree of spatial variability/locality, and this is where kernel methods shine. Kernel methods propose that you fit a model that looks like
s(x) = a_1K(x,x_1) + a_2K(x,x_2) + a_3K(x,x_3) + a_4K(x,x_4) + ... + a_NK(x,x_N)
Now, does this require you to invert an NxN matrix - absolutely. But because there is some locality in these K functions, different kernels will end up being 0 in different locations. Suppose that x_4 through x_N are very far from x_1; your interpolation in that setting looks like:
s(x) = a_1K(x,x_1) + a_2K(x,x_2) + a_3K(x,x_3) + a_4 0 + ... + a_N 0 , for x close to x_1
Notice I use the term "close" here; the decision of what close means is the hyperparameter solution that this blog post was trying to allude to, and that this other post http://blog.sigopt.com/post/134931028143/sigopt-in-depth-pro... addresses more thoroughly. These summation values are probably not identically 0 (although it could be, if you like the Wendland or Wu kernels) but it may be close enough to think of it as 0. Suppose, instead, that you are over near x_N and that only x_4 is close to x_N:
s(x) = a_1 0 + a_2 0...
> Notice that both the polynomial and kernel methods have less than N terms which are nonzero, but [...]
starts to answer my worry a little bit. It really is not the case that there are N true degrees of freedom, we just nominally started out with that many and then gave up certain ones of them.
What I would love to have is a quantitative analysis of the "number of free parameters" in my model after I have fit it, so that I can compare it to (say) a 50-parameter polynomial model that is equally good at reproducing the training set.
One short answer I could point you to regarding this question is a topic called "local interpolation" using a compactly defined Lagrange basis with a prescribed size (perhaps size 50) of the parameter space. I talk about this in the context of kernel-based finite difference methods in remark 19.8 of my book: some of the references I point to there are (http://epubs.siam.org/doi/abs/10.1137/090769570, Hangelbroek is an outstanding author) and (http://www.sciencedirect.com/science/article/pii/S0377042711...). I think it's also possible to look at this question from the framework of moving least squares/approximate approximation (http://link.springer.com/chapter/10.1007/978-3-642-56103-0_8) or maybe even quasi-interpolation, depending on how you think about it (http://link.springer.com/article/10.1007/BF01279020).
If you'd like to talk about this offline, email me (see the blog post). I'll still post any thoughts I have here.
For example, suppose that I'm interested in learning 4-body gravitational motion using 3-body training data ie. predict total energy of a system given mutual distances.
Notwithstanding the fact that this is trivial to compute directly, how do I set this up as a GP? Are there general strategies for this? What do I look up on google scholar, or what application field of ML most deals with this?
Thanks a lot for your help.
It is frequently the case that one's data lies on a (curvy) low dimensional manifold embedded within a higher dimensional space. Because of this, performing an orthogonal projection to a lower dimension is often a bad idea. Consider the famous "Swiss roll" dataset: http://www.convexoptimization.com/dattorro/manifold_learning... Although the data lies on a 2D manifold, there exists no affine subspace to project the data onto; one must instead perform a non-linear mapping.
Problems like this are tackled using a variety of methods that attempt to "unroll" the manifold into its own vector space. The rest of the problem can then be solved using the tools of linear algebra. See http://repository.upenn.edu/cgi/viewcontent.cgi?article=1001... for example.
Say I flip 50 fair coins, P(heads)=P(tails)=0.5, 10 times each a day for 100 days. For each coin we have 100 records like HTHTHHTTHT. Now I count the number of tails (5 in the example) for each day for each coin. So each coin has a record with length 100 like [5, 8, 5, ...,6]. Now I can calculate the variance of these values for each coin and end up with 50 variances. What do we expect the distribution of these variances to look like?
But I don't think that's the variance you're talking about. Do you mean the sample variance of the list of length 100? I'm sure that we can figure that out too (or at least what it should look like) but I want to confirm that's the variance you are interested in. Off the top of my head, I think this sample variance should have a roughly Chi-Square distribution, since Binom(10, .5) is close to normal (sort of). In the Wikipedia article, it explains how Chi-Square can be expressed a sum of normals squared (https://en.wikipedia.org/wiki/Chi-squared_distribution). I could maybe figure out the actual distribution if you really need it.
Then you ask what should the distribution look like? Is this answered by the answer above, or did you mean something related to the 50 samples in particular?
Did I understand well that kernel interpolation is what we'd call a kind of non-noisy kernel regression? If it's the case and dimension d of the regressors is large, multivariate non-parametric estimation will have very slow convergence, won't it?
I think it is fair to say that kernel interpolation is a non-noisy kernel regression. Of course, this would also depend on your choice of terminology ... I use the word regression to mean any sort of fitting of specific basis functions to data - the fit doesn't have to be perfect. Interpolation, at least in my mind, means that the fit does have to be perfect.
So when I say kernel interpolation I think it is fine to think of it as regression that perfectly passes through any data points.
Technically, I like to think of the fact that the data has noise as independent of the strategy that you use to fit it. So you could perfectly interpolate noisy data (a bad idea) and you could do a regression on noise-less data (also bad, but not as bad). But yeah, I think it's probably safe to say "non-noisy kernel regression".
As far as your real question, chapter 9 of my book talks about that, and the reason I mention that here is that it is a complicated question. Just to confirm, when you say dimension you mean physical dimension, not the number of pieces of data, right?
The classical convergence bounds all deal with the smoothness of the data (or rather the function that generated the data) and the smoothness of the reproducing kernel. The relevance of the smoothness is augmented by the dimension, so the dimension can be relevant. Arguably, this is why statisticians like Michael Stein would say that Gaussians are not appropriate for low dimensions (too much smooothness) but are necessary for high dimensions (because without the smoothness the convergence is too slow).
Some recent results (http://epubs.siam.org/doi/abs/10.1137/10080138X) talk about dimension independent bounds, but with some caveats, and almost entirely from the numerical analysis standpoint (which can be tough to read).
Sorry to sort of fork the answer into pieces there. Anything else I can help with?
The original source that I often return to (even before I look in my book) is Wendland's 2005 book "Scattered Data Approximation". But that book is really tough to read, even for me. Lots of tough math. Fasshauer's 2007 book "Meshfree Approximation Methods in Matlab" is much easier to read, but references Wendland's book for most of the heavy lifting.
There is a newer branch of study on convergence dealing with what are called "sampling inequalities". In many ways, the math behind these is even worse as it usually requires a bunch of polynomial theory, but the results are more easily accessible. In particular, Christian Rieger and Barbara Wolmuth have outstanding content on this which is readily available on the web. In particular, I was able to find Christian's PhD thesis (https://www.deutsche-digitale-bibliothek.de/binary/JOWXLCGSV...) which provides a good journey from the start to actual results.
Sorry I can't provide a cleaner statement about this. Maybe the most basic point I can make about convergence here would be to reference the early part of Christian's thesis where he alludes to the theorem that says the quality of an interpolant sort of looks like:
error of interpolant = O(h^{k-d})
That is a gross simplification, but gives the gist of the result. h is the "fill distance" (or grid width for structured data) k is the smoothness of the kernel and d is the dimension of the data. This indicates why smooth kernels (large k) are needed to get the convergence for high d.
You are correct that there are few texts that discuss real-world examples of reproducing kernel Hilbert spaces without a significant amount of overhead to get there. In reality, I'm not sure there is an easy way to talk about reproducing kernel Hilbert spaces without all that functional analysis, which is a tough topic. I would argue that's probably the greatest thing holding back kernels - their high barrier to entry.
I'll spend one line here to mention my book (kernel-based approximation methods in Matlab) and also say that it is pretty heavy on the theory; the second half of the book deals with reproducing kernels as used in a variety of fields, but getting there requires a lot of very tough math.
Lemme take a moment to list places (some more real-world than others) where I know that reproducing kernels pop up - and anyone else reading this please feel free to add to this list:
1) Approximation theory - My bread and butter and the most important role for kernels as far as I'm concerned. When you have scattered data (especially in higher dimensions) and you want to interpolate and predict unobserved values, reproducing kernels give you several very nice guarantees. Through this, they also pop up in computer vision or surface reconstruction from point cloud data.
2) Spatial statistics - They actually play an identical role in kriging as they do in approximation theory, albeit for totally different reasons. People who are doing remote sensing or topological mapping may use reproducing kernels (or software that uses these kernels.
3) Machine learning - Both RBF networks and support vector machines use reproducing kernels, although in two different ways. The mechanism by which they are used in support vector machines (feature maps) is difficult for me to understand which is why I don't talk about them much. They are still a very important field that uses kernels.
4) Differential equations - The same theory that makes reproducing kernels work for approximation theory also works for numerically solving differential equations. You could lump stuff regarding Green's functions in here too, although those may not be reproducing kernels.
5) Numerical integration - Especially for integration in higher dimensions, basically all the theory of convergence for quasi Monte Carlo methods relies on reproducing kernel Hilbert spaces.
Maybe the best thing I can point you to is the same place where I first learned about this stuff: math.iit.edu/~fass/590. That's an old (recently updated) class website from a class I took at the Illinois Institute of Technology when I was an undergrad there. Dr. Fasshauer does an outstanding job putting notes and slides for his class together.
If you're feeling adventurous, I would recommend our book, but a good starting point would be his class notes, which I still turn back to often and are written to introduce graduate students (who know nothing about kernels) to the topic. I'll work on putting together a list of good application-oriented sources as well.
Thanks for your interest! Did you have a specific example/topic of kernels from data analysis that piqued your interest? I could look into more targeted content.
Probably another thing holding them back is their high barrier to entry - you can't do much with the theory of kernels unless you know statistics and functional analysis (and of course numerical analysis and linear algebra,) which is one of the reasons I had difficulty finding research students. To draw a parallel to numerical analysis (where I am much more comfortable than machine learning,) this is why finite element methods are much more popular than, say, boundary element methods for PDEs, despite the fact that boundary element methods are better in many circumstances. But when implementing things, especially in industry, simplicity and robustness count for something and those are not strong points of kernels.
But kernels haven't gone anywhere yet, and they've been fundamental to analysis since at least David Hilbert, so they are not going anywhere any time soon. Deep learning ... we'll see. Maybe it has the legs to stick around or maybe it'll be swept aside by the next hot thing. 8 years ago all the research funding was going into uncertainty quantification, 5 years ago compressed sensing was the hot idea (thanks Terry Tao,) and now deep learning is going to light the world on fire. All of those are, and should be, hotter topics than kernels, but we'll see what we're talking about in 5 years. It won't be kernels (unless something crazy happens and I become the leader of the US, in which case I have a great idea for a new reality show) but I don't know that it will be deep learning either.