EconBase
← Back to paper

Minimax Optimal Kernel Operator Learning via Multilevel Training

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

50,324 characters · 14 sections · 61 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Minimax Optimal Kernel Operator Learning via Multilevel Training

abstractLearning mappings between infinite dimensional function spaces has achieved empirical success in many disciplines of machine learning, including generative modeling, machine learning solving Partial Difference Equations , functional data analysis, causal inference, and multi-agent reinforcement learning. In this paper, we study the statistical limit of learning a Hilbert-Schmidt operator between two infinite-dimensional Sobolev reproducing kernel Hilbert spaces. We establish the information-theoretic lower bound in terms of the Sobolev Hilbert-Schmidt norm and show that a regularization that learns the spectral components below the bias contour and ignores the ones that above the variance contour can achieve the optimal learning rate. At the same time, the spectral components between the bias and variance contours give us the flexibility in designing computationally feasible machine learning algorithms. Based on this observation, we develop a multilevel kernel operator learning algorithm that is optimal when learning linear operators between infinite-dimensional function spaces.

Introduction

The supervised learning of operators between two infinite-dimensional spaces has attracted attention in many machine learning applications, such as scientific computation lu2019deeponet,li2020fourier,de2021convergence,li2018organic,li2021principles, functional data analysis crambes2013asymptotics,hormann2015note,wang2020functional, learning mean-field games guo2019learning,wang2020breaking, conditional probability regression song2009hilbert,song2013kernel,muandet2017kernel and econometrics singh2019kernel,muandet2020dual,dikkala2020minimax. Despite the empirical success of operator learning, the statistical limit of learning an infinite-dimensional operator is poorly studied. In this paper, we study the problem of learning Hilbert Schmidt operators between infinite-dimensional Sobolev reproducing kernel Hilbert spaces $\mathcal{H}_K^\beta$ and $\mathcal{H}_L^\gamma$ with given kernels $k$ and $l$ respectively and $\beta,\gamma\in[0,1)$ adams2003sobolev,christmann2008support,fischer2020sobolev. Our goal is to derive the optimal sample complexity to learn the linear operator, i.e. how much data is required to achieve a certain performance level.

We first establish an information-theoretic lower bound for learning a Hilbert-Schmidt operator between Sobolev spaces respect to a general Sobolev norm. Our information-theoretic lower bound indicates that the optimal learning rate is determined by the minimum of two polynomial rates: one is purely decided by the input Sobolev reproducing kernel Hilbert space and its evaluating norm, while the other one is purely determined by the output space along with its evaluating norm. The rate is novel in the sense that all existing results fischer2020sobolev,li2022optimal,de2021convergence only establish rates that depend on the parameter of input space. The reason is all previous works talwai2022sobolev,li2022optimal,de2021convergence only consider the case of the output space as a subspace of a trace bounded reproducing kernel Hilbert space but not a general Sobolev space. We refer to (ref) for detail comparisons.

To design a learning algorithm for approximating an infinite-dimensional operator, we need to learn a finite-dimensional restriction instead of the whole operator, as the latter would result in infinite variance. The finite-dimensional selection leads to bias error but decreases the variance. A natural task is then to study the shape of regularization that can lead to the optimal bias-variance trade-off and achieve the optimal learning rate. In this paper, we consider the bias and variance contour at the scale of optimal learning. Once the regularization enables one to learn all the spectral parts above the bias contour and below the variance contour, the learning is optimal. Finally, utilizing the region between the bias contour and variance contour, we developed a multilevel training algorithm lye2021multi,li2021multilevel which first learns the mapping on low frequency and then successively fine-tunes the machine learning models to fit the high-frequency output. The intuition of our algorithm aligns with the original motivation of multilevel Monte Carlo giles2008multilevel,giles2015multilevel: we use the next level to reduce bias while keeping the variance at the same scale. We demonstrate that such a multilevel algorithm can achieve an optimal non-parametric rate for linear operator learning.

Related Work

\paragraph{Machine Learning Based PDE Solver} Solving partial differential equations (PDEs) plays a prominent role in many scientific and engineering discipline, such as physics, chemistry, operation management, macro-economy, etc. The recent deep learning breakthrough has drawn attention to solving PDEs via machine learning methods raissi2019physics,han2018solving,sirignano2018dgm,yu2018deep,khoo2019solving,chen2021solving. The statistical power and computational cost of these problem is well-studied by recent papers lu2021machine,lu2022sobolev,nickl2020convergence,nickl2020polynomial. This paper focuses on operator learning chen1995universal,long2018pde,long2019pde,feliu2020meta,khoo2021solving,lu2019deeponet,li2020fourier,kovachki2021neural,stepaniants2021learning, i.e. learning a map between two infinite dimensional function spaces. For example, one can learn a PDE solver that maps from the boundary condition to the solution or an inverse problem that maps from the boundary measurement to the coefficient field. In terms of the mathematical foundation of operator learning, liu2022deep considers the learning rate of non-parametric operator learning. However, non-parametric functional data analysis often suffers from slower-than-polynomial convergence rates mas2012lower, due to the small ball probability problem for the probability distributions in infinite dimensional spaces delaigle2010defining. The most relevant works are lin2011fast,reimherr2015functional,de2021convergence, which consider the rates for learning a linear operator. For the comparison between our work and de2021convergence, see (ref).

\paragraph{Learning with kernel.} Supervised least square regression in RKHS and its generalization capability have been thoroughly studied caponnetto2007optimal,smale2007learning,de2005learning,rosasco2010learning,mendelson2010regularization. The minimax optimality with respect to the Sobolev norm has been discussed recently in fischer2020sobolev,liu2020estimation,lu2022sobolev. Our paper is highly related to recent works schuster2020kernel,mollenhauer2020nonparametric,talwai2022sobolev,li2022optimal on identifying the Sobolev norm learning rate for the kernel mean embeddingsong2009hilbert,song2013kernel,muandet2017kernel,park2020measure,singh2020kernels, which can also formulated as learning an operator. The difference of our work see (ref).

\paragraph{Multilevel Monte Carlo} By combining biased estimators with multiple stepsizes, multilevel Monte Carlo (MLMC) giles2008multilevel,giles2015multilevel dramatically improves the rate of convergence and achieves in many settings the canonical square root convergence rate associated with unbiased Monte Carlo rhee2015unbiased,blanchet2015unbiased. Multilevel Monte Carlo can also be used for random variable with infinite variance blanchet2016malliavin,chen2020unbiased. To the best of our knowledge, this is the first paper that provides optimal sample complexity for multilevel Monte Carlo type algorithm for infinite variance problems in the non-parametric regime. Very recently, lye2021multi,li2021multilevel developed a multilevel machine learning Monte Carlo algorithm (ML2MC) / multilevel fine-tuning algorithm for learning solution maps, by first learning the map on coarsest grid and then successively fine-tuning the network on samples generated at finer grids. The authors also showed that, following the telescoping in MLMC, the multilevel training procedure can reduce the generalization error without spending more time on generating training samples. schafer2021sparse,boulle2022learning consider such multi-scale algorithm for learning Green's function. However, the statistical power of such algorithm is still under investigation. Another difference with boulle2022learning is that we consider the Green function in $H^{-1}$ norm rather than the $\ell_1$ norm used in boulle2022learning. In this paper, we qualify a specific setting where this multilevel procedure can and is necessary to achieve the minimax optimal learning rate.

Contribution

itemize• We derive a novel information-theoretic lower bound of learning a linear operator between two infinite-dimensional Sobolev reproducing kernel Hilbert spaces. The optimal learning rate is a minimum of two polynomial rates, one only dependent on the parameters of the input space while the other only on the parameters of the output space. The first rate aligns with the previous works li2022optimal, while the second lower bound is novel to the literature. • We study the shape of regularization that can lead to the optimal learning rate. One should learn all the spectral parts under the bias contour at the level of the optimal learning rate but not the spectral parts above the variance contour at the level of learning rate. This enables the estimator to enjoy an optimal balance of bias-variance. • We qualify a specific setting where a multilevel training procedure lye2021multi,li2021multilevel is necessary and capable of achieving a minimax optimal learning rate for learning a linear operator. We achieve the optimal learning rate via $O(\ln \ln n)$ ensemble of ridge regression models. This is different from finite-dimensional operator learning where a single level estimator can be optimal.

Problem Formulation

Preliminary

Let $P_K$ be a distribution over the input space $\mathcal{H}_K$ and define covariance operator $\mathcal{C}_{KK}= \mathbb{E}_{u\sim P_K} u\otimes u$. Consider its spectral decomposition $\mathcal{C}_{KK} = \sum_{i=1}^{+\infty} \mu_i^2 e_i \otimes e_i$, where $\{\mu_i^{\frac{1}{2}}e_i\}_{i=1}^{+\infty}$ is an orthogonal eigenbasis and $\{\mu_i\}$ is the corresponding eigenvalues of $\mathcal{C}_{KK}$ (here the $g\otimes h$ is an operator defined as $g\otimes h=gh^\ast:f\rightarrow \left\langle f,h\right\rangle g$). In the typical machine learning applications, the test distribution is the same as the training distribution, so we can assume that $\mathcal{H}_K = \left\{ \sum_{i} a_i \mu_i^{\frac{1}{2}}e_i: \{a_i\}_{i=1}^{\infty}\in\ell_2 \right\}$ without loss of generality. Note that this automatically holds in the context of learning the conditional mean embedding (CME) fischer2020sobolev,talwai2022sobolev,li2022optimal.

Following christmann2008support,fischer2020sobolev, we define the interpolation Sobolev space $\mathcal{H}_K^\beta=\left\{f=\sum_i a_i(\mu_i^{\frac{\beta}{2}}e_i):\{a_i\}_{i=1}^\infty\in l^2\right\}$ for any $\beta>0$, equipped with Sobolev norm defined by the inner product $\left<\sum_i a_i(\mu_i^{\beta/2}e_i),\sum_i b_i(\mu_i^{\beta/2}e_i)\right>_{\mathcal{H}_K^\beta}=\sum_i a_i b_i$. For the output space, we fix a user-specified distribution $Q_L$ and a reproducing Kernel Hilbert Space. We can similarly define the covariance operator $\mathcal{C}_{Q_L}$ and the Sobolev space $\mathcal{H}_L^{\gamma}$. Natural choices of $Q_L$ include some distribution on kernel functions $\left\{\ell(y,\cdot):y\in Y\right\}$ of $\mathcal{H}_L$ induced by some distribution $Q_L$ on $Y$, so that $\mathcal{C}_{Q_L}$ is a kernel integral operator with respect to $Q_L$ and $\mathcal{H}_L^{\gamma}$ is an interpolation space between $\mathcal{H}_L$ and $\mathcal{L}^2(Q_L)$; see (ref) for a specific choice of $Q_L$ and its practical implications.

Following li2022optimal, in this paper we consider the Hilbert-Schmidt norm between two Sobolev Spaces for all the operators, which is defined as following.

definition[$(\beta,\gamma)$-norm] Let $T:\mathcal{H}_K\mapsto\mathcal{H}_L$ be a possibly unbounded linear operator. $I_{1,\beta,P_{K}}:\mathcal{H}_K\mapsto\mathcal{H}_K^{\beta}, \beta\in(0,1)$ is the canonical embedding mapping that takes $u\in\mathcal{H}_K$ to the same element $u$ in the larger space $\mathcal{H}_K^{\beta}$, and $I_{1,\gamma,Q_L}:\mathcal{H}_L\mapsto\mathcal{H}_L^{\gamma},\gamma\in(0,1)$ is similarly defined. Then the $(\beta,\gamma)$-norm of $T$ is defined as \begin{equation} \notag \left\VertT\right\Vert_{\beta,\gamma} = \left\Vert(I_{1,\gamma,Q_{L}}^\ast)^{\dagger}\circ T \circ I_{1,\beta,P_{K}}^\ast\right\Vert_{\mathrm{HS}\left(\mathcal{H}_K^{\beta},\mathcal{H}_{L}^{\gamma}\right)}= \left\Vert\mathcal{C}_{Q_L}^{-(1-\gamma)/2}\circ T\circ \mathcal{C}_{KK}^{(1-\beta)/2}\right\Vert_{\mathrm{HS}\left(\mathcal{H}_K,\mathcal{H}_{L}\right)}, \end{equation} where we omit the dependence of $\left\Vert\cdot\right\Vert_{\beta,\gamma}$ on $P_{K}$ and $Q_L$ since it will always be clear from context.

Problem Formulation

We consider the problem of learning an unknown linear operator $\mathcal{A}_0:\mathcal{H}_K\mapsto\mathcal{H}_L$ between two reproducing kernel Hilbert spaces corresponding to kernl $k$ and $l$ respectively. We are given $N$ noisy data pairs $(u_i,v_i),1\leqslant i\leqslant N$ related by

equation[equation omitted — 81 chars of source]

where $u_i \mathrel{\overset{\mathrm{i.i.d.}}{\scalebox{2.33}[1]{$\sim$}}} P_{K}$ for some unknown distribution $P_{K}$ and ${\epsilon}_i$ is the noise drawn from some distribution with zero mean that may depend on $u_i$. We use $P_{KL}$ for the joint distribution of $(u_i,v_i)$. Denote $\mathcal{C}_{KK} = \mathbb{E}_{u\sim P_K}u\otimes u$, $\mathcal{C}_{KL}=\mathbb{E}_{(u,v)\sim P_{KL}}u\otimes v$ and its adjoint $\mathcal{C}_{LK}=\mathcal{C}_{KL}^\ast$ be uncentered cross-covariance operators associated with $P_{KL}$. Then we can reformulate the ground turth operator as $\mathcal{A}_0=\mathcal{C}_{LK}\mathcal{C}_{KK}^{\dagger}$, where $\dagger$ is the pseudo-inverse talwai2022sobolev,li2022optimal. With the goal of understanding the relative difficulty of learning different types of linear operators, we investigate the sample efficiency of learning $\mathcal{A}_0$ under certain source assumptions imposed on the data model (ref). Source condition caponnetto2007optimal,mendelson2010regularization,steinwart2009optimal,rosasco2010learning,fischer2020sobolev assumes that the learning target lies in a parameterized function class and study the learning rate for different problems with different hardness. Specifically, the source condition assume that the learning target is bounded in certain Sobolev norm. In this paper, we consider learning an operator with bounded $(\beta,\gamma)$-norm, which is the Hilbert-Schmidt norm that maps from $\mathcal{H}_K^\beta$ to $\mathcal{H}_L^\gamma$. We consider the generalization error/convergence rate under another $(\beta',\gamma')$-norm as in fischer2020sobolev,lu2022sobolev,talwai2022sobolev,li2022optimal.

remarkAlthough recent works have considered similar problems in the context of conditional mean embedding talwai2022sobolev,li2022optimal and functional data analysis de2021convergence. In all these papers, the output space is a trace bounded reproducing kernel Hilbert space de2021convergence rather than the general parameterized Sobolev space in our paper.

We then list all the assumptions imposed on the underlying kernel for our theoretical results. We follow the standard capacity assumptions and embedding properties used in kernel regression fischer2020sobolev,talwai2022sobolev,li2022optimal.

assumption[Capacity Condition of the Covariance] The eigenvalues $\{\mu_i\}_{i\geqslant 1}$ of the covariance operator $\mathcal{C}_{KK} = \mathbb{E}_{u\sim P_{K}}u\otimes u$ satisfies $\mu_i\propto i^{-\frac{1}{p}}$ for some $p\in(0,1)$. Similarly, the eigenvalues $\{\rho_i\}_{i\geqslant 1}$ of the covariance operator $\mathcal{C}_{Q_L} = \mathbb{E}_{v\sim Q_L}v\otimes v$ satisfies $\rho_i\propto i^{-\frac{1}{q}}$ for some $q\in(0,1)$.
assumption[$\ell_\infty$ Embedding Property of the Input RKHS] There exists a smallest $\alpha\in(0,1)$ such that $\left\Vert\left(I_{1,\alpha,P_{K}}^\ast\right)^{\dagger}f\right\Vert_{\mathcal{H}_K^{\alpha}} \leqslant A_1$ a.s. under $P_{K}$ for some $A_1<+\infty$.
assumption[$\ell_\infty$ Embedding Property of the Output RKHS] There exists $A_2 < +\infty$ such that $\left\Vert\mathcal{A}_0 u\right\Vert_{\mathcal{H}_{L}} \leqslant A_2$ holds for all $u \in \mathrm{supp}(P_K)$, except from a $P_K$-null set.
assumption[Moment Condition] There exists an operator $V:\mathcal{H}_{L}\mapsto\mathcal{H}_{L}$ with $\operatorname{tr}\left(V\right) \leqslant \sigma^2$ such that for every $u \in \mathrm{supp}(P_K)$, we have $\mathbb{E}_{v\sim P_{KL}(\cdot\mid u)}\left[ \left((v-\mathcal{A}_0 u) \otimes (v-\mathcal{A}_0 u)\right)^k \right] \preceq \frac{1}{2}(2k)! R^{2k-2}V$ holds for all $k\ge 2$.
assumption[Source Condition] $\mathcal{A}_0$ is bounded under $(\beta,\gamma)$-norm i.e. $\left\Vert\mathcal{A}_0\right\Vert_{\beta,\gamma} \leqslant B$ for some $B\in(0,+\infty)$.

Examples

In this section, we will introduce two examples of our theory. The first one is about learning a differential operator, for example inferring an advection-diffusion model portone2022bayesian from observations or predicting the future long2018pde,lu2019deeponet,li2020fourier,feliu2020meta,huang2021meta. The second example is about learning conditional mean embedding song2009hilbert,song2013kernel,muandet2017kernel, which represents a conditional distribution as an RKHS element. Thus conditional distribution regression can be reduced to a kernel operator learning. Our theory can also be used for linear inverse problem such as radial electrical impedance tomography (EIT) mueller2012linear and the severely ill-posed inverse boundary problem for the Helmholtz equation with unknown wave-number parameter agapiou2014bayesian. For detailed discussion, we refer to de2021convergence

example[Learning differential operators] Suppose that the ground-truth operator $\mathcal{A}_0 = \Delta^{t}$ where $\Delta$ is the Laplacian and $t \in\mathbb{Z}$. Let $\mathcal{H}_K = \mathcal{H}^{m+2t}([0,1])$ be the Sobolev space with smoothness $m+2t$ on $[0,1]$ and $\mathcal{H}_L = \mathcal{H}^{m}([0,1])$, then $\mathcal{A}_0$ is a bounded operator from $\mathcal{H}_K$ to $\mathcal{H}_L$ which corresponds to the $\beta=\gamma=1$ case. However, we will see below that we can obtain a better characterization of the learning error using our theory. Consider for example that the input has mean zero and the Mat{\'e}rn-type covariance operator $C_{KK} = \sigma^2\left(-\Delta+\tau^2 I\right)^{-s}$. Its eigenvalues satisfy $\mu_n \propto n^{-2s}$. On the other hand, we choose $Q_L$ to be a distribution supported on $\left\{\ell(y,\cdot):y\in[0,1]\right\}$ induced by a uniform distribution on $[0,1]$, where $\ell$ is the kernel function of $\mathcal{H}_L$. Then $\mathcal{C}_{Q_L}$ is essentially the kernel integral operator on $\mathcal{H}_L$ w.r.t. the uniform distribution, and its eigenvalues are $\rho_n \propto n^{-2m}$. The assumption $||\mathcal{A}_0||_{\beta,\gamma} < +\infty$ is satisfied if and only if $(1-\gamma)m < (1-\beta)s-\frac{1}{2} \Rightarrow \gamma > 1 - \frac{2(1-\beta)s-1}{2m}$.
example[Conditional mean embedding] Suppose that we would like to learn the conditional distribution $P(y\mid x)$ from a data set $\{(x_i,y_i):1\leqslant i\leqslant N\}\subset X\times Y$ where $x_i \mathrel{\overset{\mathrm{i.i.d.}}{\scalebox{2.33}[1]{$\sim$}}} P_X$. Let $\mathcal{H}_K$ and $\mathcal{H}_L$ be two RKHSs on $X$ and $Y$ respectively, with measurable kernel $k(\cdot,\cdot)$ and $\ell(\cdot,\cdot)$. Then we can define a conditional mean embedding (CME) operator $C_{Y\mid X}$ that satisfies $$ C_{Y\mid X} k(x,\cdot) = \mathbb{E}_{Y\mid x}\ell(Y,\cdot) =: \mu_{Y\mid x}, \text{ and } \mathbb{E}_{Y\mid x} \left[ g(Y) \right] = \left\langle g,\mu_{Y\mid x}\right\rangle \forall x\in X. $$ We choose $\mathcal{A}_0 = C_{Y\mid X}$. In this case, $\mathcal{C}_{KK} = \mathbb{E}_{P_X} k(X,\cdot)\otimes k(X,\cdot)$. (ref) states that $\sup_{x\in X}k^{\alpha}(x,x) = A_1$, while (ref) is equivalent to $\sup_{x\in X}\left\Vert\mu_{Y\mid x}\right\Vert \leqslant A_2$ (for simplicity we only focus on the case $\zeta = 1$). According to (ref), we assume that $\left\VertC_{Y\mid X}\right\Vert_{\beta,\gamma}\leqslant B$. The mis-specified setting where $\beta<1$ has been studied in previous work fischer2020sobolev,talwai2022sobolev,li2022optimal. However, they only consider the case $\gamma = 1$. Our results also cover the case $\gamma<1$, which allows us to obtain theoretical guarantee for computing conditional expectation of the larger function class $\mathcal{H}_L^{\gamma}$. \begin{comment} To see this, consider the well-specified setting $\beta>\alpha$. For any $g \in \mathcal{H}_L^{\gamma'}$, we have \begin{equation} \notag \begin{aligned} \left| \left\langle g,\hat{\mu}_{Y\mid x}\right\rangle_{\mathcal{H}_L}-\mathbb{E}_{Y\mid x} g(Y)| &= \left| \left\langle g,\mu_{Y\mid x}-\hat{\mu}_{Y\mid x}\right\rangle_{\mathcal{H}_L}| \\ &\leqslant \left||C_{LL}^{\frac{1-\gamma'}{2}}g||_{\mathcal{H}_L}\cdot \left||C_{YY}^{-\frac{1-\gamma'}{2}}\left(\mu_{Y\mid x}-\hat{\mu}_{Y\mid x}\right)||_{\mathcal{H}_L} \\ &\leqslant \sup_{x\in X}k^{\beta'}(x,x)\cdot\left||C_{YY}^{\frac{1-\gamma'}{2}}g||_{\mathcal{H}_L}\cdot \left||\hat{\mathcal{A}}-\mathcal{A}_0||_{\beta',\gamma'} \end{aligned} \end{equation} when the conditional expectation $\mathbb{E}_{Y\mid x}g(Y)$ is finite. When $\beta'\in(\alpha,\beta)$ we also have $\sup_{x\in X}k^{\beta'}(x,x) < +\infty$. As a result, the above gives a uniform convergence rate of the estimation error of conditional expectation in $\mathcal{H}_L^{\gamma'}$. \end{comment}

Information Theoretic Lower Bound

In this section, we provide an information-theoretic lower bound for the convergence rate of the operator learning problem formulated in (ref).

theoremSuppose that $\mathcal{H}_K$ and $\mathcal{H}_L$ are two Hilbert spaces, $P_K$ and $Q_L$ are probability distributions on $\mathcal{H}_K$ and $\mathcal{H}_L$ respectively such that (ref) hold. Then for any estimator $\mathcal{L}:\left(\mathcal{H}_K\times \mathcal{H}_L\right)^{\otimes N}\mapsto \mathrm{HS}\left(\mathcal{H}_K^{\beta},\mathcal{H}_L^{\gamma}\right)$, there exists a linear operator $\mathcal{A}_0$ and a joint data distribution $P_{KL}$ with marginal distribution $P_K$ on $\mathcal{H}_K$ satisfying (ref), such that with probability $\geqslant 0.99$ over $(u_i,v_i) \mathrel{\overset{\mathrm{i.i.d.}}{\scalebox{2.33}[1]{$\sim$}}} P_{KL}$ we have \begin{equation} \notag \left\Vert\mathcal{L}\left(\{(u_i,v_i)\}_{i=1}^N\right)-\mathcal{A}_0\right\Vert_{\beta',\gamma'}^2 \gtrsim N^{-\min\left\{\frac{\max\{\alpha,\beta\}-\beta'}{\max\{\alpha,\beta\}+p},\frac{\gamma'-\gamma}{1-\gamma}\right\}}. \end{equation}
remarkOur lower bound is composed of a minimum of two parts. The first rate $N^{-\frac{\max\{\alpha,\beta\}-\beta'}{\max\{\alpha,\beta\}+p}}$ is the minimax optimal Sobolev learning rate for kernel regression fischer2020sobolev,talwai2022sobolev,li2022optimal,lu2022sobolev and is fully determined by the parameter of the input Sobolev reproducing kernel Hilbert space. Our second rate $N^{-\frac{\gamma'-\gamma}{1-\gamma}}$ is novel to the literature. This bound shows that how the infinite dimensional problem is different from finite dimensional regression problem and is fully determined by the parameter of the output Sobolev reproducing kernel Hilbert space. Our lower bound shows that the hardness of learning a linear operator is determined by the harder part between the input and output spaces. We will explain why the lower bound has such structure in (ref) and (ref).

On the Shape of Regularization

In this section, we aim to understand the shape of regularization so that the constructed estimator $\hat{\mathcal{A}}$ based on $N$ i.i.d. data $\{(u_i,v_i)\}_{i=1}^n\sim P_{KL}^{\otimes n}$ for $1\leqslant i\leqslant N$ enjoys an optimal learning rate.

Compared with existing approaches where a regularized least-squares estimator can achieve statistical optimality fischer2020sobolev,talwai2022sobolev,li2022optimal,de2021convergence under $(\beta,1)$-norm, we study the learning rate under the $(\beta,\gamma)$-norm ($\beta'\in(0,\beta),\gamma'\in(\gamma,1)$) which is defined in (ref) as $\left\Vert\hat{\mathcal{A}}-\mathcal{A}_0\right\Vert_{\beta',\gamma'}=\left\Vert\mathcal{C}_{Q_L}^{-\frac{1-\gamma'}{2}}\left(\hat{\mathcal{A}}-\mathcal{A}_0\right)\mathcal{C}_{KK}^{\frac{1-\beta'}{2}}\right\Vert_{\mathrm{HS}(\mathcal{H}_K,\mathcal{H}_L)}$. The norm of the additional $\mathcal{C}_{Q_L}^{-\frac{1-\gamma'}{2}}$ term is unbounded which make our setting harder than the convergence in $(\beta,1)$-norm in existing works. Since $\mathcal{C}_{Q_L}^{-\frac{1-\gamma'}{2}}$ is bounded when restricted to the finite-dimensional space $\mathrm{span}\left(\rho_i^{\frac{1}{2}}f_i:1\leqslant i\leqslant n\right)$, we should also include another bias-variance trade-off via regularizing in the output shape. As a result, we are interested in answering the following question

center[center omitted — 194 chars of source]

To answer this question, we investigate the problem in the spectral space, i.e. considering the spectral representation of operator $\mathcal{A}_0=\sum_{i,j=1}^{+\infty}a_{ij}\mu_i^{\frac{\beta}{2}}e_i\otimes \rho_j^{1-\frac{\gamma}{2}}f_j$. The problem of estimating $\mathcal{A}$ then reduces to learning the coefficients “matrix" $(a_{ij})_{i,j=1}^\infty$. The source condition (ref) enforces $\sum_{i,j=1}^\infty a_{ij}^2\le B$. We show in (ref) that regularizing the basis $e_i\otimes f_j$ will introduce a bias of order $\left\Verta_{ij}\mu_i^{\frac{\beta}{2}}e_i\otimes \rho_j^{1-\frac{\gamma}{2}}f_j\right\Vert_{\beta',\gamma'}^2 = a_{ij}^2\mu_i^{\beta-\beta'}\rho_j^{\gamma'-\gamma} \propto i^{-\frac{\beta-\beta'}{p}}j^{-\frac{\gamma'-\gamma}{q}}$ under the $(\beta',\gamma')$-norm. On the other hand, when $\alpha\leqslant\beta+p$, we show in (ref) that the variance of learning $(i,j)$ from noisy data scales as $\frac{1}{N}\mu_i^{-\beta'}\rho_j^{-(1-\gamma')} \propto \frac{1}{N}i^{\frac{\beta'}{p}}j^{\frac{1-\gamma'}{q}}$. Since the variance would accumulate for a fixed $j$, learning $(i,j)$ for $i\leqslant i_{\max}$ results in a variance of $\propto \frac{1}{N}i_{\max}^{\frac{\beta'+p}{p}}j^{\frac{1-\gamma'}{q}}$. (Similar analysis can be carried out for the ${\alpha>\beta+p}$ case as well, but the variance now scales as $\frac{1}{N}i_{\max}^{\frac{\beta'+\alpha-\beta}{p}}j^{\frac{1-\gamma'}{q}}$; see (ref) for detailed derivations.) In summary, we need to make bias-variance trade off in the $(i,j)-$plane, i.e. decide whether we should learn or regularize over the basis $e_i\otimes f_j$.

commentLet $y_N = N^{\frac{q}{1-\gamma'}\max\left\{1-\frac{\beta-\beta'}{\max\{\alpha,\beta+p\}},\frac{1-\gamma'}{1-\gamma} \right\}}$. For all $1\leqslant j\leqslant y_N$, we define \begin{equation} \lambda_j = \max\left\{ \left( N^{\max\left\{ 1-\frac{\beta-\beta'}{\max\{\alpha,\beta+p\}},\frac{1-\gamma'}{1-\gamma} \right\}} j^{-\frac{1-\gamma'}{q}}\right)^{-\frac{1}{\beta'+\max\{\alpha-\beta,p\}}}, c_0 \left(\frac{N}{\log N}\right)^{-\frac{1}{\alpha}}\right\} \end{equation} where $c_0$ is a small constant. Given data $\{(u_i,v_i):1\leqslant i\leqslant N\}$, we can compute the following empirical quantities \begin{equation} \notag \hat{\mathcal{C}}_{KK} = \frac{1}{N}\sum_{i=1}^N u_i \otimes u_i,\quad \hat{\mathcal{C}}_{XY} = \frac{1}{N}\sum_{i=1}^N v_i \otimes u_i. \end{equation} Then our estimator is defined as \begin{equation} \hat{\mathcal{A}} = \sum_{j=1}^{y_N} \left( \rho_j^{\frac{1}{2}}f_j \otimes \rho_j^{\frac{1}{2}} f_j\right) \hat{\mathcal{C}}_{XY} \left(\hat{\mathcal{C}}_{KK}+\lambda_j I\right)^{-1}. \end{equation} The main idea behind (ref) and (ref) will be elaborated in (ref). \begin{remark}[Computational aspects of our estimator] Computing the quantities $\hat{\mathcal{C}}_{KK}$ and $\hat{C}_{XY}$ only requires the knowledge of the given data and the underlying Hilbert space $\mathcal{H}_K$. On the other hand, since the distribution $Q_Y$ on $\mathcal{H}_L$ is user specified, we can assume that its eigenbasis $\left\{\rho_j^{\frac{1}{2}}f_j:j\geqslant 1\right\}$ is known. \end{remark}

Intuitive Explanation

{In this section, we provide an intuitive explanation of our lower bound (Theorem (ref)). As the previous paragraph explains, learning an operator is equivalent to learning an "infinite" size matrix with larger variance and smaller bias on the right upper corner. The core proof of this paper is considering proper bias-variance trade-off. Suppose one wants to construct an estimator with $N^{-\theta}$ learning rate. They need to learn every spectral component below the bias counter at the level of $N^{-\theta}$. Otherwise, the bias itself will become larger than $N^{-\theta}$. At the same time, they still need not learn any spectral component above the variance counter at the level of $N^{-\theta}$. Otherwise, the variance itself will become larger than $N^{-\theta}$. Thus, to enable an feasible estimator achieves learning rate at $N^{-\theta}$, the variance counter at the level of $N^{-\theta}$ should always be above the bias counter at the level of $N^{-\theta}$ (Figure (ref)). Depending on the hyperparameters, there are two different ways to achieve this goal, as shown in Figure (ref). Each situation is mapped to the rate depending only on the input space, and the rate depends only on the output space. In Section (ref), we demonstrated how a multilevel algorithm could be used to satisfy this requirement.}

Regularization via variance contour

The underlying idea of regularization is that some components are intrinsically hard to learn due to large variance; these components are then neglected by adding regularization and are counted as bias. The remaining components are easy to learn due to controllable variance. This intuition works well when the estimation error results from the noise of the data and is well-studied in a line of works fischer2020sobolev,de2021convergence,talwai2022sobolev,li2022optimal. This idea still works in our setting, but we need to re-evaluate the bias and variance of each component. Since we work with the Hilbert-Schmidt norm, this can be done in a coordinate-wise manner, meaning that we can look at each $a_{ij}$ separately and decide whether to neglect it (contribute to bias) or to learn it from data (contribute to variance).

figure*[figure* omitted — 732 chars of source]

Since the variance term measures the hardness of learning, we naturally introduce the notion of variance contour, which is a curve on the $\mathbb{R}_{+}^2$ plane on which all points induce the same order of variance (here we work with real coordinates for convenience, although we only care about integer points). Formally, we fix an arbitrary constant $C > 0$ and define

equation[equation omitted — 184 chars of source]

A reasonable regularization scheme is then to learn all coordinates $(i,j)\in\mathbb{Z}_+^2$ below the curve $\ell_{C,\mathtt{var}}$ and `regularize out' the remaining coordinates that are difficult to learn due to large variances. This can gives us the estimator with smallest estimator at give variance level. This observation motivates us to construct our estimator as

equation[equation omitted — 227 chars of source]

where {$\hat C_{LK} = \frac{1}{N}\sum_{i=1}^N v_i \otimes u_i$,} $\lambda_j (1\leqslant j\leqslant y_N= C^{\frac{q}{1-\gamma'}})$ are the regularization coefficients imposed on different dimensions of the output space. According to (ref) and noting that $\mu_i \propto i^{-\frac{1}{p}}$, we define

equation[equation omitted — 300 chars of source]

with $C = N^{\max\left\{ 1-\frac{\beta-\beta'}{\max\{\alpha,\beta+p\}},\frac{1-\gamma'}{1-\gamma} \right\}}$ in (ref). The additional $N^{-\frac{1}{\alpha}}$ term in (ref) is needed for controlling the error of approximating $\mathcal{C}_{KK}$ via $\hat{\mathcal{C}}_{KK}$ (cf. (ref)) which is standard in the Sobolev learning literature fischer2020sobolev,talwai2022sobolev,lu2022sobolev. The following theorem describes the convergence rate of our estimator defined by (ref).

theoremConsider the estimator $\hat{\mathcal{A}}$ defined by (ref). Suppose that (ref) hold, then there exists a universal constant $C$ such that with probability $\geqslant 1-e^{-\tau}$, we have \begin{equation} \notag \left\Vert\hat{\mathcal{A}}-\mathcal{A}_0\right\Vert_{\beta',\gamma'}^2 \leqslant C \tau^2 \left(\frac{N}{\log N}\right)^{-\min\left\{\frac{\beta-\beta'}{\max\{\alpha,\beta+p\}},\frac{\gamma'-\gamma}{1-\gamma}\right\}} \log^2 N. \end{equation}
remarkCompared with (ref), our upper bound is optimal up to logarithmic factors when $\alpha\leqslant\beta$. The optimal learning rate in the $\alpha > \beta$ regime is an outstanding problem for decades, even without the additional problem-dependent parameters $\gamma,\gamma'$ (see e.g. the discussions following fischer2020sobolev). In this paper, we do not address this problem either.
figure*[figure* omitted — 1,041 chars of source]
comment\begin{proof-sketch} For simplicity, we assume $\alpha\leqslant\beta+p$ here. The case $\alpha>\beta+p$ is similar and deferred to the appendix. Consider the bias-variance decomposition $||\hat{\mathcal{A}}-\mathcal{A}_0||_{\beta',\gamma'} \leqslant ||\hat{\mathcal{A}}-\mathcal{A}_{\lambda}||_{\beta',\gamma'}+||\mathcal{A}_{\lambda}-\mathcal{A}_0||_{\beta',\gamma'}$, where $\mathcal{A}_{\lambda} = \sum_{j=1}^{y_N} \left( \rho_j^{\frac{1}{2}}f_j \otimes \rho_j^{\frac{1}{2}} f_j\right) \mathcal{C}_{XY} \left(\mathcal{C}_{KK}+\lambda_j I\right)^{-1}$. Let $\mathcal{A}_0=\sum_{i,j=1}^{+\infty}a_{ij}\mu_i^{\frac{\beta}{2}}e_i\otimes \rho_j^{1-\frac{\gamma}{2}}f_j$ so that (ref) implies $\sum_{i,j=1}^{+\infty} a_{ij}^2 \leqslant B^2$. For the bias term, we show that $||\mathcal{A}_{\lambda}-\mathcal{A}_0||_{\beta',\gamma'}^2 \lesssim \sum_{j=1}^{y_N} j^{-\frac{\gamma'-\gamma}{q}} \lambda_j^{-(\beta-\beta')}\sum_{i=1}^{+\infty} a_{ij}^2 \lesssim B^2 \max _{1 \leqslant j \leqslant y_{N}} j^{-\frac{\gamma^{\prime}-\gamma}{q}} \lambda_{j}^{-\left(\beta-\beta^{\prime}\right)}$. The choice of $\lambda_j$ (ref) then implies that $||\mathcal{A}_{\lambda}-\mathcal{A}_0||_{\beta',\gamma'}^2 \lesssim N^{-\min\left\{\frac{\beta-\beta'}{\beta+p},\frac{\gamma'-\gamma}{1-\gamma}\right\}}$. For the variance term, we show that $||\hat{\mathcal{A}}-\mathcal{A}_{\lambda}||_{\beta',\gamma'}^2 \lesssim \sum_{j=1}^{y_{N}} j^{\frac{1-\gamma'}{q}} ||G_j||^2 \lambda_j^{-\beta'}|| U_j\rho_j^{\frac{1}{2}}f_j||^2$, where $G_j = \left(\mathcal{C}_{K K}+\lambda_{j} I\right)^{\frac{1}{2}}\left(\hat{\mathcal{C}}_{K K}+\lambda_{j} I\right)^{-1}\left(\mathcal{C}_{K K}+\lambda_{j} I\right)^{\frac{1}{2}} = \mathcal{O}(1)$ with high probability and $U_j = \left(\mathcal{C}_{KK}+\lambda_j I\right)^{-\frac{1}{2}}\left[ \hat{\mathcal{C}}_{XY} - \left(\hat{\mathcal{C}}_{KK}+\lambda_j I\right) \left( \mathcal{C}_{KK}+\lambda_j I\right)^{-1} \mathcal{C}_{XY} \right]$. We then apply Bernstein's inequality in Hilbert spaces to obtain a high-probability bound $\sum_{j=1}^{y_{N}} j^{\frac{1-\gamma'}{q}} \lambda_j^{-\beta'}|| U_j\rho_j^{\frac{1}{2}}f_j||^2 \lesssim \frac{1}{N} j^{\frac{1-\gamma}{q}}\lambda_{j}^{-(\beta'+\max\{\alpha-\beta,p\})} \lesssim N^{-\min\left\{\frac{\beta-\beta'}{\beta+p},\frac{\gamma'-\gamma}{1-\gamma}\right\}}$. \end{proof-sketch}

Regularization via Bias Contour

We have showed that if we learn all the spectral components under certain variance contour and regularize all other component can achieve optimal rate. In this section, we introduce another scheme to design the optimal estimator via learning all the spectral component under a certain bias contour. Specifically, we consider deciding the regularization strength according to the spectral elements induce a certain level of bias i.e. the bias contour $\ell_{C',\mathtt{bias}} = \left\{ (x,y)\in\mathbb{R}_+^2: x^{\frac{\beta-\beta'}{p}}y^{\frac{\gamma'-\gamma}{q}}= C'\right\}$. does not coincide with $\ell_{C,\mathtt{var}}$ for any $C'$ up to constant scaling. Thus, there exists a point $(x^*,y^*)$ on the variance contour with maximal contribution to bias. Naturally, we can also construct our estimator using a bias contour that passes through $(x^*,y^*)$. In this case, we may define $\lambda_j = \max\left\{ \left(j^{-\frac{\gamma'-\gamma}{q}}N^{\min\left\{ \frac{\beta-\beta'}{\max\{\alpha,\beta+p\}},\frac{\gamma'-\gamma}{1-\gamma} \right\}}\right)^{-\frac{1}{\beta-\beta'}}, c_0\left(\frac{N}{\log N}\right)^{-\frac{1}{\alpha}} \right\}$ for similar reasons as (ref), which also yields optimal rate as stated in (ref) below.

remark[On the optimal shape of regularization] The discussion in (ref) reveals another understanding of our information theoretic lower bound. Firstly, we should learn all the spectral components under the bias contour otherwise the bias will exceed the lower bound. Secondly, we should not learn any spectral component over the variance contour since otherwise the variance will exceed the lower bound. Thus \alertinline{the bias contour should always be under the variance contour}, otherwise no estimator can be designed. The bias and variance contours at the level of optimal learning rate are plotted in (ref). They only meet at $(x^*,y^*)$ with $x^*=1$ or $y^*=1$, which has the largest contribution to the bias (resp. variance) among all points on the variance (resp. bias) contour, thus dominating the estimation error. When the two curves meet at $y^*=1$, it reduces to the original kernel regression case. When the two curves meet at $x^*=1$, it leads to our new rate that depends on the output space.
theoremConsider the estimator $\hat{\mathcal{A}}$ defined by (ref) with $\lambda_j$ defined above. Suppose that (ref) hold, then there exists a universal constant $C$, such that \begin{equation} \notag \left\Vert\hat{\mathcal{A}}-\mathcal{A}_0\right\Vert_{\beta',\gamma'}^2 \leqslant C \tau^2 \left(\frac{N}{\log N}\right)^{-\min\left\{\frac{\beta-\beta'}{\max\{\alpha,\beta+p\}},\frac{\gamma'-\gamma}{1-\gamma}\right\}} \log^2 N \end{equation} holds with probability $\geqslant 1-e^{-\tau}$.

MultiLevel Kernel Operator Learning

In this section, we study a multilevel machine learning algorithm lye2021multi,li2021multilevel,boulle2022learning but at each level we consider a cost-accuracy trade-off de2022cost to control the variance at a proper scale. We show that the multilevel level algorithm can cover all the spectral component below the bias contour and achieve the optimal learning rate. Our idea is similar to the multilevel Monte Carlo giles2008multilevel,giles2015multilevel, which reduces bias from multilevel algorithm. Our multilevel estimator differs from the DeepONet lu2019deeponet and the PCA-Net bhattacharya2020model since we add different regularizations for each level. Our theory indicates that the multilevel approach outperforms previous ones and achieves the optimal learning rate.

figure*[figure* omitted — 810 chars of source]

The basic idea is to design a minimum number of machine learning estimators that cover all the spectral elements under the bias contour but do not exceed the variance contour at the same time. To achieve this, we choose sequences $\{x_i\}$ and $\{y_i\}$ for $1\leqslant i\leqslant L_N$ where $y_i$ denotes the $i$-th level and $x_i$ controls the corresponding regularization via the regularization coefficient $\lambda_i^{(K)}=x_i^{-\frac{1}{p}}$. The sequences are chosen in a staircase manner as plotted in (ref) (for formal definitions see (ref)). The eigenbasis $\left\{\rho_j^{\frac{1}{2}}f_j\right\}$ of the output space is divided into defferent levels by $\{y_i\}$. The main idea behind our multilevel method is that different levels of the output need to be learned with different regularization. Formally, we define our multilevel estimator as

equation[equation omitted — 283 chars of source]

The following theorem shows that the estimator (ref) can achieve the optimal convergence rate with $L_N=\mathcal{O}(\ln\ln N)$ when $\frac{\beta-\beta^{\prime}}{\max \{\alpha, \beta+p\}} \neq \frac{\gamma^{\prime}-\gamma}{1-\gamma}$. We also show that $O(\ln N)$ estimator is needed for the case when $\frac{\beta-\beta'}{\max\{\alpha,\beta+p\}}=\frac{\gamma'-\gamma}{1-\gamma}$ ((ref) Right) in (ref).

theoremSuppose that (ref) hold, then there exists a sequence $\{y_i\}_{1\leqslant i\leqslant L_N}$ with $L_N = \mathcal{O}(\ln N)$ when $\frac{\beta-\beta'}{\max\{\alpha,\beta+p\}}=\frac{\gamma'-\gamma}{1-\gamma}$ and $\mathcal{O}(\ln\ln N)$ otherwise, such that the estimator $\hat{\mathcal{A}}_{\mathtt{ml}}$ satisfies $\left\Vert\hat{\mathcal{A}}_{\mathtt{ml}}-\mathcal{A}_0\right\Vert_{\beta',\gamma'}^2 \leqslant C \tau^2 \left(\frac{N}{\log N}\right)^{-\min\left\{\frac{\beta-\beta'}{\max\{\alpha,\beta+p\}},\frac{\gamma'-\gamma}{1-\gamma}\right\}} \log^2 N$ with probability $\geqslant 1-e^{-\tau}$, where $C$ is a universal constant.
remark{ Our multilevel algorithm first apply the regression algorithm on low-frequency projections of the output samples with small regularization and then successively fine-tune the regression model on high-frequency projections of the output samples with stronger regularization, which matches the empirical use li2021multilevel,lye2021multi.}

Conclusion and Discussion

We considered the sample complexity of learning an operator between two infinite-dimensional Sobolev kernel Hilbert spaces. We provided an information theoretical lower bound for this problem along with a multi-level machine learning algorithm. Our lower bound is determined by the harder rate of two polynomial rates: one is fully determined by the hardness of input space while the other is fully controlled by the hardness of the output space. The second rate is new to the literature. We explained the our bound from the viewpoint of variance and bias counters in (ref) and (ref). The optimal estimator should learn all the spectral element under the bias contour but learn no information above the variance contour. To meet this requirement, we combined the idea of multilevel Monte Carlo with kernel operator learning, using successive levels to fit higher frequency information while keeping the variance at the same scale in order to reduce the bias. Our paper is the first paper on the non-parametric statistical optimality for multilevel algorithms. We leave estimation from discretely observe functional covariates with noise as future work zhoufunctional,zhou2022theory.

acknowledgementJikai Jin is partially supported by the elite undergraduate training program of School of Mathematical Sciences in Peking University. Yiping Lu is supported by the Stanford Interdisciplinary Graduate Fellowship (SIGF). Jose Blanchet is supported in part by the Air Force Office of Scientific Research under award number FA9550-20-1-0397. Lexing Ying is supported is supported by National Science Foundation under award DMS-2208163.

\printbibliography