Keywords: operator learning, Gaussian Sobolev operators, sample complexity, principal component analysis, least-squares approximation, Christoffel sampling

Corresponding author: maier@ins.uni-bonn.de (Gregor Maier)

1 Introduction↩︎

Many natural phenomena, despite their diverse individual nature, can be commonly described on an abstract level as learning an operator \[F : \mathcal{X}\to \mathcal{Y},\] which maps between infinite-dimensional spaces \(\mathcal{X}\) and \(\mathcal{Y}\). In computational science and engineering (CSE), for example, one is often interested in the solution operator of a continuum model which is governed by partial differential equations (PDEs) or variational inequalities. In this context, \(F\) typically maps inputs, such as parameters, boundary conditions, or coefficient functions, to outputs, such as states, fields, or observables. Prototypical examples with high and low operator regularity include holomorphic parameter-to-solution mappings in parametric PDEs and obstacle-to-solution mappings in obstacle problems, respectively. The latter exhibit only finite Sobolev regularity with respect to Gaussian measures. While there are many results on the approximation of smooth operators in the literature, comparatively few works are available for operators with only finite regularity which discuss both theoretical as well as practical approximation properties. To bridge this gap, we present in this work a practical, fully data-driven algorithm for learning Gaussian Sobolev operators from noisy, pointwise samples termed ‘Hermite-PCA approximation’. It combines empirical principal component analysis (PCA) for dimension reduction with weighted least-squares approximation based on Hermite polynomials. Our main contribution is a full error analysis of Hermite-PCA approximation – taking account of all sources of error in the problem – for the learning of operators with Gaussian Sobolev regularity. As we demonstrate, our algorithm is spectral – the smoother the operator, the faster the convergence – and, in the absence of discretization errors, its convergence rate is near-optimal. Moreover, we provide numerical experiments which empirically validate our theoretical results.

1.1 Motivations↩︎

A key question for practical applications is how to design suitable surrogate operators, i.e., approximations to \(F\), that achieve high accuracy with feasible computational costs. In Operator Learning, methods from machine learning, especially (deep) neural networks (NNs), are incorporated into the surrogate operator design in order to meet these demands. Many different designs, typically referred to as neural operators, along with theoretical error analyses have been proposed in recent years, often achieving impressive performances in numerical experiments. Among those are, for example, PCA-Net [1], [2], DeepONet [3], [4], FNO [5], [6], general neural operators [7], Poseidon [8], and variants thereof. We also refer to the reviews [9]–[11] and references therein.

However, neural operators face two drawbacks which are critical for their reliable deployment in CSE applications. First, they lack interpretability. Second, it is generally unclear whether a neural operator with desired performance guarantees can be obtained from practical optimization routines. For this reason, it is important to compare neural operators to simpler NN-free operator surrogates that alleviate these shortcomings. In [12], the authors present a framework for learning operators based on the theory of operator-valued reproducing kernel Hilbert spaces and Gaussian processes along with convergence guarantees and rigorous a priori error bounds. They provide numerical experiments which show that their kernel method is competitive to neural operators across various benchmark problems. In [13], the authors empirically compare the performance of neural and polynomial operator surrogates in extensive numerical experiments. Their key conclusion is that there is no universally best surrogate model. Polynomial surrogates seemingly perform significantly better than neural operators for smooth inputs, whereas neural operators tend to be superior for rough input data. However, the former are typically much cheaper computationally to train, as they rely on simple procedures such as linear least-squares.

These findings motivate to put more focus on the design of NN-free operator surrogates. This work contributes to this agenda, in that we provide a polynomial surrogate model based on generalized Wiener-Hermite polynomial chaos expansions together with a non-intrusive training procedure via weighted-least squares approximation with carefully chosen pointwise training samples.

1.2 Contributions↩︎

We consider a Lipschitz continuous operator \(F : \mathcal{X}\to \mathcal{Y}\) between real separable Hilbert spaces \(\mathcal{X}, \mathcal{Y}\). We equip \(\mathcal{X}\) with an unknown Gaussian measure \(\mu\) and denote the (unknown) PCA eigenvalues of its covariance operator by \(\lambda_1 \geq \lambda_ 2 \geq \cdots > 0\). We assume pointwise access to \(F\), with operator evaluations corrupted by additive noise with noise level \(\sigma \geq 0\). Our surrogate operator design has the standard encoder-decoder structure, \(\widehat{F}= \widehat{\mathcal{D}}_{\mathcal{Y}}\circ \widehat{f}\circ \widehat{\mathcal{E}}_{\mathcal{X}}\), consisting of an encoder \(\widehat{\mathcal{E}}_{\mathcal{X}}: \mathcal{X}\to \mathbb{R}^{d_{\mathcal{X}}}\), a decoder \(\widehat{\mathcal{D}}_{\mathcal{Y}}: \mathbb{R}^{d_{\mathcal{Y}}} \to \mathcal{Y}\), and a latent function \(\widehat{f}: \mathbb{R}^{d_{\mathcal{X}}} \to \mathbb{R}^{d_{\mathcal{Y}}}\), see Figure 1.

We present an algorithm to construct \(\widehat{F}\), which has two main parts. First, we learn the encoder and decoder by empirical PCA based on \(N_{\mathcal{X}}\) and \(N_{\mathcal{Y}}\) samples, respectively. Second, we construct a suitable \(s\)-dimensional approximation space of vector-valued linear combinations of Hermite polynomials, from which to select \(\widehat{f}\). The Hermite coefficients of \(\widehat{f}\) are then learned by weighted least-squares approximation based on \(M\) labeled training samples. For this, we employ Christoffel sampling methods to identify an optimal distribution of the input training data.

Note that the encoder induces a Gaussian measure \(\hat{\varrho} : = \widehat{\mathcal{E}}_{\mathcal{X}} \sharp \mu\) on \(\mathbb{R}^{d_{\mathcal{X}}}\). We can now state an informal version of our main result.

In fact, the previous result follows from a more general error bound which holds for arbitrary Lipschitz operators and does not require any additional Sobolev regularity, see Theorem 2.

We now summarize our main contributions.

• Practicality: We present a fully data-driven algorithm to construct a polynomial surrogate operator for Lipschitz operators with Sobolev regularity from pointwise, noisy samples. The only assumption in the problem setup is that \(\mu\) is a Gaussian measure, but its covariance structure is unknown and needs to be learned from data. Our algorithm employs linear least-squares fitting, and is therefore very computationally efficient.

• Weak assumptions: Our measure \(\mu\) has unbounded support. This complicates the analysis but is the most adequate support condition under minimal assumptions on the input data distribution. This situation has, to best of our knowledge, not been adequately addressed in the literature so far, as existing results typically require bounded support assumptions. In addition, we also allow for arbitrary decay of the PCA eigenvalues \(\lambda_i\), not just algebraic as typically considered in the literature.

• Full error analysis: The bound ?? takes into account all sources of errors originating from the model design and learning procedure. It bounds the overall approximation error by (empirical) PCA projection errors, an approximation error in latent space, and a noise error. It explicitly shows the effect of the (hyper-)parameters of the problem to each error term, thus enabling optimal design choices for practical applications.

• Optimal approximation rates: In the absence of the PCA error terms and noise, we show that the error bound ?? achieves optimal approximation rates. In fact, we argue that the approximation error admits a matching lower bound. We provide explicit examples how \(w_s\) decays in \(s\) in the case of an infinite-dimensional input space, \(\dim(\mathcal{X}) = \infty\), and algebraically and exponentially decaying \(\lambda_i\).

• Spectral approximation property: The bound ?? reveals that the approximation error is determined by the weight \(w_s^k\) in the case where the \(\widehat{\mathcal{E}}_{\mathcal{Y}}\circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}}\) is \(k\)-Sobolev regular. As the algorithm itself is independent of \(k\), it thus automatically achieves faster convergence rates for higher regular Sobolev operators. In particular, for infinitely-smooth operators the error decays faster than \(w_s^k\) as \(s \rightarrow \infty\) for any \(k\).

• Universality: Our algorithm chooses a certain polynomial space for Sobolev regular operators. However, it can be easily adapted to other regularity, e.g., mixed regularity, as well. In fact, our main theoretical contributions in their most general form consider arbitrary polynomial spaces and impose no regularity assumption on the operator.

• Numerical experiments: We provide empirical validation of our theoretical results by applying our algorithm to learn (i) the obstacle-to-solution operator of an obstacle problem and (ii) functionals of varying Sobolev regularity.

1.3 Related work↩︎

PCA methodologies for dimensionality reduction are prominent among classical reduced basis methods [14]. In the context of operator learning, they are employed in [15], [16], though only for the decoder in the output space. The work [1] introduces PCA-Net, which applies empirical PCA for both the encoder and decoder in the input and output space, respectively, and is thus closest to the encoder-decoder scheme used in this work. Its theoretical properties are further analyzed in [2]. PCA-Net models the latent function \(\widehat{f}\) as a deep NN. In contrast, we use Hermite polynomials to obtain a fully NN-free approximation scheme. We mention in passing that there are other neural operator designs, which are not based on NN approximations, but instead, e.g., on the random feature model [17]–[20], or polynomial chaos expansions [21].

A long list of recent works established that holomorphic (nonlinear) operators can be approximated efficiently, that is, with algebraic or even exponential convergence rates. This is both in terms of NN expression rate bounds [22]–[24] as well as sample complexity estimates [25]–[30]. Such operators arise, e.g., as data-to-solution operators in parametric PDEs and in problems of uncertainty quantification [31]. Algebraic convergence results are also available for approximating operators of special structure, such as linear operators [32]–[34], and pseudo-differential operators [35]. Complexity estimates for non-smooth operators, as we consider in this work, are given in [36]–[39]. They reveal intrinsic curses of parametric or sample complexities for non-smooth operators on infinite-dimensional domains: no approximation algorithm can achieve an algebraically decaying worst-case error in terms of NN parameters or employed data samples. The experiments in the present work are, to the best of our knowledge, the first numerical evidence of this phenomenon, and ours is the first fully-data driven algorithm that achieves such rates, while also tackling practical issues such as noise and unknown \(\mu\). We also mention in passing the recent work [40], which establishes upper (subalgebraic) expression rate bounds with respect to a subgaussian input measure for the approximation of Fréchet differentiable operators with PCA-Net. A concise overview of current results in operator learning theory is given in [41].

The recovery of operators from finite noisy samples is studied in [42]–[44]. While the first two works only derive upper bounds with a focus on Lipschitz and holomorphic operators, respectively, the authors of [44] prove information-theoretic upper and matching (or near-matching) lower bounds for the minimax risk of Lipschitz and Hölder operators. Their findings again confirm in a minimax sense a curse of sample complexity for finitely regular operators. However, in contrast to our results, all these works require uniform boundedness assumptions for the objective operator. Moreover, they do not establish practical approximation algorithms nor provide any numerical results, as we do in this work.

Weighted least-squares approximation and Christoffel sampling techniques, originally developed for scalar-valued approximations to achieve near-optimal sample complexities [45], have recently been generalized to the Hilbert-valued setting [30], [31], [46], [47]. An operator extension of Christoffel sampling is considered in [48]. But this work does not address the practical scenario where the measure \(\mu\) is unknown, as we do in this work. It also does not consider algorithms that achieve concrete (and, as we establish, near-optimal) approximation rates for Sobolev operators, which is a major focus of this paper.

The Gaussian setting in this paper is also adopted in [22] and [38], [39] to derive expression rate and sample complexity estimates for operator learning, respectively, and in [29] to study problems in uncertainty quantification with Gaussian random field inputs. However, these works are just theoretical studies and do not provide numerical evidence. In [49] the authors adopt the Gaussian setting to analyze approximation errors in derivative-informed operator learning strategies and also perform numerical experiments. However, in all the works listed above, the authors assume the Gaussian measure \(\mu\) to be known. In contrast, we assume \(\mu\) to be unknown and approximate it from data – a major asset of our Hermite-PCA method making it fully data-driven.

1.4 Outline↩︎

The rest of the paper is organized as follows. Section 2 introduces our setup and main problem in detail. In Section 3, we describe the surrogate operator construction and discuss the resulting Hermite-PCA algorithm. Section 4 contains our main result, Theorem 1, which is a detailed error analysis of the Hermite-PCA approximation for Gaussian Sobolev operators and which we empirically validate in numerical experiments in Section 5. In Section 6, we present a general error bound for the Hermite-PCA method, see Theorem 2. Together with an \(\ell^2\)-characterization of Gaussian Sobolev spaces in Section 7 we use this result to prove Theorem 1 in Section 8. We finally conclude in Section 9 with limitations and future work. The appendix provides a number of further definitions and results which are used in the main part of the paper.

2 Setup and main problem↩︎

We commence by presenting our setup in detail and state the main problem.

2.1 Setup↩︎

Let \(F : \mathcal{X}\to \mathcal{Y}\) be an \(L\)-Lipschitz continuous operator between two separable Hilbert spaces for some \(L > 0\), that is, \[{\left\|F(X) - F(X')\right \|}_{\mathcal{Y}} \leq L {\left\|X - X'\right \|}_{\mathcal{X}}, \quad \forall X, X' \in \mathcal{X}.\] We assume that we can access \(F\) only through noisy, pointwise evaluation, \[Y = F(X) + \sigma E,\] where the scalar \(\sigma \geq 0\) denotes the noise level, and \(E\) is a \(\mathcal{Y}\)-valued random noise variable.

Recall that a random vector \(Z\) in a separable Hilbert space \(\mathcal{H}\) is subgaussian with parameter \(K > 0\) if \({\left\|Z\right \|}_{\mathcal{H}}\) is subgaussian in the conventional sense, that is, \[\label{eq:32subgaussian32moment32condition} \mathbb{E}[ {\left\|Z\right \|}_{\mathcal{H}}^p ]^{1/p} \leq K \sqrt{p}, \quad \forall p \geq 1.\tag{1}\] We call the distribution of a subgaussian random variable with parameter \(K\) a subgaussian distribution with parameter \(K\). As an example, the measure \(\mu\) on \(\mathcal{X}\) is subgaussian. Indeed, by [50], we have \[\label{eq:32tail32bound32for32mu} \mathbb{P}_{X \sim \mu} [{\left\|X\right \|}_{\mathcal{X}} > t] = \mu\left( \left\{X \in \mathcal{X}: {\left\|X\right \|}_{\mathcal{X}} > t \right\} \right) \leq 4 \exp\left(- \frac{t^2}{8 \mathrm{tr}(\Sigma)} \right), \quad \forall t \geq 0.\tag{2}\] A standard computation (see, e.g., the proof of [51]) shows that this is equivalent to 1 for some parameter \(K_{\mu} > 0\), which differs from \(\mathrm{tr}(\Sigma)\) only by an absolute constant.

We write \(\nu = \nu(\sigma)\) for the distribution of \(Y = F(X) + \sigma E\), where \(X \sim \mu\) and \(E \sim \rho\) are mutually independent. It follows directly from the moment condition 1 that \(\nu\) is subgaussian with parameter \(K_{\nu} = {\left\|F(0)\right \|}_{\mathcal{Y}} + L K_{\mu} + \sigma\). Notice that in the noiseless case, we have \(\nu(\sigma = 0) = F \sharp \mu\).

2.2 Main problem↩︎

We now consider the approximation of \(F\). We follow the standard operator learning approach and construct the surrogate operator \(\widehat{F}\) from three parts:

1. an encoder for the input space, \(\widehat{\mathcal{E}}_{\mathcal{X}} : \mathcal{X}\to \mathbb{R}^{d_{\mathcal{X}}}\),

2. a decoder for the output space, \(\widehat{\mathcal{D}}_{\mathcal{Y}} : \mathcal{Y}\to \mathbb{R}^{d_{\mathcal{Y}}}\),

3. and a latent space function \(\widehat{f} : \mathbb{R}^{d_{\mathcal{X}}} \to \mathbb{R}^{d_{\mathcal{Y}}}\).

Here, \(d_{\mathcal{X}}, d_{\mathcal{Y}} \in \mathbb{N}\), with \(d_{\mathcal{X}} \leq \dim(\mathcal{X})\), \(d_{\mathcal{Y}} \leq \dim(\mathcal{Y})\), denote the encoding and decoding dimension, respectively. We then define the approximation \(\widehat{F} \approx F\) as \[\widehat{F} := \widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{f} \circ \widehat{\mathcal{E}}_{\mathcal{X}},\] see Figure 1. In order to compute this approximation, we require three datasets.

1. To compute \(\widehat{f}\), we generate \(M\) labeled data samples \((X_i, Y_i) \in \mathcal{X}\times \mathcal{Y}\) with \[\label{eq:32X95i} X_i \sim_{\mathrm{i.i.d.}} \mu_{\mathrm{samp}}, \qquad i = 1, \dots, M\tag{3}\] where \(\mu_{\mathrm{samp}}\) is a custom sampling measure, which we will specify later on, and \[\label{eq:32Y95i} Y_i = F(X_i) + \sigma E_i, \quad E_i \sim_{\mathrm{i.i.d.}} \rho, \qquad i = 1, \dots, M,\tag{4}\] are the corresponding noisy operator evaluations. We assume that the \(E_i\) and the \(X_j\) are mutually independent so that the \(Y_i\) are independent and identically distributed.

2. To compute \(\widehat{\mathcal{E}}_{\mathcal{X}}\) we require an additional amount of \(N_{\mathcal{X}} \geq d_{\mathcal{X}}\) unlabeled data points \[\widehat{X}_i \sim_{\mathrm{i.i.d.}} \mu, \qquad i = 1, \dots, N_{\mathcal{X}}.\]

3. To compute \(\widehat{\mathcal{D}}_{\mathcal{Y}}\) we require \(d_{\mathcal{Y}} \leq N_{\mathcal{Y}} \leq N_{\mathcal{X}}\) noisy labels \[\widehat{Y}_i = F(\widehat{X}_i) + \sigma \widehat{E}_i, \quad \widehat{E}_i \sim_{\mathrm{i.i.d.}} \rho, \qquad i = 1, \dots, N_{\mathcal{Y}}.\] Similarly as above, we assume that the \(\widehat{E}_i\) and the \(\widehat{X}_j\) are mutually independent so that the \(\widehat{Y}_i\) are i.i.d. with respect to \(\nu\).

Throughout, we measure the accuracy of \(\widehat{F}\) using the \(L^2_{\mu}(\mathcal{X}; \mathcal{Y})\)-norm. We write \(L_{\mu}^2(\mathcal{X}; \mathcal{Y})\) for the Lebesgue-Bochner space of (equivalence classes of) strongly measurable operators \(F: \mathcal{X}\to \mathcal{Y}\) with finite Bochner norm \[{\left\|F\right \|}_{L_{\mu}^2(\mathcal{X}; \mathcal{Y})} := \left( \int_{\mathcal{X}} {\left\|F(X)\right \|}_{\mathcal{Y}}^2 \,\mathrm{d}\mu(X) \right)^{1/2}.\]

3 Construction of the approximation↩︎

We now describe the construction of \(\widehat{F}\), which we divide into the three components (i), (ii), and (iii) from Section 2.2.

3.1 Encoder and decoder construction: PCA and empirical PCA↩︎

We compute \(\widehat{\mathcal{E}}_{\mathcal{X}}\) and \(\widehat{\mathcal{D}}_{\mathcal{Y}}\) using empirical PCA, which we now describe. Recall that \(\Sigma = \Sigma_{\mu}\) is the covariance operator of \(\mu\). We write \(\lambda_1^{\mu} \geq \lambda_2^{\mu} \geq \cdots > 0\) for the corresponding PCA eigenvalues and \(\{\phi^{\mu}_i\}_{i = 1}^{\dim(\mathcal{X})}\) for the PCA eigenvectors, which form an orthonormal basis of \(\mathcal{X}\). To simplify notation, we write \(\lambda_i^{\mu} = \lambda_i\) and \(\phi^{\mu}_i = \phi_i\) if not specified otherwise. Given an encoding dimension \(d_{\mathcal{X}}\), we define the true PCA encoder and decoder for \(\mathcal{X}\) as, respectively, \[\label{eq:32true32PCA32encoder32decoder32for32X} \mathcal{E}_{\mathcal{X}}(X) := (\left\langle X, \phi_1 \right\rangle_{\mathcal{X}}, \dots, \left\langle X, \phi_{d_{\mathcal{X}}} \right\rangle_{\mathcal{X}}), \quad \forall X \in \mathcal{X}, \quad \mathrm{ and } \quad \mathcal{D}_{\mathcal{X}}(\boldsymbol{x}) := \sum_{i = 1}^{d_{\mathcal{X}}} x_i \phi_i, \quad \forall \boldsymbol{x} \in \mathbb{R}^{d_{\mathcal{X}}}.\tag{5}\] The construction for the output space \(\mathcal{Y}\) is based on the distribution \(\nu\) of the random variable \(Y = F(X) + E\) with \(X \sim \mu\), \(E \sim \rho\). Its covariance operator is given by \[\Sigma_{\nu} := \mathbb{E}[(Y - \mathbb{E}[Y]) \otimes (Y - \mathbb{E}[Y])],\] with eigenvalues \(\lambda_1^{\nu} \geq \lambda_2^{\nu} \geq \dots \geq 0\) and corresponding eigenvectors \(\{\psi_i\}_{i = 1}^{\dim(\mathcal{Y})}\), which form an orthonormal basis of \(\mathcal{Y}\), see Remark 1. The associated true encoder and decoder mappings are given by, respectively, \[\mathcal{E}_{\mathcal{Y}}(Y) := (\left\langle Y, \psi_1 \right\rangle_{\mathcal{Y}}, \dots, \left\langle Y, \psi_{d_{\mathcal{Y}}} \right\rangle), \quad \forall Y \in \mathcal{Y}, \quad \mathrm{ and } \quad \mathcal{D}_{\mathcal{Y}}(\boldsymbol{y}) := \sum_{i = 1}^{d_{\mathcal{Y}}} y_i \psi_i, \quad \forall \boldsymbol{y} \in \mathbb{R}^{d_{\mathcal{Y}}}.\]

Note, however, that \(\mathcal{E}_{\mathcal{X}}\) and \(\mathcal{D}_{\mathcal{Y}}\) cannot be used for practical computations, as the measure \(\mu\) is assumed to be unknown and hence, so are the basis functions \(\phi_i\) and \(\psi_i\). Instead, we use empirical PCA. Given \(N\) unlabeled data points \(\widehat{X}_i\) in (b) above, we first define the empirical measure and the associated empirical covariance operator by, respectively, \[\hat{\mu} := \frac{1}{N_{\mathcal{X}}} \sum_{i = 1}^{N_{\mathcal{X}}} \delta_{\widehat{X}_i} \quad \mathrm{ and } \quad \Sigma_{\hat{\mu}} := \frac{1}{N_{\mathcal{X}}} \sum_{i = 1}^{N_{\mathcal{X}}} \widehat{X}_i \otimes \widehat{X}_i.\] The latter operator has eigenvalues \(\lambda_1^{\hat{\mu}} \geq \lambda_2^{\hat{\mu}} \geq \dots \geq 0\) with \(\lambda_i^{\hat{\mu}} = 0\) for all \(N_{\mathcal{X}} < i \leq \dim(\mathcal{X})\) and corresponding eigenvectors \(\{\widehat{\phi}_i\}_{i = 1}^{\dim(\mathcal{X})}\), which form an orthonormal basis of \(\mathcal{X}\). If not specified otherwise, we write \(\Sigma_{\hat{\mu}} = \widehat{\Sigma}\) and \(\lambda_i^{\hat{\mu}} = \hat{\lambda}_i\) in what follows. Given this, we define the empirical PCA encoder and decoder for \(\mathcal{X}\) as, respectively, \[\label{eq:32empirical32PCA32encoder32decoder32for32X} \widehat{\mathcal{E}}_{\mathcal{X}}(X) := (\langle X, \widehat{\phi}_1 \rangle_{\mathcal{X}}, \dots, \langle X, \widehat{\phi}_{d_{\mathcal{X}}} \rangle_{\mathcal{X}}), \quad \forall X \in \mathcal{X}, \quad \mathrm{ and } \quad \widehat{\mathcal{D}}_{\mathcal{X}}(\boldsymbol{x}) := \sum_{i = 1}^{d_{\mathcal{X}}} x_i \widehat{\phi}_i, \quad \forall \boldsymbol{x} \in \mathbb{R}^{d_{\mathcal{X}}}.\tag{6}\] The construction for \(\mathcal{Y}\) is analogous. Given \(N_{\mathcal{Y}}\) data points \(\widehat{Y}_i\) as in (c) above, we define the empirical covariance operator \[\Sigma_{\hat{\nu}} := \frac{1}{N_{\mathcal{Y}}} \sum_{i = 1}^{N_{\mathcal{Y}}} (\widehat{Y}_i - \mathbb{E}[\widehat{Y}_i]) \otimes (\widehat{Y}_i - \mathbb{E}[\widehat{Y}_i]),\] with eigenvalues \(\lambda_1^{\hat{\nu}} \geq \lambda_2^{\hat{\nu}} \geq \dots \geq 0\) with \(\lambda_i^{\hat{\nu}} = 0\) for all \(N_{\mathcal{Y}} < i \leq \dim(\mathcal{Y})\) and corresponding eigenvectors \(\{\widehat{\psi}_i\}_{i = 1}^{\dim(\mathcal{Y})}\), which form an orthonormal basis of \(\mathcal{Y}\). The empirical PCA encoder and decoder for \(\mathcal{Y}\) are defined as, respectively, \[\widehat{\mathcal{E}}_{\mathcal{Y}}(Y) := (\langle Y, \widehat{\psi}_1 \rangle_{\mathcal{Y}}, \dots, \langle Y, \widehat{\psi}_{d_{\mathcal{Y}}} \rangle_{\mathcal{Y}}), \quad \forall Y \in \mathcal{Y}, \quad \mathrm{ and } \quad \widehat{\mathcal{D}}_{\mathcal{Y}}(\boldsymbol{y}) := \sum_{i = 1}^{d_{\mathcal{Y}}} y_i \widehat{\psi}_i, \quad \forall \boldsymbol{y} \in \mathbb{R}^{d_{\mathcal{Y}}}.\]

3.2 Latent space approximation construction: Least-squares approximation with Hermite polynomials↩︎

Our latent space approximation \(\widehat{f}\) is constructed using certain orthogonal polynomials. We now describe this construction, along with the construction of the sampling measure \(\mu_{\mathrm{samp}}\). For \(n \in \mathbb{N}_0\), we define the \(n\)th normalized (probabilist’s) Hermite polynomial on \(\mathbb{R}\) by \[H_n: \mathbb{R}\to \mathbb{R},H_n(x) := \frac{(-1)^n}{\sqrt{n!}} \exp \left( \frac{x^2}{2} \right) \frac{d^n}{d x^n} \exp \left( -\frac{x^2}{2} \right).\] The higher-dimensional Hermite polynomials are defined as products of the one-dimensional ones. Given \(d \in \mathbb{N}\), a sequence \(\boldsymbol{\lambda} = (\lambda_i)^{d}_{i=1}\) with \(\lambda_i > 0\), and a multi-index \(\boldsymbol{\gamma} = (\gamma_i)^{d}_{i=1} \in \mathbb{N}^d_0\), we define \[H_{\boldsymbol{\gamma},\boldsymbol{\lambda}} : \mathbb{R}^d \to \mathbb{R},\quad H_{\boldsymbol{\gamma},\boldsymbol{\lambda}} (\boldsymbol{x}) := \prod_{i = 1}^{d} H_{\gamma_i} \left( \frac{x_i}{\sqrt{\lambda_i}} \right).\] Notice that the family \(\{ H_{\boldsymbol{\gamma},\boldsymbol{\lambda} } \}_{\boldsymbol{\gamma} \in \mathbb{N}^d_0 }\) forms an orthonormal basis of \(L^2_{\varrho}(\mathbb{R}^d)\), where \(\varrho = \varrho_{\boldsymbol{\lambda}} : = \mathcal{N}(0,\boldsymbol{\lambda})\) is the Gaussian measure with mean zero and diagonal covariance with \(i\)th diagonal entry \(\lambda_i\).

Let \(S \subset \mathbb{N}^{d_{\mathcal{X}}}_0\) be a set of size \(|S| = s\). We will make a specific choice of \(S\) later. Then we consider a latent space approximation \(\widehat{f}\) as \[\widehat{f} = \sum_{\boldsymbol{\gamma} \in S} \boldsymbol{c}_{\boldsymbol{\gamma}} H_{\boldsymbol{\gamma},\boldsymbol{\hat{\lambda}}},\] where \(\boldsymbol{\hat{\lambda}} = (\hat{\lambda}_i)^{d_{\mathcal{X}}}_{i=1}\) is the vector of empirical PCA eigenvalues and \(\boldsymbol{c}_{\boldsymbol{\gamma}} \in \mathbb{R}^{d_{\mathcal{Y}}}\) are vector-valued coefficients. We compute these via (weighted) least squares. With the above notation, let \(\hat{\varrho} = \hat{\varrho}_{\boldsymbol{\hat{\lambda}}} := \mathcal{N}(0,\boldsymbol{\hat{\lambda}})\) and define the subspace \[\label{eq:32hat-f32subspace32def} \mathcal{P}= \mathcal{P}_{\mathbb{R}^{d_{\mathcal{Y}}}} := \left \{ \sum_{\boldsymbol{\gamma} \in S} \boldsymbol{c}_{\boldsymbol{\gamma}} H_{\boldsymbol{\gamma},\boldsymbol{\hat{\lambda}}} : \boldsymbol{c}_{\boldsymbol{\gamma}} \in \mathbb{R}^{d_{\mathcal{Y}}} \right \} \subseteq L^2_{\hat{\varrho} } (\mathbb{R}^{d_{\mathcal{X}} } ; \mathbb{R}^{d_{\mathcal{Y}}} ).\tag{7}\] Now define the sampling measure on \(\mathbb{R}^{d_{\mathcal{X}}}\) \[\label{eq:32upsilon-samp32def} \,\mathrm{d}\varrho_{\mathrm{samp}} : = w^{-1} \,\mathrm{d}\hat{\varrho}, ~~\text{where}~~ w(\boldsymbol{x}) := \left ( \frac{1}{s} \sum_{\boldsymbol{\gamma} \in S} H_{\boldsymbol{\gamma},\boldsymbol{\hat{\lambda}}}(\boldsymbol{x})^2 \right )^{-1},\quad \boldsymbol{x} \in \mathbb{R}^{d_{\mathcal{X}}},\tag{8}\] and the corresponding sampling measure on \(\mathcal{X}\) \[\mu_{\mathrm{samp}} := \widehat{\mathcal{D}}_{\mathcal{X}} \sharp \varrho_{\mathrm{samp}}.\] Next, consider \(M\) labeled data samples \((X_i,Y_i)\) as in (a), and given by 3 4 . Then we define \(\widehat{f}\) as a minimizer of the weighted least-squares fit \[\label{eq:32f-hat32def} \widehat{f} \in {\underset{p \in \mathcal{P}}{\operatorname{argmin}}} \frac{1}{M} \sum^{M}_{i=1} w(\widehat{\mathcal{E}}_{\mathcal{X}}(X_i)) {\|\widehat{\mathcal{E}}_{\mathcal{Y}}(Y_i) - p(\widehat{\mathcal{E}}_{\mathcal{X}}(X_i))\|}^2_2.\tag{9}\]

In practice, 9 can be solved by solving a set of \(d_{\mathcal{Y}}\) algebraic least-squares problems of size \(M \times s\). Indeed, it is a short exercise to show that \[\label{eq:32L-S32approximation32rewritten} \widehat{f}(\boldsymbol{x}) = \left( \sum_{\boldsymbol{\gamma} \in S} c^k_{\boldsymbol{\gamma}} H_{\boldsymbol{\gamma}, \boldsymbol{\hat{\lambda}} }(\boldsymbol{x}) \right )^{d_{\mathcal{Y}}}_{k=1} , \qquad \boldsymbol{c^k} : = (c^k_{\boldsymbol{\gamma}} )_{\boldsymbol{\gamma} \in S} \in {\underset{\boldsymbol{c} \in \mathbb{R}^s}{\operatorname{argmin}}}{{\|\boldsymbol{A} \boldsymbol{c} - \boldsymbol{b^k}\|}_2^2}, \quad k = 1, \dots, d_{\mathcal{Y}},\tag{10}\] where the weighted measurement matrix and vectors are given by \[\label{eq:32A4432b94k} \boldsymbol{A} = \left( \sqrt{\frac{w(\boldsymbol{x_i})}{M}} H_{\boldsymbol{\gamma}, \boldsymbol{\hat{\lambda}} }(\boldsymbol{x_i}) \right)_{i \in [M],\boldsymbol{\gamma} \in S} \in \mathbb{R}^{M \times s}, \qquad \boldsymbol{b^k} = \left( \sqrt{\frac{w(\boldsymbol{x_i})}{M}} \left\langle Y_i, \widehat{\psi}_i \right\rangle_{\mathcal{Y}} \right)_{i \in [M]} \in \mathbb{R}^M,\tag{11}\] and \(\boldsymbol{x}_i = \widehat{\mathcal{E}}_{\mathcal{X}}(X_i)\) for \(i \in [M]\). Here and elsewhere, we use the notation \([d] = \{1, \dots, d\}\) and \([\infty] = \mathbb{N}\).

It remains to specify the multi-index set \(S\). This should be chosen to give as good an approximation as possible from the resulting subspace \(\mathcal{P}\), which naturally depends on the regularity assumptions placed on \(F\). In this work, we assume \(F\) has Sobolev regularity. As we shall show, in this case, a suitable choice of \(S\) is defined as follows. Suppose that \(\hat{\lambda}_{d_{\mathcal{X}}} > 0\), introduce the sequence \[\label{eq:32hat40v4195gamma} \hat{v}_{\boldsymbol{\gamma}} := \left( 1 + \sum_{i = 1}^{d_{\mathcal{X}}} \frac{\gamma_i}{\hat{\lambda}_{i}}\right)^{-1/2}, \quad \boldsymbol{\gamma} \in \mathbb{N}_0^{d_{\mathcal{X}}},\tag{12}\] and let \[\label{eq:32hat40tau41} \hat{\tau} : \mathbb{N}\to \mathbb{N}_0^{d_\mathcal{X}}\tag{13}\] be a bijection that gives a nonincreasing rearrangement, i.e., \(\hat{v}_{\boldsymbol{\hat{\tau}(1)}} \geq \hat{v}_{\boldsymbol{\hat{\tau}(2)}} \geq \cdots > 0\). Then we set \[\label{eq:32S32hat-tau32def} S = \{ \boldsymbol{\hat{\tau}(1)},\ldots,\boldsymbol{\hat{\tau}(s)} \}.\tag{14}\] This specific choice for \(S\) comes from the analysis of the approximation error. Later, in Section 4.3 we explain that it leads to an approximation error that decays with optimal rates in \(s\) under the stipulated regularity assumptions.

3.3 Algorithm and practical aspects↩︎

With all the pieces in place, we summarize the computation of \(\widehat{F}\) in Algorithm 2. A key facet of this algorithm is that the main steps can also be performed numerically. We now discuss these steps in more detail.

Step 1 assumes we are given \(N\) unlabeled samples from the unknown measure \(\mu\). This is a reasonable assumption in practice. In Step 2 we generate \(N_{\mathcal{Y}} \leq N_{\mathcal{X}}\) noisy evaluations of the target operator \(F\). In Step 3, we use the results of Steps 1 and 2 to compute the empirical encoders and decoders. This is done by first forming the empirical covariance operators \(\Sigma_{\hat{\mu}}\) and \(\Sigma_{\hat{\nu}}\) and second computing, respectively, their first \(d_{\mathcal{X}}\) and \(d_{\mathcal{Y}}\) eigenvalues and eigenvectors. These are operators on the Hilbert spaces \(\mathcal{X}\) and \(\mathcal{Y}\), respectively. In practice, these spaces are normally discretized first using, for example, finite elements. In which case this computation involves computing the first \(d_{\mathcal{X}}\) and \(d_{\mathcal{Y}}\) eigenvalues and eigenvectors of matrices of size \(D_{\mathcal{X}} \times D_{\mathcal{X}}\) and \(D_{\mathcal{Y}} \times D_{\mathcal{Y}}\), respectively, where \(D_{\mathcal{X}},D_{\mathcal{Y}}\) are the sizes of the discretizations.

Step 4 requires the computation of the nonincreasing rearrangement \(\hat{\tau}\) in 13 . This can be done efficiently by using a Dijkstra-like algorithm on the weighted lattice \(\mathbb{N}_0^{d_{\mathcal{X}}}\), where all edges in direction \(i \in [d_{\mathcal{X}}]\) carry the weight \(1 / \lambda_i\), to iteratively select those nodes which can be reached from the origin with the smallest costs.

Step 5 involves sampling \(M\) points i.i.d.from \(\varrho_{\mathrm{samp}}\). This measure can be expressed as an additive mixture (with weights \(1/s\)) of the probability measures \[H_{\boldsymbol{\gamma},\boldsymbol{\hat{\lambda}}}(\boldsymbol{x})^2 \,\mathrm{d}\hat{\varrho}(\boldsymbol{x}) = \bigotimes^{d_{\mathcal{X}}}_{i=1} H_{\gamma_i} \Bigg ( \frac{x_i}{\sqrt{\hat{\lambda}_i} } \Bigg )^2 \frac{1}{\sqrt{2 \pi \hat{\lambda}_i}} e^{-x^2_i / (2 \hat{\lambda}_i)} \,\mathrm{d}x_i,\] which are tensor products of one-dimensional probability measures. Thus, to sample \(\boldsymbol{x} \sim \varrho_{\mathrm{samp}}\), one proceeds as follows. First, draw an index \(\boldsymbol{\gamma} \in S\) uniformly at random from all possible \(|S| = s\) indices. Then, draw \(\boldsymbol{x} \sim H_{\boldsymbol{\gamma},\boldsymbol{\hat{\lambda}}}(\boldsymbol{x})^2 \,\mathrm{d}\hat{\varrho}(\boldsymbol{x})\) by drawing its components \(x_i\) independently from the corresponding univariate measures. Note that this latter step can be efficiently carried out via, for example, inverse transform sampling [52].

Having done this, Step 6 just involves sampling \(F\) at the sample points \(X_1,\ldots, X_M\). Finally, in Step 7 we compute \(\widehat{f}\) by solving the algebraic least-squares problems 10 11 . These can be solved efficiently via, e.g., conjugate gradients (notice that our theoretical results also guarantee that \(\boldsymbol{A}\) is well-conditioned) in \(\mathcal{O}(s M d_{\mathcal{Y}})\) flops.

4 Hermite-PCA approximation of Sobolev operators↩︎

We now state and discuss our main result on the approximation of Sobolev operators via the Hermite-PCA method, Algorithm 2.

4.1 Main result↩︎

To state this result, we need several further concepts. Let \(\Gamma := \{ \boldsymbol{\gamma} \in \mathbb{N}^{\mathbb{N}}_0 : |\mathrm{supp}(\boldsymbol{\gamma}) |< \infty \},\) where \(\mathrm{supp}(\boldsymbol{\gamma}) = \{ i \in \mathbb{N}: \gamma_i \neq 0 \}\), be the set of infinite multi-indices with only finitely-many nonzero entries. Define the weights \[\label{eq:32v95gamma} v_{\boldsymbol{\gamma}} := \left( 1 + \sum_{i = 1}^{\dim(\mathcal{X})} \frac{\gamma_i}{\lambda_{i}}\right)^{-1/2}, \quad \boldsymbol{\gamma} \in \mathbb{N}_0^{\dim(\mathcal{X})},\tag{15}\] and let \[\label{eq:32tau} \tau : \mathbb{N}\to \mathbb{N}_0^{\dim(\mathcal{X})}\tag{16}\] be a bijection that gives a nonincreasing rearrangement, i.e., \(v_{\boldsymbol{\tau(1)}} \geq v_{\boldsymbol{\tau(2)}} \geq \cdots > 0\). If \(\dim(\mathcal{X}) = \infty\), we replace \(\mathbb{N}_0^{\dim(\mathcal{X})}\) by \(\Gamma\). As we see in a moment, this rearranged sequence determines the error due to approximating \(F\) in a finite-dimensional subspace of Hermite polynomials.

4.2 Discussion of the error bound↩︎

Theorem 1 decomposes the error \(F - \hat{F}\) into four main terms \[\label{eq:32error32decomposition} {\left\|F - \widehat{F}\right \|}_{L_{\mu}^2(\mathcal{X}; \mathcal{Y})} \lesssim \mathrm{Err}_{\mathrm{true-PCA}} + \mathrm{Err}_{\mathrm{emp-PCA}} + \mathrm{Err}_{\mathrm{approx}} +\mathrm{Err}_{\mathrm{noise}}.\tag{17}\] These are, respectively, a true PCA projection error \[\label{eq:32true32PCA32projection32error} \mathrm{Err}_{\mathrm{true-PCA}} : = L \sqrt{\sum^{\dim(\mathcal{X})}_{i = d_{\mathcal{X}}+1} \lambda_i} + \sqrt{\sum^{\dim(\mathcal{Y})}_{i = d_{\mathcal{Y}}+1} \lambda^{F \sharp \mu}_i},\tag{18}\] an empirical PCA projection error \[\mathrm{Err}_{\mathrm{emp-PCA}} : = ({\left\|F(0)\right \|}_{\mathcal{Y}} + L (1+K_{\mu}) + \sigma ) \left [ \left ( \frac{d_{\mathcal{X}} \log(12/\epsilon)}{N_{\mathcal{X}}} \right )^{\frac{1}{4}} + \left ( \frac{d_{\mathcal{Y}} \log(12/\epsilon)}{N_{\mathcal{Y}}} \right )^{\frac{1}{4}} \right ]\] an approximation error \[\mathrm{Err}_{\mathrm{approx}} : = \frac{1+{\left\|F(0)\right \|}_{\mathcal{Y}} + L }{\sqrt{\epsilon}} \left( \frac{3}{2} \right)^k v^{k}_{\boldsymbol{\tau(s+1)}} {\left\| \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}} \right \|}_{H_{\hat{\varrho}}^{k}(\mathbb{R}^{d_{\mathcal{X}}}; \mathbb{R}^{d_{\mathcal{Y}}})}\] and a noise error \[\label{eq:32noise32error} \mathrm{Err}_{\mathrm{noise}} : = \frac{1+{\left\|F(0)\right \|}_{\mathcal{Y}} + L }{\sqrt{\epsilon}} \sqrt{\frac{s}{M}} \sigma .\tag{19}\]

4.2.0.1 The term \(\mathrm{Err}_{\mathrm{true-PCA}}\).

This is the error due to the truncation \(d_{\mathcal{X}} \leq \mathrm{dim}(\mathcal{X})\) and \(d_{\mathcal{Y}} \leq \mathrm{dim}(\mathcal{Y})\) incurred by considering finitely-many PCA terms. As expected from standard results [4], it behaves like the sum of the omitted PCA eigenvalues. Note that the eigenvalues \(\lambda_i^{F \sharp \mu}\) might exhibit very slow decay. However, if we know that \(F\) maps to some function space of higher regularity (e.g., \(H^1\)), then it can be shown that the PCA truncation error decays at least polynomially in \(d_{\mathcal{Y}}\), see [2].

4.2.0.2 The term \(\mathrm{Err}_{\mathrm{emp-PCA}}\).

This term accounts for computing the PCA eigenvalues and eigenfunctions from samples and scales like \((N_{\mathcal{X}})^{-1/4}\) and \((N_{\mathcal{Y}})^{-1/4}\), where \(N_{\mathcal{X}}\) is the amount of unlabeled samples \(\widehat{X}_i \sim \mu\) used to compute the PCA basis of \(\mathcal{X}\) and \(N_{\mathcal{Y}}\) is the amount of labeled samples \(\widehat{Y}_i = F(\widehat{X}_i) + \sigma \widehat{E}_i\) used to compute the PCA basis of \(\mathcal{Y}\). This scaling is expected from general bounds for empirical PCA, see [2].

Notice that \(N_{\mathcal{X}}\) and \(N_{\mathcal{Y}}\) must also satisfy ?? and ?? , respectively. The latter is quite reasonable, as it scales linearly in \(d_{\mathcal{Y}}\), the dimension of the PCA truncation. The former condition is more stringent. In particular, it suggests a scaling in \(s\) that is at least quartic. However, in practice (see Section 5) a log-linear scaling appears to suffice. We consequently believe this scaling is an artefact of our analysis. Improving it is a topic for future work.

4.2.0.3 The term \(\mathrm{Err}_{\mathrm{approx}}\).

This term is controlled by the parameter \(s\) and the weight \(v_{\boldsymbol{\tau(s+1)}}\), the former being related log-linearly to \(M\) via ?? . Put another way, with log-linear oversampling of \(s\) in the amount of labeled samples, one achieves an approximation rate specified by the behaviour \(v_{\boldsymbol{\tau(s+1)}}\) as \(s \rightarrow \infty\).

The precise decay rate of \(v_{\boldsymbol{\tau(s+1)}}\) is determined by the decay of the PCA eigenvalues \(\lambda_i\). We elaborate more on this behavior in the next section. A key point is that the approximation error behaves like \(v_{\boldsymbol{\tau(s+1)}}\) raised to the power \(k\), where \(k\) denotes the Sobolev regularity. This holds for any \(k\), thereby demonstrating the spectral approximation properties of Hermite-PCA approximation. Notice, in particular, that Algorithm 2 is independent of \(k\).

The sequence \(v_{\boldsymbol{\gamma}}\) and rearrangement \(\tau\) are closely related to the corresponding finite-dimensional quantities \(\hat{v}_{\boldsymbol{\gamma}}\) and \(\hat{\tau}\) defined in 12 and 13 that are used to construct the approximation \(\hat{F}\) via 14 . One can consider \(\hat{v}_{\boldsymbol{\gamma}}\) and \(\hat{\tau}\) as approximations to the ‘true’ terms \(v_{\boldsymbol{\gamma}}\) and \(\tau\) stemming from, firstly, the dimension truncation \(d_{\mathcal{X}} \leq \mathrm{dim}(\mathcal{X})\) and, secondly, the approximate computation of the true PCA eigenvalues \(\lambda_i\) via the empirical PCA eigenvalues \(\hat{\lambda}_i\). In particular, these quantities coincide when \(d_{\mathcal{X}} = \mathrm{dim}(\mathcal{X})\) and \(\hat{\lambda}_i = \lambda_i\), \(\forall i\). Crucially, Theorem 1 estimates the error in terms of the ‘true’ weight \(v_{\boldsymbol{\tau(s+1)}}\), with an additional additive term that accounts for the empirical PCA error. As a result, our method achieves near-optimal approximation rates for any level of empirical PCA error. We discuss the near-optimality of these rates further in Section 4.3 below.

4.2.0.4 The term \(\mathrm{Err}_{\mathrm{noise}}\).

This term scales like the noise standard deviation \(\sigma\) multiplied by the factor \(\sqrt{s/M}\). In particular, it tends to zero as \(M \rightarrow \infty\) for fixed \(s\), thereby demonstrating the denoising properties of Hermite-PCA approximation in the setting of statistical noise considered in this work. For other works addressing statistical noise in operator learning, see [43], [44].

4.3 Optimal approximation of Sobolev operators↩︎

As noted, in the absence of noise and empirical PCA error, the error \(F - \widehat{F}\) is determined by the quantity \(v^{k}_{\boldsymbol{\tau(s+1)}}\), which is itself controlled by the Sobolev regularity \(k\) and the true PCA eigenvalues \(\lambda_i\). It transpires that this quantity is fundamental: no algorithm that uses \(s\) linear measurements (which may or may not be pointwise evaluations) can achieve a faster rate of approximation in the \(L^2_{\mu}\)-norm uniformly for all Sobolev operators. More specifically, \(v^{k}_{\boldsymbol{\tau(s+1)}}\) is a lower bound for the (adaptive) \(s\)-width for the class of Sobolev operators \(H^k_{\mu}(\mathcal{X}; \mathcal{Y})\). We present this result in Appendix 11. It follows from adapting arguments in [39], which considered the case \(k = 1\) only.

A key consequence is that Hermite-PCA approximation achieves the optimal worst-case rates for learning Sobolev operators from \(M\) samples, up to the log-linear oversampling amount ?? . This also shows that pointwise samples, sampled randomly from an appropriate distribution, constitute near-optimal information for this problem. Moreover, as noted above, Hermite-PCA achieves this in a spectral fashion: it attains these near-optimal rates for any \(k\) and without any knowledge of \(k\).

4.4 Concrete rates of approximation↩︎

In applications, e.g., in the context of (functional) PCA [53], [54], the PCA eigenvalues \(\lambda_i\) are typically assumed to exhibit algebraic or exponential decay. The following bounds 20 , 21 follow from [39] in the case of infinite-dimensional domains. By Lemma 6, they hold for finite-dimensional domains as well. If \(\lambda_i = i^{- \alpha}\) for some \(\alpha > 1\), then for every \(\epsilon > 0\), \[\label{eq:32error32decay5932alg32EVs} v_{\boldsymbol{\tau(s+1)}} \lesssim \log(s)^{-\frac{\alpha}{2} + \epsilon}, \quad s \to \infty.\tag{20}\] If \(\lambda_i = e^{-\alpha i^{\beta}}\) for some \(\alpha, \beta > 0\), then for every \(\epsilon > 0\), \[\label{eq:32error32decay5932exp32EVs} v_{\boldsymbol{\tau(s+1)}} \lesssim e^{-\frac{1}{2} \alpha^{\frac{1}{\beta + 1}} (\beta' \log(s))^{1 / \beta'} + \epsilon}, \quad s \to \infty \quad \text{ with } \quad \beta' := 1 + \frac{1}{\beta}.\tag{21}\] These bounds yield estimates for the approximation error in terms of the approximation dimension \(s\). More concretely, if \(\lambda_i = i^{-\alpha}\), then \[\label{eq:32approx32error5932alg32EVs} \mathrm{Err}_{\mathrm{approx}} \lesssim_{F,k,\epsilon} \log(s)^{-\frac{k \alpha}{2} + \epsilon}, \quad s \to \infty\tag{22}\] subject to the log-linear scaling ?? . The underlying logarithmic decay is imposed by the algebraic decay of the eigenvalues, while the decay rate \(\alpha\) of the eigenvalues and Sobolev regularity \(k\) of the target operator determine the precise exponent: faster decay and/or higher Sobolev regularity ensure faster convergence. Similarly, if \(\lambda_i = e^{-\alpha i^{\beta}}\), then \[\label{eq:32approx32error5932exp32EVs} \mathrm{Err}_{\mathrm{approx}} \lesssim_{F,k,\epsilon} e^{-\frac{k}{2} \alpha^{\frac{1}{\beta + 1}} (\beta' \log(s))^{1 / \beta'} + \epsilon}, \quad s \to \infty.\tag{23}\] Note that in both cases, the decay is only subalgebraic in \(M\). In fact, this phenomenon is independent of the decay of the \(\lambda_i\) and constitutes an intrinsic curse of sample complexity on infinite-dimensional domains: Regardless of the decay rate of the PCA eigenvalues \(\lambda_i\), the quantity \(v_{\boldsymbol{\tau(s+1)}}\) can only decay subalgebraically as \(M \to \infty\). For further details, we refer to [39].

4.5 The regularity condition↩︎

Theorem 1 assumes the Sobolev regularity \(\widehat{\mathcal{E}}_{\mathcal{Y}}\circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}}\in H^k_{\hat{\varrho}}(\mathbb{R}^{d_{\mathcal{X}}} ; \mathbb{R}^{d_{\mathcal{Y}}} )\). This may be less than desirable, as it pertains to the latent space function \(\widehat{\mathcal{E}}_{\mathcal{Y}}\circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}}\) and the measure \(\hat{\varrho} = \mathcal{N}(0,\boldsymbol{\hat{\lambda}})\) is defined in terms of the empirical PCA eigenfunctions \(\hat{\lambda}_i\), as opposed to their true counterparts \(\lambda_i\). Ideally, one would consider Sobolev regularity \(F \in H^k_{\mu}(\mathcal{X}; \mathcal{Y})\) on the underlying operator \(F\) with respect to the underlying Gaussian measure \(\mu\). However, this does not appear straightforward within this framework, where regularity of the encoded function appears critically important to address the errors stemming from the approximation of the true PCA eigenvalues and eigenfunctions. It is a short argument to show that \(F \in H^k_{\mu}(\mathcal{X}; \mathcal{Y}) \not\Rightarrow \widehat{\mathcal{E}}_{\mathcal{Y}}\circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}}\in H^k_{\hat{\varrho}}(\mathbb{R}^{d_{\mathcal{X}}} ; \mathbb{R}^{d_{\mathcal{Y}}} )\), in general.

This raises the question of what kinds of regularity assumptions on \(F\) imply that \(\widehat{\mathcal{E}}_{\mathcal{Y}}\circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}}\in H^k_{\hat{\varrho}}(\mathbb{R}^{d_{\mathcal{X}}} ; \mathbb{R}^{d_{\mathcal{Y}}} )\). Fortunately, assuming a small amount of additional regularity is sufficient. In Appendix 12 we show that the class \(C^k_{\mu\text{-adm}}(\mathcal{X}; \mathcal{Y})\) consisting of \(C^k\)-operators whose derivatives do not grow too fast automatically satisfy \(\widehat{\mathcal{E}}_{\mathcal{Y}}\circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}}\in H^k_{\hat{\varrho}}(\mathbb{R}^{d_{\mathcal{X}}} ; \mathbb{R}^{d_{\mathcal{Y}}} )\). As a special case, this includes the class \(C^k_{\mathrm{Lip}}(\mathcal{X}; \mathcal{Y})\) of \(C^k\)-operators whose derivatives are Lipschitz continuous.

On the other hand, it remains an open problem to determine whether \(H^k_{\mu}(\mathcal{X}; \mathcal{Y})\)-regularity is sufficient to achieve optimal approximation rates in the practical setting studied in this paper: namely, where \(\mu\) is Gaussian but its covariance structure is unknown and therefore can only be estimated from samples.

4.6 Further discussion↩︎

Several other aspects of Theorem 1 warrant further discussion. First, notice that the result require \(F\) to be \(L\)-Lipschitz. This is a standard assumption used in encoding-decoding approaches to operator learning to bound errors stemming from approximate encoders and decoders. Second, Theorem 1 presents an error bound holding in high probability. One could also derive error bounds in expectation using standard techniques from the analysis of least-squares approximation from random samples. See [47]. Third, notice that the approximation and noise error terms scale like \(1/\sqrt{\epsilon}\) in the failure probability \(\epsilon\). This stems from the application of Markov’s inequality in the analysis of the least-squares problem 9 . While there are ways to mitigate this scaling [47], these generally require faster, algebraic and/or \(L^{\infty}\)-norm convergence of best approximations in \(\mathcal{P}\) to the target operator. Neither generally holds in the case of Sobolev operators, as elaborated previously. Fourth, we mention that there are ways to reduce the log-linear \(s\)-scaling of \(M\) to a linear scaling [47]. Incorporating them into the Hermite-PCA method is left to future work.

5 Numerical experiments↩︎

We present two different experiments which illustrate the validity of the theoretical error bound in Theorem 1 in practice. In the first experiment, we consider an obstacle problem on the one-dimensional unit interval. More specifically, we approximate the obstacle-to-solution operator associated to minimizing the Dirichlet energy of functions under a unilateral side constraint. We show how each of the error terms 18 –19 directly influences the overall approximation error via 17 , empirically validating our theoretical findings. In the second experiment, we demonstrate the spectral property of our algorithm by approximating Gaussian Sobolev functionals of varying regularity.

In each experiment, we use \(M = s \log(s)\) labeled training samples \((X_i, Y_i) \in \mathcal{X}\times \mathcal{Y}\) as in 3 , 4 for the least-squares fit and \(2000\) unseen labeled test samples from \(\mu \otimes \nu\) to obtain a Monte Carlo estimate for the relative test error \({\| F - \widehat{F}\|}_{L_{\mu}^2(\mathcal{X}; \mathcal{Y})} / {\left\| F \right \|}_{L_{\mu}^2(\mathcal{X}; \mathcal{Y})}.\) We run each experiment \(10\) times with different random seeds, but same test set, and plot the mean approximation error along with its one-standard-deviation-band (as shaded region) on a log-log-scale.

5.1 Obstacle problem in 1D↩︎

We consider the problem of minimizing the Dirichlet energy functional \(I(u) := \frac{1}{2} \int_0^1 (u'(x))^2 \,\mathrm{d}x\) among all functions \(u\) which belong to the set \(\mathcal{A}:= \{ w \in H_0^1(0,1) : w \geq v \text{ a.e. in } (0,1) \}\), where \(v \in H_0^1(0,1)\) is a given obstacle function. Equivalently, we seek the solution \(u \in \mathcal{A}\) of the variational inequality \[\int_0^1 u' (w' - u') \,\mathrm{d}x \geq 0, \quad \forall w \in \mathcal{A}.\] It is easy to see that such a solution exists and is unique. Hence, the obstacle-to-solution operator \[F : H_0^1(0,1) \to H_0^1(0,1), \quad v \mapsto F(v) = u,\] is well-defined. It is well-known from the regularity theory of obstacle problems that \(F\) is Lipschitz continuous, but has no higher global regularity, see [55], [56]. We can thus conclude that the latent space function \(\widehat{\mathcal{E}}_{\mathcal{Y}}\circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}}\) belongs to \(H_{\hat{\varrho}}^1(\mathbb{R}^{d_{\mathcal{X}}}; \mathbb{R}^{d_{\mathcal{Y}}})\), see Appendix 12, but can generally not expect higher Sobolev regularity.

Next, we describe our experimental setup. We parametrize the input and output functions by \(10\) sine basis functions in \(H_0^1(0,1)\). We solve the obstacle problem computationally by applying the projected Gauss-Seidel method on a mesh with \(257\) nodes to solve the corresponding Hamilton-Jacobi equation. We take \(\lambda_i \asymp i^{-2}\) and normalize, such that \(\sum_{i=1}^{10} \lambda_i = 1\). We apply both true and empirical PCA to encode (decode) the input (output) functions. With true PCA, there is no empirical projection PCA error, i.e., \(\mathrm{Err}_{\mathrm{emp-PCA}} = 0\), and we can individually study the influence of the remaining error terms \(\mathrm{Err}_{\mathrm{true-PCA}}\), \(\mathrm{Err}_{\mathrm{approx}}\), and \(\mathrm{Err}_{\mathrm{noise}}\) on the overall Hermite-PCA approximation error. In this case, the PCA basis functions in \(\mathcal{X}\) (\(\mathcal{Y}\)) are just given by the first \(d_{\mathcal{X}}\) (\(d_{\mathcal{Y}}\)) sine basis functions. In real-world applications, we of course have no full knowledge of the input and output functions and the true PCA basis functions. To model this situation, we apply empirical PCA and discretize the input and output functions via \(128\) piecewise linear finite element functions. See Figure 3 for evaluations of the learned obstacle-to-solution operator at a test obstacle.

True PCA: Approximation and noise error. We first set \(d_{\mathcal{X}} = d_{\mathcal{Y}} = 10\). In this case, there is no PCA projection error, i.e., \(\mathrm{Err}_{\mathrm{true-PCA}} = 0\), and the overall Hermite-PCA approximation error is controlled by \(\mathrm{Err}_{\mathrm{approx}} + \mathrm{Err}_{\mathrm{noise}}\). In the noiseless case (\(\sigma = 0\)) so that \(\mathrm{Err}_{\mathrm{noise}} = 0\), we expect by 22 the error term \(\mathrm{Err}_{\mathrm{approx}}\) to be of order \(\mathcal{O}(\log(s)^{-1})\). This is empirically confirmed by the results in Figure [fig:true-pca-approx]. In the presence of noise (\(\sigma > 0\)) our choice of \(M\) yields \(\mathrm{Err}_{\mathrm{noise}} \asymp \log(s)^{-1/2}\), see 19 , thus dominating \(\mathrm{Err}_{\mathrm{approx}}\). This is confirmed by Figure [fig:true-pca-noise], where we set \(\sigma = 0.2\) and chose the noise distribution to be the standard normal distribution \(\rho = \mathcal{N}(0, I_{10})\).

True PCA: Projection errors. The true PCA projection error term \(\mathrm{Err}_{\mathrm{true-PCA}}\), as described in 18 , is composed of a projection error in \(\mathcal{X}\), given by \(\sqrt{\sum_{i = d_{\mathcal{X}} + 1}^{\dim(\mathcal{X})} \lambda_i }\), and a projection error in \(\mathcal{Y}\), given by \(\sqrt{\sum_{i = d_{\mathcal{Y}} + 1}^{\dim(\mathcal{Y})} \lambda_i^{F \sharp \mu} }\). To study both terms separately, we set \(d_{\mathcal{X}} = 5\), \(d_{\mathcal{Y}} = 10\) and \(d_{\mathcal{X}} = 10\), \(d_{\mathcal{Y}} = 5\), respectively. To avoid any noise errors, we set \(\sigma = 0\). As \(\mathrm{Err}_{\mathrm{true-PCA}}\) is independent of \(s\), it dominates the Hermite-PCA approximation error for sufficiently large \(s\), which is reflected by the flattening of the error curves in Figure 5.

Empirical PCA. In order to study the empirical PCA projection error \(\mathrm{Err}_{\mathrm{emp-PCA}}\) and avoid other error sources, we again set \(d_{\mathcal{X}} = d_{\mathcal{Y}} = 10\) and \(\sigma = 0\). We check the hypothesis formulated in Section 4.2 that the at least quartic \(s\)-scaling for \(N_{\mathcal{X}}\) in ?? is pessimistic in practice. To this end, we construct the empirical PCA encoder and decoder using \(N_{\mathcal{X}} = d_{\mathcal{X}} s \log(s)\) samples from \(\mu\) and \(N_{\mathcal{Y}} = d_{\mathcal{Y}} \log(s)\) corresponding labels from \(F \sharp \mu\). The results in Figure 6 suggest that this log-linear \(s\)-scaling of \(N_{\mathcal{X}}\) and logarithmic \(s\)-scaling of \(N_{\mathcal{Y}}\) is sufficient in practice to obtain an empirical PCA projection error of order \(\mathcal{O}(\log(s)^{-1})\) – the same order as the approximation error term \(\mathrm{Err}_{\mathrm{approx}}\).

In summary, the findings in this section illustrate the excellent match between the empirically observed convergence of the Hermite-PCA approximation algorithm and our theoretical error bounds. In fact, the algorithm seemingly performs even better in practice – showing near-optimal scaling in all hyperparameters – than predicted by our theory.

5.2 Approximation of functionals↩︎

To empirically confirm the influence of the Sobolev order \(k\) on the approximation error, as predicted in 22 , 23 , we approximate parametric functionals whose parameter allows to control the Sobolev regularity. Let \(\varrho = \varrho_{\boldsymbol{\lambda}} := \bigotimes_{i = 1}^{d_{\mathcal{X}}} \mathcal{N}(0, \lambda_i)\). For \(\boldsymbol{a} \in \ell^2(\mathbb{N})\) and \(b > 0\) define the functional \[F_b : \ell^2(\mathbb{N}) \to \mathbb{R}, \quad F_b(\boldsymbol{x}) = \left|\left\langle \boldsymbol{a}, \boldsymbol{x} \right\rangle_{\ell^2}\right|^{b},\] as well as the input-truncated function \(f_b(\boldsymbol{x}) := |\sum_{i = 1}^{d_{\mathcal{X}}} a_i x_i|^b\) for \(\boldsymbol{x} \in \mathbb{R}^{d_{\mathcal{X}}}\). The next result shows how the parameter \(b\) influences the Sobolev regularity of \(f_b\).

In our experiments, we set \(a_i = i^{-1}\), \(b = k - \frac{1}{2} + 0.001\), and consider the cases \(k = 1, 2, 3\). We truncate the input sequences after \(10\) elements and use true PCA for the encoder in \(\mathcal{X}\) with \(d_{\mathcal{X}} = 10\). As \(d_{\mathcal{Y}} = 1\), there is no decoding in \(\mathcal{Y}\). To avoid noise errors, we set \(\sigma = 0\). Consequently, the only non-vanishing error term is the approximation error \(\mathrm{Err}_{\mathrm{approx}}\). We take \(\lambda_i \asymp e^{-i}\) and normalize such that \(\sum_{i=1}^{10} \lambda_i = 1\). By 21 , we expect the approximation error to be of order \(\mathcal{O}(e^{-k \sqrt{\log(s) / 2}})\). This is reflected by the results in Figure 7.

An interesting observation, most notable in Figures [fig:functional952], [fig:functional953], is the step-like decrease of the mean error curve. This is expected behavior, as the algorithm constructs the surrogate operator based on the index set \(S = \{ \boldsymbol{\hat{\tau}(1)}, \dots, \boldsymbol{\hat{\tau}(s)} \}\) corresponding to the \(s\) largest weights \(v_{\boldsymbol{\hat{\tau}(i)}}\). This choice of \(S\) is asymptotically optimal in the sense that it yields the optimal error \(v_{\boldsymbol{\tau(s + 1)}}^k\) as \(s \to \infty\). Locally, however, consecutive weights may have values close to each other. In this case, increasing \(s\) only leads to a slight decrease of the error and therefore to the formation of steps in the error curve.

Overall, the findings in this section confirm a major facet of the Hermite-PCA approximation algorithm: its spectral approximation property. The algorithm itself is independent of the Sobolev order \(k\), but it automatically yields faster convergence the smoother the objective operator is. The empirically observed decay rates match the rates predicted by our theory.

6 A general error bound for Hermite-PCA approximation↩︎

The next three sections of this paper are devoted to the proof of Theorem 1. First, in this section we establish a general error bound, Theorem 2, for the Hermite-PCA approximation defined in Algorithm 2, without making any regularity assumptions on \(F\) and also allowing for an arbitrary multi-index set \(S\) in the definition of the subspace \(\mathcal{P}\) in 7 . This derives the desired error bound, up to a best approximation error of the encoded function \(\widehat{\mathcal{E}}_{\mathcal{Y}} \circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}}\) in \(\mathcal{P}\). To analyze this term, we require an \(\ell^2\)-characterization of Gaussian Sobolev spaces in terms of Hermite polynomial coefficients, in tandem with several estimates for the decay of the weight sequences that arise in this characterization. This is the topic of Section 7. Finally, in Section 8 we combine these two pieces to establish Theorem 1.

6.1 Overview of the proof of Theorem 2↩︎

The remainder of this section is devoted to the proof of Theorem 2. By the triangle inequality, we can split the total error into three terms, \[\begin{align} \begin{aligned} \label{eq:32total32error32splitting} {\left\|\widehat{F}- F\right \|}_{L_{\mu}^2(\mathcal{X}; \mathcal{Y})} &\leq {\left\|\widehat{F} - \widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}} \circ \widehat{\mathcal{E}}_{\mathcal{X}}\right \|}_{L_{\mu}^2(\mathcal{X}; \mathcal{Y})} \\ &\quad + {\left\|\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}} \circ \widehat{\mathcal{E}}_{\mathcal{X}} - \widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F\right \|}_{L_{\mu}^2(\mathcal{X}; \mathcal{Y})} \\ &\quad + {\left\|\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F - F\right \|}_{L_{\mu}^2(\mathcal{X}; \mathcal{Y})} \\ &=: \mathrm{Err}_{\mathcal{A}} + \mathrm{Err}_{\mathcal{X}} + \mathrm{Err}_{\mathcal{Y}}. \end{aligned} \end{align}\tag{24}\] The right-hand side consists of the the approximation error \(\mathrm{Err}_{\mathcal{A}}\) and empirical PCA projection errors \(\mathrm{Err}_{\mathcal{X}}\) and \(\mathrm{Err}_{\mathcal{Y}}\) on \(\mathcal{X}\) and \(\mathcal{Y}\), respectively. Recall that \(F\) is \(L\)-Lipschitz and observe that the encoders and decoders are all \(1\)-Lipschitz by construction. Therefore, we can simplify these error terms as follows: \[\begin{align} &\mathrm{Err}_{\mathcal{X}} \leq L {\|\widehat{\mathcal{D}}_{\mathcal{X}} \circ \widehat{\mathcal{E}}_{\mathcal{X}} - \mathcal{I}_{\mathcal{X}}\|}_{L_{\mu}^2(\mathcal{X}; \mathcal{X})}, \tag{25} \\ &\mathrm{Err}_{\mathcal{Y}} = {\left\|\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} - \mathcal{I}_{\mathcal{Y}}\right \|}_{L_{F \sharp \mu}^2(\mathcal{Y}; \mathcal{Y})}, \tag{26} \\ &\mathrm{Err}_{\mathcal{A}} = {\|\widehat{F} - \widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}} \circ \widehat{\mathcal{E}}_{\mathcal{X}}\|}_{L_{\mu}^2(\mathcal{X}; \mathcal{Y})}. \tag{27} \end{align}\]

Here, \(\mathcal{I}_{\mathcal{X}}\), \(\mathcal{I}_{\mathcal{Y}}\) denote the identity operators on \(\mathcal{X}\) and \(\mathcal{Y}\), respectively. Since \(\widehat{\mathcal{E}}_{\mathcal{X}}\) is bounded and linear, \(\widehat{\mathcal{E}}_{\mathcal{X}} \sharp \mu\) in 27 is a Gaussian measure on \(\mathbb{R}^{d_{\mathcal{X}}}\).

In the next three subsections, we estimate these terms separately.

6.2 True PCA projection errors↩︎

Since \(\widehat{\mathcal{E}}_{\mathcal{X}} , \widehat{\mathcal{D}}_{\mathcal{X}}\) are approximations to the true PCA encoders and decoders \(\mathcal{E}_{\mathcal{X}} , \mathcal{D}_{\mathcal{X}}\), and likewise for \(\widehat{\mathcal{E}}_{\mathcal{Y}} , \widehat{\mathcal{D}}_{\mathcal{Y}}\), our first task is to study the true PCA projection errors. The former is straightforward. By [4], we have \[\label{eq:32encoding32projection32error} {\left\|\mathcal{D}_{\mathcal{X}}\circ \mathcal{E}_{\mathcal{X}}- \mathcal{I}_{\mathcal{X}}\right \|}_{L_{\mu}^2(\mathcal{X}; \mathcal{X})} = \sqrt{\sum_{i = d_{\mathcal{X}} + 1}^{\dim(\mathcal{X})} \lambda_i}.\tag{28}\] Analogously, \[\label{eq:32decoding32projection32error} {\|\mathcal{D}_{\mathcal{Y}}\circ \mathcal{E}_{\mathcal{Y}}- \mathcal{I}_{\mathcal{Y}}\|}_{L_{\nu}^2(\mathcal{Y}; \mathcal{Y})} = \sqrt{\sum_{i = d_{\mathcal{Y}} + 1}^{\dim(\mathcal{Y})} \lambda_i^{\nu}},\tag{29}\] where we recall from Section 2.1 that \(\nu\) is the distribution of \(Y = F(X) + \sigma E\), where \(X \sim \mu\) and \(E \sim \rho\) are mutually independent. However, the decoding error 26 is given in terms of the measure \(F \sharp \mu\), which only coincides with \(\nu\) in the noiseless case \(\sigma = 0\). To handle the noisy case, we require the following lemma.

Note that this bound is noise-consistent in the sense that if \(\sigma \to 0^+\), then we recover the noiseless PCA projection error in \(\mathcal{Y}\), given by \(\sqrt{\sum_{i = d_{\mathcal{Y}} + 1}^{\dim(\mathcal{Y})} \lambda_i^{F \sharp \mu}}\).

6.3 Empirical PCA projection errors↩︎

We now estimate \(\mathrm{Err}_{\mathcal{X}}\) and \(\mathrm{Err}_{\mathcal{Y}}\). For this we rely on the following result, which relates the empirical and true PCA projection errors. To state the result in suitable generality, let us consider a generic separable Hilbert space \(\mathcal{H}\) equipped with a subgaussian measure \(\tilde{\mu}\) with parameter \(K\). We define the true PCA encoders and decoders \(\mathcal{E}_{\mathcal{H}}, \mathcal{D}_{\mathcal{H}}\) analogously as in 5 . Furthermore, we draw \(n\) samples \(Z_i \sim_{\mathrm{i.i.d.}} \tilde{\mu}\), and define the empirical PCA encoder and decoder \(\widehat{\mathcal{E}}_{\mathcal{H}}, \widehat{\mathcal{D}}_{\mathcal{H}}\) with dimension \(d_{\mathcal{H}}\) analogously as in 6 .

We can thus estimate the empirical PCA projection errors with high probability via the bounds for the true PCA projection error, which we derived in Section 6.2, and sampling error terms depending on the number of samples \(N_{\mathcal{X}}\), \(N_{\mathcal{Y}}\) and the encoding and decoding dimension \(d_{\mathcal{X}}\), \(d_{\mathcal{Y}}\). We summarize this in the following result.

6.4 Error bounds for weighted least-squares approximation↩︎

We now consider the term \(\mathrm{Err}_{\mathcal{A}}\) defined in 27 . Since \(\widehat{f}\) is a solution of a weighted least-squares problem, in this subsection we develop some general results on weighted least-squares approximations to operators. In particular, we show how to choose the sampling measure \(\mu_{\mathrm{samp}}\) in such a way to obtain near-optimal sample complexity. It is based on the Christoffel function of a suitable chosen approximation space. The resulting Christoffel sampling method was originally developed in [45] for scalar-valued problems. See also [47] for extensions, variations and further topics.

Let \(\varpi\) be a probability measure on \(\mathcal{X}\). Fix \(s \in \mathbb{N}\) and let \(\mathcal{P}\subset L_{\varpi}^2(\mathcal{X}) \cap C(\mathcal{X})\) be a vector space of dimension \(s\) with orthonormal basis \(\{\xi_i\}_{i = 1}^{s}\). Given a vector subspace \(\mathcal{Y}' \subset \mathcal{Y}\), the \(\mathcal{Y}'\)-lift of \(\mathcal{P}\) is defined as the vector space \[\mathcal{P}_{\mathcal{Y}'} := \left \{ \sum_{i = 1}^s Y_i \xi_i : Y_i \in \mathcal{Y}' \right \} \subset L_{\varpi}^2(\mathcal{X}; \mathcal{Y}) \cap C(\mathcal{X}; \mathcal{Y}).\] Next, let \(w : \mathcal{X}\to (0, \infty)\) be a weight function, which we will specify later on, such that \(\int_{\mathcal{X}} 1 / w \,\mathrm{d}\varpi = 1\), and set \[\,\mathrm{d}\varpi_{\mathrm{samp}} := w^{-1} \,\mathrm{d}\varpi.\] Now draw \(X_1,\ldots,X_M \sim_{\mathrm{i.i.d.}} \varpi_{\mathrm{samp}}\) and \(E_1,\ldots,E_M \sim_{\mathrm{i.i.d.}} \rho\), where \(\rho\) satisfies Assumption 2 and set \(Y_i = F(X_i) + \sigma E_i\), \(i = 1,\ldots, M\). The weighted least-squares approximation of \(F\) in \(\mathcal{P}_{\mathcal{Y}'}\) is now given by \[\label{eq:32L-S32approximation} \widehat{F}= \widehat{F}(X_1, \dots, X_M) \in {\underset{P \in \mathcal{P}_{\mathcal{Y}'}}{\operatorname{argmin}}}{\frac{1}{M} \sum_{i = 1}^M w(X_i) {\left\|P(X_i) - Y_i\right \|}_{\mathcal{Y}}^2}.\tag{34}\] with \(X_i, Y_i\) as in 3 and 4 , respectively. The (reciprocal) Christoffel function of \(\mathcal{P}\) is defined as \[K(\mathcal{P}): \mathcal{X}\to \mathbb{R}, \qquad K(\mathcal{P})(X) := \sup\left\{ \frac{\left|p(X)\right|^2}{{\left\|p\right \|}_{L_{\varpi}^2(\mathcal{X})}^2} : p \in \mathcal{P}, p \neq 0 \right\},\] and we set \[\kappa_w (\mathcal{P}) := {\left\|w K(\mathcal{P})\right \|}_{L_{\varpi}^{\infty}(\mathcal{X})}.\] It can be readily checked that \(K(\mathcal{P}) = \sum_{i = 1}^s \xi_i^2\) for any orthonormal basis \(\{\xi_i\}_{i = 1}^s\) of \(\mathcal{P}\).

6.5 Application to the Hermite-PCA approximation↩︎

We now use Theorem 4 to bound \(\mathrm{Err}_{\mathcal{A}}\). Let \(S \subseteq \mathbb{N}^{d_{\mathcal{X}}}_0\) be an arbitrary index set of size \(|S| = s\). Then let \[\begin{align} \varpi & = \widehat{\mathcal{D}}_{\mathcal{X}} \sharp\hat{\varrho}, & \mathcal{P}&= \mathrm{span}\left\{ H_{\boldsymbol{\gamma}, \boldsymbol{\hat{\lambda}}} \circ \widehat{\mathcal{E}}_{\mathcal{X}}: \boldsymbol{\gamma} \in S \right\} \\ \mathcal{Y}' &= \widehat{\mathcal{D}}_{\mathcal{Y}}(\mathbb{R}^{d_{\mathcal{Y}}}), & w &= \left( \frac{1}{s} \sum_{\boldsymbol{\gamma} \in S} (H_{\boldsymbol{\gamma}, \boldsymbol{\hat{\lambda}}} \circ \widehat{\mathcal{E}}_{\mathcal{X}})^2 \right)^{-1}. \end{align}\] We now consider the weighted least-squares approximation \[\label{eq:32least-squares32for32hat-F} \widehat{F} \in {\underset{P \in \mathcal{P}_{\mathcal{Y}'}}{\operatorname{argmin}}}{\frac{1}{M} \sum^{M}_{i=1} w(X_i) {\left\|P(X_i) - \widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}}(Y_i)\right \|}_{\mathcal{Y}}^2}.\tag{36}\] Notice that this formulation uses the projected samples \(\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}}(Y_i)\). This is crucial, in that it allows us to relate \(\widehat{F}\) to the latent space approximation \(\widehat{f}\).

We now seek to apply Theorem 4. For this, we require two observations. First, notice that \[\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}}(Y_i) = \widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F(X_i) + \sigma \widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}}(E_i)\] and the noise terms \(\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}}(E_i)\) are subgaussian with parameter \(K = 1\) since the random variable \(E_i\) are subgaussian with \(K_{\rho} = 1\) and the map \(\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}}\) is an orthogonal projection. Second, observe that the functions \(\{ H_{\boldsymbol{\gamma} , \boldsymbol{\hat{\lambda}} } \circ \widehat{\mathcal{E}}_{\mathcal{X}} \}_{\boldsymbol{\gamma} \in S }\) are an orthonormal basis for \(\mathcal{P}\) in \(L^2_{\varpi}(\mathcal{X})\), where we recall that \(\varpi = \widehat{\mathcal{D}}_{\mathcal{X}} \sharp\hat{\varrho}\). Indeed, using the fact that \(\widehat{\mathcal{E}}_{\mathcal{X}} \circ \widehat{\mathcal{D}}_{\mathcal{X}} = \mathcal{I}_{\mathbb{R}^{d_{\mathcal{X}}}}\) once more, we have \[\begin{align} \int_{\mathcal{X}} H_{\boldsymbol{\gamma} , \boldsymbol{\hat{\lambda}} } \circ \widehat{\mathcal{E}}_{\mathcal{X}} (X) H_{\boldsymbol{\gamma'} , \boldsymbol{\hat{\lambda}} } \circ \widehat{\mathcal{E}}_{\mathcal{X}}(X) \,\mathrm{d}\widehat{\mathcal{D}}_{\mathcal{X}} \sharp\hat{\varrho}(X) = \int_{\mathbb{R}^{d_{\mathcal{X}}}} H_{\boldsymbol{\gamma} , \boldsymbol{\hat{\lambda}} } (\boldsymbol{x}) H_{\boldsymbol{\gamma'} , \boldsymbol{\hat{\lambda}} } (\boldsymbol{x}) \,\mathrm{d}\hat{\varrho}(\boldsymbol{x}) = \delta_{\boldsymbol{\gamma},\boldsymbol{\gamma'}}. \end{align}\] Therefore \(w^{-1}\) is, up to the scalar multiple \(s\), precisely the reciprocal Christoffel function of \(\mathcal{P}\). This implies that \(\kappa_{w}(\mathcal{P}) = s\) and therefore Theorem 4 with \(F\) replaced by \(\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F\) gives that \[\begin{align} {\left\|\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F - \widehat{F} \right \|}_{L^2_{\varpi}(\mathcal{X}; \mathcal{Y})} \leq & \left ( 1 + \sqrt{\frac{3 s}{M \epsilon}} \frac{1}{1-\delta} \right ) \inf_{P \in \mathcal{P}_{\mathcal{Y}'}} {\left\|\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F - P \right \|}_{L^2_{\varpi}(\mathcal{X}; \mathcal{Y})} + \sigma \sqrt{\frac{6 s}{M \epsilon} } \frac{1}{1-\delta}, \end{align}\] with probability at least \(1-\epsilon\), provided \(M \geq C_{\delta} s \log(6 s / \epsilon).\) To relate this to the error \(\mathrm{Err}_{\mathcal{A}}\), we use the fact that \(\widehat{\mathcal{E}}_{\mathcal{X}} \circ \widehat{\mathcal{D}}_{\mathcal{X}} = \mathcal{I}_{\mathbb{R}^{d_{\mathcal{X}}}}\) for a third time and write \[\begin{align} {\left\|\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F - \widehat{F} \right \|}_{L^2_{\varpi}(\mathcal{X}; \mathcal{Y})} & = {\left\|\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}} - \widehat{F} \circ \widehat{\mathcal{D}}_{\mathcal{X}} \right \|}_{L^2_{\hat{\varrho}}(\mathbb{R}^{d_{\mathcal{X}}} ; \mathcal{Y})} \\ & = {\left\|\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}} \circ \widehat{\mathcal{E}}_{\mathcal{X}} \circ \widehat{\mathcal{D}}_{\mathcal{X}} - \widehat{F} \circ \widehat{\mathcal{D}}_{\mathcal{X}} \right \|}_{L^2_{\hat{\varrho}}(\mathbb{R}^{d_{\mathcal{X}}} ; \mathcal{Y})} \\ & = {\left\|\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}} \circ \widehat{\mathcal{E}}_{\mathcal{X}} - \widehat{F} \right \|}_{L^2_{\varpi}(\mathcal{X}; \mathcal{Y})} \end{align}\] Thus we get \[\label{eq:32least32squares32error32bound} \begin{align} {\left\|\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}} \circ \widehat{\mathcal{E}}_{\mathcal{X}} - \widehat{F} \right \|}_{L^2_{\varpi}(\mathcal{X}; \mathcal{Y})} \leq & \left ( 1 + \sqrt{\frac{3 s}{M \epsilon}} \frac{1}{1-\delta} \right ) \inf_{P \in \mathcal{P}_{\mathcal{Y}'}} {\left\|\widehat{\mathcal{D}}_{\mathcal{Y}} \circ \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F - P \right \|}_{L^2_{\varpi}(\mathcal{X}; \mathcal{Y})} \\ & + \sigma \sqrt{\frac{6 s}{M \epsilon} } \frac{1}{1-\delta}, \end{align}\tag{37}\] with probability at least \(1-\epsilon\), provided \(M \geq C_{\delta} s \log(6 s / \epsilon)\). However, there remains an issue. The left-hand side is close to the approximation error \(\mathrm{Err}_{\mathcal{A}}\) defined in 27 , but the measure is \(\varpi = \widehat{\mathcal{D}}_{\mathcal{X}} \sharp\hat{\varrho}\) as opposed to \(\mu\). In the next subsection, we perform the necessary switch of measure.

6.6 Bounding the \(L^2_{\mu}(\mathcal{X}; \mathcal{Y})\)-norm error in terms of the \(L^2_{\varpi}(\mathcal{X}; \mathcal{Y})\)-norm error↩︎

We now show how a change of measure from \(\mu\) to \(\varpi\) affects the overall error bound. The following lemma assumes closeness of the true PCA eigenvalues \(\lambda_i\) and their approximate counterparts \(\hat{\lambda}_i\). Later, we quantify the error between the two, which leads to the main result of this subsection, Theorem 6. In the following and throughout, \({\left\| \cdot \right \|}_{\infty}\) denotes the operator norm and \(\Gamma(s,x) = \int_{x}^{\infty} t^{s-1} e^{-t} \,\mathrm{d}t\) is the upper incomplete Gamma function.

Lemma 3 yields a bound for switching between the measure \(\mu\) and the measure \(\varpi\), assuming sufficient closeness of the true and approximate PCA eigenvalues, as quantified by the uniform norm of the term \(\Delta = \Sigma - \widehat{\Sigma}\). We now estimate this quantity, which leads to a more explicit bound. For this, we require the following theorem.

Next, we notice that Lemma 3 introduces two parameters \(\kappa\) and \(R\). The final step in this subsection is to optimize these parameters. For this we require the following two lemmas.

With this in hand, we can now present the main result of this subsection.

6.7 Final arguments↩︎

We are now, finally, ready to finish the proof of Theorem 2. This involves combining Theorem 3, which bounds the empirical PCA projection errors 25 and 26 , with Theorem 6, which bounds the approximation error 27 .

7 \(\ell^2\)-characterization of Gaussian Sobolev spaces↩︎

In this and the next section, we aim to use Theorem 2 to prove the main result, Theorem 1. The key step is the analysis of the best approximation error \[\inf_{p \in \mathcal{P}} {\left\|\widehat{\mathcal{E}}_{\mathcal{Y}} \circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}}- p \right \|}_{L^2_{\hat{\varrho}}(\mathbb{R}^{d_{\mathcal{X}}} ; \mathbb{R}^{d_{\mathcal{Y}}})}\] and, in tandem, the derivation of the choice 14 for the index set \(S\). This analysis relies crucially on, firstly, the \(\ell^2\)-characterization of the Sobolev spaces \(H^k_{\varrho}(\mathbb{R}^d ; \mathbb{R}^{d'} )\) in terms of the Hermite polynomial coefficients and, secondly, analysis of the decay of the associated weight sequences that arise in the former characterization. We tackle these aspects of the analysis in this section, before completing the proof of Theorem 1 in the next.

We fix dimensions \(d, d' \in \mathbb{N}\). Recall the notation \([d] = \{1, \dots, d\}\) and \([\infty] = \mathbb{N}\), and, as before, let \(\boldsymbol{\lambda} = (\lambda_i)_{i=1}^d\) and \(\varrho = \varrho_{\boldsymbol{\lambda}} := \mathcal{N}(0, \boldsymbol{\lambda})\) be a centered, nondegenerate Gaussian measure on \(\mathbb{R}^{d}\) with diagonal covariance matrix with diagonal entries \(\lambda_i\).

7.1 \(\ell^2\)-characterization of \(H^k_{\varrho}(\mathbb{R}^d ; \mathbb{R}^{d'} )\)↩︎

We first require some further notation. For \(i \in \mathbb{N}\) and \(\boldsymbol{\gamma} \in \mathbb{N}_0^{d}\), we set \(\boldsymbol{\gamma^{(i)}} = \boldsymbol{0}\) if \(\gamma_i = 0\). If \(\gamma_i > 0\), we set \[\gamma_k^{(i)} = \begin{cases} \gamma_k - 1 & \mathrm{ if } k = i, \\ \gamma_k & \mathrm{ otherwise}. \end{cases}\] Recursively, we define \(\boldsymbol{\gamma^{(i,j)}} := \boldsymbol{(\gamma^{(i)})^{(j)}}\). For \(\alpha_1, \dots, \alpha_n \in \mathbb{N}\), we set \(\boldsymbol{\gamma^{(\alpha_1, \dots, \alpha_{j-1})}} := \boldsymbol{\gamma}\) if \(j = 1\). Note that we can extend this definition to \(d = \infty\) by replacing \(\mathbb{N}_0^d\) with \(\Gamma\).

Next, fix a multiindex \(\boldsymbol{\alpha} \in \mathbb{N}^n\). By ?? , the weight \(\prod_{j = 1}^n \gamma_{\alpha_j}^{(\alpha_1, \dots, \alpha_{j-1})} / \lambda_{\alpha_j}\) corresponds to the derivative \(D^{\boldsymbol{\alpha}} f\). In view of the definition of the \(H_{\varrho}^{k}\)-norm in Appendix 10, we define the weights \[\label{eq:32weight32formula} v_{\boldsymbol{\gamma}, (k), d} := \left( 1 + \sum_{j = 1}^k \sum_{\boldsymbol{\alpha} \in [d]^j} \prod_{i = 1}^j \frac{\gamma_{\alpha_i}^{(\alpha_1, \dots, \alpha_{i-1})}}{\lambda_{\alpha_i}} \right)^{-1/2}\tag{53}\] for \(k \in \mathbb{N}\) and \(\boldsymbol{\gamma} \in \mathbb{N}_0^{d}\). We extend this definition to \(d = \infty\) in which case we replace \(\mathbb{N}_0^d\) by \(\Gamma\).

7.2 Weight formula↩︎

Next, we derive an equivalent expression for the weights \(v_{\boldsymbol{\gamma}, (k), d}\), which we will need later to derive upper bounds in Proposition 1, but which might be of independent interest as well. To this end, let \(\Pi(n)\) denote the set of partitions of \([n]\), \(n \in \mathbb{N}\), and set \(\Pi(0) := \emptyset\).

7.3 Upper bounds for the weights↩︎

We now derive an upper bound for the decay of the weights \(v_{\boldsymbol{\gamma}, (k), d}\) in 53 . We denote by \(\tau_{(k),d} : \mathbb{N}\to \mathbb{N}_0^{d}\) a bijection that gives a nonincreasing rearrangement of the weights \(v_{\boldsymbol{\gamma}, (k), d}\), i.e., \(v_{\boldsymbol{\tau_{(k), d}(1)}, (k), d} \geq v_{\boldsymbol{\tau_{(k), d}(2)}, (k), d} \geq \cdots\). This mapping is unique up to permutations of weights of the same value. To avoid notational redundancy, we write \(v_{\boldsymbol{\tau_{(k), d}(i)}, (k), d} = v_{\boldsymbol{\tau_{(k), d}(i)}}\). For brevity, we further introduce the following notational conventions. If \(d = \dim(\mathcal{X})\), we usually write \(v_{\boldsymbol{\gamma}, (k), \dim(\mathcal{X})} = v_{\boldsymbol{\gamma}, (k)}\) and \(\tau_{(k), \dim(\mathcal{X})} = \tau_{(k)}\). If \(k = 1\), we usually write \(v_{\boldsymbol{\gamma}, (1), d} = v_{\boldsymbol{\gamma}, d}\) and \(\tau_{(1), d} = \tau_{d}\). In particular, we usually write \(\tau_{(1), \dim(\mathcal{X})} = \tau\). If \(d = \infty\), the same notation applies with \(\mathbb{N}_0^d\) replaced by \(\Gamma\). Note that this notation is consistent with 15 , 16 . If the \(\lambda_i\) are replaced by their empirical versions \(\hat{\lambda}_i\), we define the weights \(\hat{v}_{\boldsymbol{\gamma}, (k), d}\) together with a nonincreasing rearrangement \(\hat{\tau}_{(k), d}\) analogously. The only difference is that we write \(\hat{v}_{\boldsymbol{\gamma}, (k), d_{\mathcal{X}}} = \hat{v}_{\boldsymbol{\gamma}, (k)}\) and \(\hat{\tau}_{(k), d_{\mathcal{X}}} = \hat{\tau}_{(k)}\). That is, we drop the dimension \(d\) in the index if \(d= d_{\mathcal{X}}\) instead of \(d = \dim(\mathcal{X})\) as before. Note that this is consistent with the notation in 12 , 13 .

We proceed with two lemmas that we will later need in Proposition 1 to bound the best approximation error in terms of the weight \(v_{\boldsymbol{\tau(s + 1)}}^k\).

Finally, we conclude this section by showing that in finite dimensions the weights to Sobolev order \(k\) are bounded by the \(k\)th power of corresponding weights to Sobolev order one. Crucially, this holds for the rearranged sequences subject to the nonincreasing rearrangement \(\tau_{(1),d}\) corresponding to the order-one Sobolev weights. This is a key step in establishing the spectral convergence properties of the Hermite-PCA approximation, as it facilitates the choice 14 of the index set \(S\) that is independent of the Sobolev order \(k\).

8 Proof of Theorem 1↩︎

In this section, we prove the main result of the paper, Theorem 1.

8.1 Analysis of the best approximation error↩︎

With the results of the previous section in hand, we first bound the best approximation error \[\inf_{p \in \mathcal{P}} {\left\|\widehat{\mathcal{E}}_{\mathcal{Y}} \circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}}- p \right \|}_{L^2_{\hat{\varrho}}(\mathbb{R}^{d_{\mathcal{X}}} ; \mathbb{R}^{d_{\mathcal{Y}}})}\] where \(\mathcal{P}= \mathcal{P}_{\mathbb{R}^{d_{\mathcal{Y}}}} := \left \{ \sum_{\boldsymbol{\gamma} \in S} \boldsymbol{c}_{\boldsymbol{\gamma}} H_{\boldsymbol{\gamma},\boldsymbol{\hat{\lambda}}} : \boldsymbol{c}_{\boldsymbol{\gamma}} \in \mathbb{R}^{d_{\mathcal{Y}}} \right \} \subseteq L^2_{\hat{\varrho} } (\mathbb{R}^{d_{\mathcal{X}} } ; \mathbb{R}^{d_{\mathcal{Y}}} )\) and \(S = \hat{\tau}([s])\) are as in 7 and 14 , respectively, and \(\hat{\tau}\) is as in 13 . Crucially, this lemma bounds the error in terms of the true weights \(v_{\boldsymbol{\gamma}}\).

8.2 Final arguments↩︎

9 Conclusions↩︎

In this article, we presented a fully data-driven algorithm, termed Hermite-PCA approximation, for the recovery of Sobolev operators from pointwise, noisy samples. It is based on PCA for dimension reduction and employs weighted linear least-squares fitting, making it computationally efficient. It has several important properties. First, it achieves near-optimal sample complexity under weak assumptions – we only assume the input measure \(\mu\) to be Gaussian. Second, it does so in a spectral fashion, achieving faster convergence rates the smoother the objective operator is. Third, it can be adapted to any other regularity, e.g., mixed regularity, while maintaining its properties. We present a full error analysis of our algorithm, thus allowing for optimal choices of all hyperparameters. Alongside, we provide numerical results which closely match our theoretical findings. In particular, to the best of our knowledge, our experiments numerically illustrate for the first time the curse of sample complexity which is inherent to the approximation of finitely regular operators.

There are several future research directions. First, our algorithm requires a custom sampling distribution \(\mu_{\mathrm{samp}}\) of the input training data. It is based on the Christoffel function of the learning problem and is thus intrinsic to using Hermite polynomials [31]. In practice, Monte Carlo sampling from the underlying Gaussian input distribution \(\mu\) would be more favorable. It is an open question whether our error bounds also hold in the case of i.i.d. samples from \(\mu\).

Second, the bound ?? indicates an at least quartic \(s\)-scaling of the number \(N_{\mathcal{X}}\) of unlabeled data samples used in the empirical PCA encoding step. We believe that this is an artifact of our analysis. In fact, our empirical results suggest that it can be significantly improved to an at least log-linear scaling. This is topic of future work.

Third, in Appendix 12, we identify conditions for the operator \(F\) which imply its latent space representation to be Sobolev regular. It can be shown that \(H_{\mu}^k(\mathcal{X}; \mathcal{Y})\)-regularity of \(F\) itself is not sufficient for the analysis of the algorithm presented in this work. It is an open question whether the Hermite-PCA algorithm can attain optimal rates for all \(H_{\mu}^k(\mathcal{X}; \mathcal{Y})\)-operators based on only finitely many samples from \(\mu\), or if different algorithmic approaches are required. We mention [65] for recent results in finite dimensions.

Finally, as noted above, our algorithm can be readily adapted to other types of regularity, such as anisotropic or mixed smoothness or otherwise. An in-depth study of operator learning under other smoothness types is an interesting topic for future work.

Acknowledgments↩︎

BA acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC) through grants RGPIN/2026-04531. MG and GM acknowledge support from the Hausdorff Center for Mathematics (HCM) in Bonn, funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2047/2 – 390685813.

10 Definition of \(H^k_{\mu}(\mathcal{X}; \mathcal{Y})\)↩︎

In this appendix, we present a formal definition of the Sobolev spaces \(H^k_{\mu}(\mathcal{X}; \mathcal{Y})\) along with several results needed in the main text. We commence with some preliminary notions.

Recall that \(\mathrm{HS}_k(\mathcal{X}; \mathcal{Y})\) denotes the set of all \(k\)-linear operators \(F : \mathcal{X}^k \to \mathcal{Y}\) with finite Hilbert-Schmidt norm \[{\left\|F\right \|}_{\mathrm{HS}_k(\mathcal{X}; \mathcal{Y})} = \left( \sum_{i_1, \dots, i_k = 1}^{\infty} {\left\|F(\boldsymbol{e_{i_1}}, \dots, \boldsymbol{e_{i_k}})\right \|}_{2}^2 \right)^{1/2},\] where \(\boldsymbol{e_i}\) denotes the \(i\)th standard unit vector in \(\mathbb{N}_0^{\dim(\mathcal{X})}\) if \(\dim(\mathcal{X}) < \infty\) and in \(\Gamma\) if \(\dim(\mathcal{X}) = \infty\).

The weak Gaussian derivatives \(D^j F\) of \(F \in H_{\mu}^{k}(\mathcal{X}; \mathcal{Y})\) are well-defined for \(1 \leq j \leq k\) in the sense that if two sequences \(\{\varphi_i\}_i, \{ \widetilde{\varphi}_i\}_i\) from \(\mathcal{F}C_b^k(\mathcal{X}; \mathcal{Y})\) are Cauchy in \(H_{\mu}^{k}(\mathcal{X}; \mathcal{Y})\) and converge to \(F\) in \(L_{\mu}^2(\mathcal{X}; \mathcal{Y})\), then the sequences \(\{ D^j \varphi_i \}_i\), \(\{ D^j \widetilde{\varphi}_i \}_i\) have equal limits (denoted by \(D^j F\)) in \(L^2_{\mu}(\mathcal{X}; \mathrm{HS}_j(\mathcal{X}; \mathcal{Y}))\). We refer to [66] for details. We remark that our definition of \(H_{\mu}^{k}(\mathcal{X}; \mathcal{Y})\) deviates from the one therein in that we consider the Fréchet derivative along the full space \(\mathcal{X}\) whereas Bogachev considers the derivative along the Cameron-Martin space \(H(\mu)\) of \(\mu\). The proof for well-definedness of the weak derivatives, however, is analogous for both cases. If \(\dim(\mathcal{X}) < \infty\), then both spaces coincide, as in this case \(H(\mu) = \mathcal{X}\).

In finite dimensions, Gaussian Sobolev functions can be characterized by classical Sobolev functions. To this end, as in the main text, write \(\boldsymbol{\lambda} = (\lambda_i)_{i = 1}^d\), \(d \in \mathbb{N}\), with \(\lambda_i > 0\), and let \(\varrho = \varrho_{\boldsymbol{\lambda}} := \mathcal{N}(0, \boldsymbol{\lambda})\) denote a centered, nondegenerate Gaussian measure on \(\mathbb{R}^{d}\) with diagonal covariance matrix with diagonal entries \(\lambda_i\).

11 Optimal approximation of Gaussian Sobolev operators↩︎

We show that the worst-case approximation error for \(H_{\mu}^k(\mathcal{X}; \mathcal{Y})\)-operators based on \(s\) potentially adaptively chosen linear samples is bounded from below by the quantity \(v_{\boldsymbol{\tau(s + 1)}}^k\). To this end, we first recall the notion of the adaptive \(s\)-width and then adapt arguments from [39].

The adaptive \(s\)-width of a set \(\mathcal{K}\) describes the smallest worst-case error that can be achieved when we reconstruct all operators in \(\mathcal{K}\) by a reconstruction mapping \(\mathcal{T}\) from \(s\) samples that have been generated by an adaptive Hilbert-valued sampling operator \(\mathcal{L}\). We are interested in a specific choice for \(\mathcal{K}\), namely the \(k\)-Sobolev unit ball \[B_{\mu}^k (\mathcal{X}; \mathcal{Y}) := \left\{ F \in H_{\mu}^k(\mathcal{X}; \mathcal{Y}): {\left\|F\right \|}_{H_{\mu}^k(\mathcal{X}; \mathcal{Y})} \leq 1 \right\}.\] We write \(\Gamma_d = \mathbb{N}_0^d\), \([d] = \{1, \dots, d\}\) for \(d \in \mathbb{N}\) and \(\Gamma_{\infty} = \Gamma\), \([\infty] = \mathbb{N}\). Recall from 53 the weights \[v_{\boldsymbol{\gamma}, (k), d} := \left( 1 + \sum_{j = 1}^k \sum_{\boldsymbol{\alpha} \in [d]^j} \prod_{i = 1}^j \frac{\gamma_{\alpha_i}^{(\alpha_1, \dots, \alpha_{i-1})}}{\lambda_{\alpha_i}} \right)^{-1/2}, \quad k\in \mathbb{N}, d \in \mathbb{N}\cup \{ \infty \}, \boldsymbol{\gamma} \in \Gamma_d.\] together with a nonincreasing rearrangement \(\tau_{(k),d} : \mathbb{N}\to \Gamma_d\). Also recall the notational conventions introduced in the beginning of Section 7.3.

For \(k = 1\) and \(\dim(\mathcal{X}) = \infty\), this result follows from Theorem 5.4 in [39]. The proof holds verbatim in the case \(\dim(\mathcal{X}) < \infty\) as well. It is based on defining a suitable operator in the \(1\)-Sobolev unit ball which allows one to reduce the continuous approximation problem over \(H_{\mu}^1\)-operators to a discrete approximation problem over finite weighted \(\ell^2\)-sequences. The error in the latter can then be analyzed explicitly via the Kolmogorov width. The generalization of the argument to arbitrary Sobolev orders \(k \in \mathbb{N}\) is essentially the same. See also [65] for a proof in the general case for scalar-valued functions. For these reasons, we omit the proof of Theorem 7.

The next result relates the weight \(v_{\boldsymbol{\tau_{(k), d}(s)}}\) to the weight \(v_{\boldsymbol{\tau_{(1), d}(s)}}\). For \(k \in \mathbb{N}\), \(d \in \mathbb{N}\cup \{\infty\}\), and \(T > 0\), we define the set \[S_{k, d}(T) := \left\{ \boldsymbol{\gamma} \in \Gamma_d : v_{\boldsymbol{\gamma}, (k), d}^{-2} \leq 1 + T \right\} = \left\{ \boldsymbol{\gamma} \in \Gamma_d : \sum_{j = 1}^k \sum_{\boldsymbol{\alpha} \in [d]^j} \prod_{i = 1}^j \frac{\gamma_{\alpha_i}^{(\alpha_1, \dots, \alpha_{i-1})}}{\lambda_{\alpha_i}} \leq T \right\}.\]

In summary, Theorem 7 and Proposition 1 constitute a lower bound for the best approximation error which matches the upper bound in Theorem 1 up to constants.

12 \(C^k\)-operators with admissible growth at infinity↩︎

We introduce a set of operators \(F\), for which the latent space function \(\widehat{f} := \widehat{\mathcal{E}}_{\mathcal{Y}} \circ F \circ \widehat{\mathcal{D}}_{\mathcal{X}}\) belongs to \(H_{\hat{\varrho}}^{k}(\mathbb{R}^{d_{\mathcal{X}}}; \mathbb{R}^{d_{\mathcal{Y}}})\) and the norm can be bounded independently of the empirical PCA eigenvalues \(\hat{\lambda}_i\). This can be achieved if one has control over the behavior of \(F\) which prevents if from growing too fast at infinity. We then discuss relevant subsets of such operators, most importantly the set of \(C^k\)-operators whose derivatives are Lipschitz continuous.

12.1 Definitions and properties↩︎

Let \(F \in C_{\mu\mathrm{-adm}}^k(\mathcal{X}; \mathcal{Y})\) with corresponding \(\mu\)-admissible growth function \(\omega\). Repeated application of the chain rule yields \[D^k \widehat{f}(\boldsymbol{x})(\boldsymbol{z_1}, \dots, \boldsymbol{z_k}) = \widehat{\mathcal{E}}_{\mathcal{Y}} \left( D^k F(\widehat{\mathcal{D}}_{\mathcal{X}}(\boldsymbol{x})) (\widehat{\mathcal{D}}_{\mathcal{X}}(\boldsymbol{z_1}), \dots, \widehat{\mathcal{D}}_{\mathcal{X}}(\boldsymbol{z_k}) \right), \quad \boldsymbol{z_1}, \dots, \boldsymbol{z_k} \in \mathbb{R}^{d_{\mathcal{X}}}\] for all points \(\boldsymbol{x} \in \mathbb{R}^{d_{\mathcal{X}}}\) where this derivative exists. Note that the set of these points is a \(\widehat{\mathcal{E}}_{\mathcal{X}}\sharp \mu\)-null set in \(\mathbb{R}^{d_{\mathcal{X}}}\). Since \(\widehat{\mathcal{E}}_{\mathcal{Y}}\) and \(\widehat{\mathcal{D}}_{\mathcal{X}}\) are contractions and \(\omega\) is monotone, we may compute for \(j \in [k]\) \[\begin{gather} \sum_{i_1, \dots, i_j = 1}^{d_{\mathcal{X}}} {\left\| D^j f(\boldsymbol{x})(\boldsymbol{e_{i_1}}, \dots, \boldsymbol{e_{i_j}}) \right \|}_2^2 =\sum_{i_1, \dots, i_j = 1}^{d_{\mathcal{X}}} {\left\| \widehat{\mathcal{E}}_{\mathcal{Y}} \left( D^j F(\widehat{\mathcal{D}}_{\mathcal{X}}(\boldsymbol{x})) (\widehat{\mathcal{D}}_{\mathcal{X}}(\boldsymbol{e_{i_1}}), \dots, \widehat{\mathcal{D}}_{\mathcal{X}}(\boldsymbol{e_{i_j}})) \right) \right \|}_2^2 \\ \leq \sum_{i_1, \dots, i_j = 1}^{d_{\mathcal{X}}} {\left\|D^j F(\widehat{\mathcal{D}}_{\mathcal{X}}(\boldsymbol{x})) (\widehat{\phi}_{i_1}, \dots, \widehat{\phi}_{i_j})\right \|}_{\mathcal{Y}}^2 \leq \omega({\|\widehat{\mathcal{D}}_{\mathcal{X}}(\boldsymbol{x})\|}_{\mathcal{X}})^2 \leq \omega({\left\|\boldsymbol{x}\right \|}_2)^2. \end{gather}\] This shows that \(\widehat{f} \in C_{\widehat{\mathcal{E}}_{\mathcal{X}} \sharp \mu \mathrm{-adm}}^k(\mathbb{R}^{d_{\mathcal{X}}}; \mathbb{R}^{d_{\mathcal{Y}}})\), provided that \(\omega\) is \(\widehat{\mathcal{E}}_{\mathcal{X}} \sharp \mu\)-admissible. In this case it follows from Proposition 1 that \(\widehat{f} \in H_{\hat{\varrho}}^{k}(\mathbb{R}^{d_{\mathcal{X}}}; \mathbb{R}^{d_{\mathcal{Y}}})\) and the corresponding strong and weak Gaussian derivatives coincide almost everywhere. Collecting all derivatives and integrating over \(\mathbb{R}^{d_{\mathcal{X}}}\) gives \[{\| \widehat{f} \|}_{H_{\hat{\varrho}}^{k}(\mathbb{R}^{d_{\mathcal{X}}}; \mathbb{R}^{d_{\mathcal{Y}}})}^2 \leq (k + 1) \int_{\mathbb{R}^{d_{\mathcal{X}}}} \omega({\left\|\boldsymbol{x}\right \|}_2)^2 \,\mathrm{d}\hat{\varrho}(\boldsymbol{x}).\] In conclusion, we can bound \({\| \widehat{f} \|}_{H_{\hat{\varrho}}^{k}}\) independently of the \(\hat{\lambda}_i\) if we can do so for \(\int_{\mathcal{X}} \omega({\left\|\boldsymbol{x}\right \|}_2)^2 \,\mathrm{d}\hat{\varrho}(\boldsymbol{x})\). We discuss three cases for which this is possible.

• No growth. If \(\omega \equiv a\), \(a > 0\), is constant, then trivially, \[\left( \int_{\mathcal{X}} \omega({\left\|\boldsymbol{x}\right \|}_2)^2 \,\mathrm{d}\hat{\varrho}(\boldsymbol{x}) \right)^{1/2} = a.\]

• Affine linear growth. Suppose that \(\omega(t) = a + bt\) with \(a, b > 0\). Suppose that the true and empirical PCA eigenvalues are sufficiently close, say \(\max_{i \in [d_{\mathcal{X}}]} |\lambda_i - \hat{\lambda}_i| \leq \lambda_{d_{\mathcal{X}}}\). Then \[\left( \int_{\mathcal{X}} \omega({\left\|\boldsymbol{x}\right \|}_2)^2 \,\mathrm{d}\hat{\varrho}(\boldsymbol{x}) \right)^{1/2} \leq a + b \sqrt{\sum_{i=1}^{d_{\mathcal{X}}} \hat{\lambda}_i} \leq a + b \sqrt{2 \sum_{i=1}^{d_{\mathcal{X}}} \lambda_i} \leq a + b \sqrt{2},\] where we used the assumption \(\sum_{i=1}^{\dim(\mathcal{X})} \lambda_i = 1\). Note that closeness of the \(\lambda_i\) and \(\hat{\lambda}_i\) can be guaranteed with high probability via Weyl’s inequality 39 if one draws sufficiently many (unlabeled) empirical PCA data samples, see Theorem 5.

• Exponential growth. Suppose that \(\omega(t) = a \exp(b t)\) with \(a, b > 0\). Then, under the same assumptions as in (ii), \[\begin{align} \left( \int_{\mathcal{X}} \omega({\left\|\boldsymbol{x}\right \|}_2)^2 \,\mathrm{d}\hat{\varrho}(\boldsymbol{x}) \right)^{1/2} &\leq a \left( \prod^{d_{\mathcal{X}}}_{i=1} \mathbb{E}_{X \sim \mathcal{N}(0,\hat{\lambda}_i)}[\exp(b X)] \right)^{1/2} = a \exp \left( \frac{b^2}{4} \sum^{d_{\mathcal{X}}}_{i=1} \hat{\lambda}_i \right ) \leq a \exp(b^2 / 2). \end{align}\]

Hence, we can have very fast growth and still maintain control over \({\| \widehat{f} \|}_{H_{\hat{\varrho}}^{k}}\) independently of the \(\hat{\lambda}_i\).

12.2 Examples↩︎

We discuss two examples of important subsets of admissible operators. First, consider the set of all boundedly differentiable \(C^k\)-operators. If \(F \in C_b^k(\mathcal{X}; \mathcal{Y})\), then \(F \in C_{\mu \mathrm{-adm}}^k(\mathcal{X}; \mathcal{Y})\) with constant admissible growth function \(\omega \equiv {\left\|F\right \|}_{C_b^k(\mathcal{X}; \mathcal{Y})}\) and by (i) above, we have \[{\| \widehat{f} \|}_{H_{\hat{\varrho}}^{k}(\mathbb{R}^{d_{\mathcal{X}}}; \mathbb{R}^{d_{\mathcal{Y}}})} \leq \sqrt{k + 1} {\left\|F\right \|}_{C_b^k(\mathcal{X}; \mathcal{Y})}.\] For the second example we introduce the following notation.

If \(F \in C_{\mathrm{Lip}}^k(\mathcal{X}; \mathcal{Y})\), then \(F \in C_{\mu \mathrm{-adm}}^k(\mathcal{X}; \mathcal{Y})\) with admissible growth function \[\omega(t) := [F]_0 + [F]_{C_{\mathrm{Lip}}^k(\mathcal{X}; \mathcal{Y})} t, \quad t \geq 0.\] It follows that \(\widehat{f} \in H_{\hat{\varrho}}^{k}(\mathbb{R}^{d_{\mathcal{X}}}; \mathbb{R}^{d_{\mathcal{Y}}})\), and in the setting of (ii) above, we have \[{\| \widehat{f} \|}_{H_{\hat{\varrho}}^{k}(\mathbb{R}^{d_{\mathcal{X}}}; \mathbb{R}^{d_{\mathcal{Y}}})} \leq \sqrt{k + 1} \left( [F]_0 + \sqrt{2} [F]_{C_{\mathrm{Lip}}^k(\mathcal{X}; \mathcal{Y})} \right) \leq \sqrt{2 (k + 1)} {\left\|F\right \|}_{C_{\mathrm{Lip}}^k(\mathcal{X}; \mathcal{Y})}.\]

References↩︎