{"id":1700,"date":"2023-12-01T17:25:39","date_gmt":"2023-12-01T17:25:39","guid":{"rendered":"https:\/\/www.ethanepperly.com\/?p=1700"},"modified":"2023-12-01T17:25:41","modified_gmt":"2023-12-01T17:25:41","slug":"five-interpretations-of-kernel-quadrature","status":"publish","type":"post","link":"https:\/\/www.ethanepperly.com\/index.php\/2023\/12\/01\/five-interpretations-of-kernel-quadrature\/","title":{"rendered":"Five Interpretations of Kernel Quadrature"},"content":{"rendered":"\n

I’m excited to share that my paper Kernel quadrature with randomly pivoted Cholesky<\/a><\/em>, joint with Elvira Moreno<\/a>, has been accepted to NeurIPS 2023 as a spotlight.<\/p>\n\n\n\n

Today, I want to share with you a little about the kernel quadrature problem. To avoid this post getting too long, I’m going to write this post assuming familiarity with the concepts of reproducing kernel Hilbert spaces<\/a> and Gaussian processes<\/a>.<\/p>\n\n\n\n

Integration and Quadrature<\/h2>\n\n\n\n

Integration is one of the most widely used operations in mathematics and its applications. As such, it is a basic problem of wide interest to develop numerical methods for evaluating integrals<\/a>.<\/p>\n\n\n\n

In this post, we will consider a quite general integration problem. Let \"\Omega\subseteq \real^d\" be a domain and let \"\mu\" be a (finite<\/a> Borel<\/a>) measure<\/a> on \"\Omega\". We consider the task of evaluating

  <\/span>   <\/span>\"\[I[f] = \int_\Omega f(x) g(x) \, \mathrm{d}\mu(x).\]\"<\/p>One can imagine that \"g\", \"\mu\", and \"\Omega\" are fixed, but we may want to evaluate this same integral \"I[f]\" for multiple different functions \"f\".<\/p>\n\n\n\n

To evaluate, we will design a quadrature approximation to the integral \"I[f]\":

  <\/span>   <\/span>\"\[\hat{I}_{w,s}[f] = \sum_{i=1}^n w_i f(s_i) \approx I[f].\]\"<\/p>Concretely, we wish to find real numbers \"w = (w_1,\ldots,w_n) \in \real^n\" and points \"s = (s_1,\ldots,s_n) \in \Omega^n\" such that the approximation \"\hat{I}_{w,s}[f] \approx I[f]\" is accurate.<\/p>\n\n\n\n

Smoothness and Reproducing Kernel Hilbert Spaces<\/h2>\n\n\n\n

As is frequently the case in computational mathematics<\/a>, the accuracy we can expect for this integration problem depends on the smoothness of the integrand \"f\". The more smooth \"f\" is, the more accurately we can expect to compute \"I[f]\" for a given budget of computational effort.<\/p>\n\n\n\n

In this post, will measure smoothness using the reproducing kernel Hilbert space<\/a> (RKHS) formalism. Let \"\mathcal{H}\" be an RKHS with norm<\/a> \"\norm{\cdot}\". We can interpret the norm as assigning a roughness \"\norm{f}\" to each function \"f\". If \"\norm{f}\" is large, then \"f\" is rough; if \"\norm{f}\" is small, then \"f\" is smooth.<\/p>\n\n\n\n

Associated to the RKHS \"\mathcal{H}\" is the titular reproducing kernel<\/a> \"k\". The kernel is a bivariate function \"k:\Omega\times\Omega\to\real\". It is related to the RKHS inner product<\/a> \"\langle\cdot,\cdot\rangle\" by the reproducing property

  <\/span>   <\/span>\"\[f(x)=\langle f, k(x,\cdot)\rangle \quad \text{for every }f\in\mathcal{H},x\in\Omega.\]\"<\/p>Here, \"k(x,\cdot)\" represents the univariate function obtained by setting the first input of \"k\" to be \"x\".<\/p>\n\n\n\n

The Ideal Weights<\/h2>\n\n\n\n

To design a quadrature rule, we have to set the nodes \"s = (s_1,\ldots,s_n) \in \Omega^n\" and weights \"w = (w_1,\ldots,w_n)\in\real^n\". Let\u2019s first assume that the nodes \"s\" are fixed, and talk about how to pick the weights \"w\".<\/p>\n\n\n\n

There\u2019s one choice of weights \"w^\star\" that we’ll called the ideal weights<\/em>. There (at least) are five equivalent ways of characterizing the ideal weights. We\u2019ll present all of them. As an exercise, you can try and convince yourself that these characterizations are equivalent, giving rise to the same weights.<\/p>\n\n\n\n

Interpretation 1: Exactness<\/h3>\n\n\n\n

A standard way of designing quadrature rules is to make them exact<\/em><\/a> (i.e., error-free) for some class of functions. For instance, many classical<\/a> quadrature<\/a> rules<\/a> are exact for polynomials of degree up to \"n-1\".<\/p>\n\n\n\n

For kernel quadrature, it makes sense to design the quadrature rule to be exact for the kernel function at the selected nodes. That is, we require

  <\/span>   <\/span>\"\[\hat{I}_{w_\star,s}[k(s_i,\cdot)]=I[k(s_i,\cdot)] \quad \text{for } i=1,2,\ldots,n.\]\"<\/p>Enforcing exactness gives us \"n\" linear equations for the \"n\" unknowns \"w^\star_1,\ldots,w^\star_n\":

  <\/span>   <\/span>\"\[\sum_{j=1}^n k(s_i,s_j)w^\star_j = \int_\Omega k(s_i,x) g(x)\,\mathrm{d}\mu(x) \quad \text{for }i=1,2,\ldots,n.\]\"<\/p>Under mild conditions, this system of linear equations is uniquely solvable, and the solution \"w^\star\in\real^n\" is the ideal weights.<\/p>\n\n\n\n

Interpretation 2: Interpolate and Integrate<\/h3>\n\n\n\n

Here\u2019s another very classical way of designing a quadrature rule. First, interpolate<\/a> the function values \"(s_i,f(s_i))\" at the nodes, obtaining an interpolant \"\hat{f}\". Then, obtain an approximation to the integral by integrating the interpolant:

  <\/span>   <\/span>\"\[\hat{I}_{w^\star,s}[f] \coloneqq \int_\Omega \hat{f}(x) g(x) \, \mathrm{d}\mu(x).\]\"<\/p>
In our context, the appropriate interpolation method is kernel interpolation<\/em>.1<\/a><\/sup>Kernel interpolation is also called Gaussian process regression<\/a> or kriging<\/a> though (confusingly) these terms can also refer to slightly different methods. It is the regularization-free limit of kernel ridge regression<\/a>.<\/span> The kernel interpolant is defined to be the minimum-norm function \"\hat{f}\" that interpolates the data:

  <\/span>   <\/span>\"\[\hat{f} = \argmin \{ \norm{h} : h(s_i) = f(s_i) \text{ for } i=1,\ldots,n\}.\]\"<\/p>Remarkably, this infinite-dimensional problem has a tractably computable solution. In fact, \"\hat{f}\" is the unique function of the form

  <\/span>   <\/span>\"\[\hat{f} = \sum_{i=1}^n \alpha_i k(\cdot,s_i)\]\"<\/p>that agrees with \"f\" on the points \"s_1,\ldots,s_n\".With a little algebra, you can show that the integral of \"\hat{f}\" is

  <\/span>   <\/span>\"\[I[\hat{f}] = \sum_{i=1}^n w^\star_i f(s_i),\]\"<\/p>where \"w^\star\" are the ideal weights.<\/p>\n\n\n\n

Interpretation 3: Minimizing the Worst-Case Error<\/h3>\n\n\n\n

Define the worst-case error of weights \"w\" and nodes \"s\" to be

  <\/span>   <\/span>\"\[\operatorname{Err}(w,s)=\sup_{\norm{f}\le 1}\left| I[f] - \hat{I}_{w,s}[f]\right|.\]\"<\/p>The quantity \"\operatorname{Err}(w,s)\" is the highest possible quadrature error for a function \"f\in\mathcal{H}\" of norm at most 1.<\/p>\n\n\n\n

Having defined the worst-case error, the ideal weights are precisely the weights that minimize this quantity

  <\/span>   <\/span>\"\[w^\star=\operatorname*{argmin}_{w\in\real^n}\operatorname{Err}(w,s).\]\"<\/p><\/p>\n\n\n\n

Interpretation 4: Minimizing the Mean-Square Error<\/a><\/h3>\n\n\n\n

The next two interpretations of the ideal weights will adopt a probabilistic framing. A Gaussian process<\/a><\/em> is a random function \"f\" such that \"f\"\u2019s values at any collection of points are (jointly<\/a>) Gaussian random variables<\/a>. We write \"f\sim \operatorname{GP}(0,k)\" for a mean-zero Gaussian process with covariance function<\/a> \"k\":

  <\/span>   <\/span>\"\[\Cov(f(x),f(y))=k(x,y)\quad \text{for every } x,y\in\Omega.\]\"<\/p>
Define the mean-square<\/a> quadrature error of weights \"w\" and nodes \"s\" to be

  <\/span>   <\/span>\"\[\operatorname{MSE}(w,s)\coloneqq \expect_{f\sim\operatorname{GP}(0,k)} \left( I[f] - \hat{I}_{w,s}[f] \right)^2.\]\"<\/p>The mean-square error reports the expected squared quadrature error over all functions \"f\" drawn from a Gaussian process with covariance function<\/a> \"k\".<\/p>\n\n\n\n

Pleasantly, the mean-square error<\/a> is equal ro the square of the worst-case error

  <\/span>   <\/span>\"\[\operatorname{MSE}(w,s)=\operatorname{Err}(w,s)^2.\]\"<\/p>As such, the ideal weights also minimize the mean-square error

  <\/span>   <\/span>\"\[w^\star=\operatorname*{argmin}_{w\in\real^n}\operatorname{MSE}(w,s).\]\"<\/p><\/p>\n\n\n\n

Interpretation 5: Conditional Expectation<\/a><\/h3>\n\n\n\n

For our last interpretation, again consider a Gaussian process \"h\sim \operatorname{GP}(0,k)\". The integral of this random function \"I[h]\" is a random variable. To numerically integrate a function \"f\", compute the expectation of \"I[h]\" conditional<\/a> on \"h\" agreeing with \"f\" at the quadrature nodes:

  <\/span>   <\/span>\"\[\hat{I}_{w^\star,s}[f]\coloneqq \expect_{h\sim\operatorname{GP}(0,k)}[I[h] \mid h(s_i)=f(s_i) \text{ for } i=1,\ldots,n].\]\"<\/p>One can show that this procedure yields the quadrature scheme with the ideal weights.<\/p>\n\n\n\n

Conclusion<\/h3>\n\n\n\n

We\u2019ve just seen five sensible ways of defining the ideal weights for quadrature in a general reproducing kernel Hilbert space. Remarkably, all five lead to exactly the same choice of weights. To me, these five equivalent characterizations give me more confidence that the ideal weights really are the \u201cright\u201d or \u201cnatural\u201d choice for kernel quadrature.<\/p>\n\n\n\n

That said, there are other reasonable requirements that we might want to impose on the weights. For instance, if \"\mu\" is a probability measure <\/a>and \"g\equiv 1\", it is reasonable to add an additional constraint that the weights \"w\" lie in the probability simplex<\/a>

  <\/span>   <\/span>\"\[w\in\Delta\coloneqq \left\{ p\in\real^n_+ : \sum_{i=1}^n p_i = 1\right\}.\]\"<\/p>With this additional stipulation, a quadrature rule can be interpreted as integrating \"f\" against a discrete probability measure \"\ \hat{\mu}=\sum_{i=1}^n w_i\delta_{s_i}\"; thus, in effect, quadrature amounts to approximating one probability measure \"\mu\" by another \"\hat{\mu}\". Additional constraints such as these can easily be imposed when using the optimization characterizations 3 and 4 of the ideal weights. See this paper<\/a> for details.<\/p>\n\n\n\n

What About the Nodes?<\/h2>\n\n\n\n

We\u2019ve spent a lot of time talking about how to pick the quadrature weights, but how should we pick the nodes \"s\in\Omega^n\"? To pick the nodes, it seems sensible to try and minimize the worst-case error \"\operatorname{Err}(w^\star,s)\" with the ideal weights \"w^\star\". For this purpose, we can use the following formula:

  <\/span>   <\/span>\"\[\operatorname{Err}(w^\star,s) = \norm{\int_\Omega (k(\cdot,x) - \hat{k}_s(\cdot,x)) g(x) \, \mathrm{d}\mu(x)}.\]\"<\/p>Here, \"\hat{k}_s\" is the Nystr\u00f6m approximation<\/a><\/em> to the kernel \"k\" induced by the nodes \"s\", defined to be

  <\/span>   <\/span>\"\[\hat{k}_s(x,y) = k(x,s) k(s,s)^{-1} k(s,y).\]\"<\/p>We have written \"k(s,s)\" for the kernel matrix<\/em> with \"ij\" entry \"k(s_i,s_j)\" and \"k(x,s)\" and \"k(s,y)\" for the row and column vectors with \"i\"th entry \"k(x,s_i)\" and \"k(s_i,y)\".<\/p>\n\n\n\n

I find the appearance of the Nystr\u00f6m approximation in this context to be surprising and delightful. Previously on this blog<\/a>, we’ve seen (column) Nystr\u00f6m approximation in the context of matrix low-rank approximation. Now, a continuum analog of the matrix Nystr\u00f6m approximation has appeared in the error formula for numerical integration.<\/p>\n\n\n\n

The appearance of the Nystr\u00f6m approximation in the kernel quadrature error \"\operatorname{Err}(w^\star,s)\" also suggests a strategy for picking the nodes.<\/p>\n\n\n\n

\n

Node selection strategy.<\/strong> We should pick the nodes \"s\" to make the Nystr\u00f6m approximation \"\hat{k}_s \approx k\" as accurate as possible.<\/p>\n<\/blockquote>\n\n\n\n

The closer \"\hat{k}_s\" is to \"k\", the smaller the function \"k(\cdot,x) - \hat{k}_s(\cdot,x)\" is and, thus, the smaller the error

  <\/span>   <\/span>\"\[\operatorname{Err}(w^\star,s) = \norm{\int_\Omega (k(\cdot,x) - \hat{k}_s(\cdot,x)) g(x) \, \mathrm{d}\mu(x)}.\]\"<\/p>Fortunately, we have randomized matrix algorithms for picking good nodes for matrix Nystr\u00f6m approximation such as randomly pivoted Cholesky<\/a>, ridge leverage score sampling<\/a>, and determinantal point process sampling<\/a>; maybe these matrix tools can be ported to the continuous kernel setting?<\/p>\n\n\n\n

Indeed, all three of these algorithms\u2014randomly pivoted Cholesky<\/a>, ridge leverage score sampling<\/a>, and determinantal<\/a> point<\/a> process<\/a> sampling<\/a>\u2014have been studied for kernel quadrature. The first of these algorithms, randomly pivoted Cholesky<\/a>, is the subject of our paper. We show that this simple, adaptive sampling method produces excellent nodes for kernel quadrature. Intuitively, randomly pivoted Cholesky is effective because it is repulsive<\/em>: After having picked nodes \"s_1,\ldots,s_i\", it has a high probability of placing the next node \"s_{i+1}\" far from the previously selected nodes.<\/p>\n\n\n\n

The following image shows 20 nodes selected by randomly pivoted Cholesky in a crescent-shaped region. The cyan\u2013pink shading denotes the probability distribution for picking the next node. We see that the center of the crescent does not have any nodes, and thus is most likely to receive a node during the next round of sampling.<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

In our paper<\/a>, we demonstrate\u2014empirically<\/a> and theoretically<\/a>\u2014that randomly pivoted Cholesky produces excellent nodes for quadrature. We also discuss<\/a> efficient rejection sampling<\/a> algorithms for sampling nodes with the randomly pivoted Cholesky distribution. Check out the paper<\/a> for details!<\/p>\n","protected":false},"excerpt":{"rendered":"

I’m excited to share that my paper Kernel quadrature with randomly pivoted Cholesky, joint with Elvira Moreno, has been accepted to NeurIPS 2023 as a spotlight. Today, I want to share with you a little about the kernel quadrature problem. To avoid this post getting too long, I’m going to write this post assuming familiarityRead more<\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[6,7],"tags":[],"class_list":["post-1700","post","type-post","status-publish","format-standard","hentry","category-expository","category-research"],"_links":{"self":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1700","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/comments?post=1700"}],"version-history":[{"count":9,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1700\/revisions"}],"predecessor-version":[{"id":1713,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/posts\/1700\/revisions\/1713"}],"wp:attachment":[{"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/media?parent=1700"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/categories?post=1700"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.ethanepperly.com\/index.php\/wp-json\/wp\/v2\/tags?post=1700"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}