A Data-Adaptive Prior for Bayesian Learning
of Kernels in Operators

We aim to learn the kernel $\phi$ in the operator $R_{\phi}:\mathbb{X}\to\mathbb{Y}$ in the following model

$$R_{\phi}[u]+\eta+\xi=f,$$

(1.1)

by fitting the model to the data consisting of $N$ input-output pairs:

$$\mathcal{D}=\{(u^{k},f^{k})\}_{k=1}^{N},\quad(u^{k},f^{k})\in\mathbb{X}\times\mathbb{Y}.$$

(1.2)

Here $\mathbb{X}$ is a Banach space, $\mathbb{Y}$ is a Hilbert space, the measurement noise $\eta$ is a $\mathbb{Y}$-valued white noise in the sense that $\mathbb{E}[\langle{\eta,f}\rangle_{\mathbb{Y}}^{2}]=\sigma_{\eta}^{2}\langle{f,f}\rangle_{\mathbb{Y}}$ for any $f\in\mathbb{Y}$. The term $\xi$ represents unknown model errors such as model misspecification or computational error due to incomplete data, and it may depend on the input data $u$.

The operator $R_{\phi}$ depends non-locally on the kernel $\phi$ in the form

$$R_{\phi}[u](y)=\int_{\Omega}\phi(y-x)g[u](x,y)\mu(dx),\quad\forall y\in\Omega,$$

(1.3)

where $(\Omega,\mu)$ is a measure space that can be either a domain in the Euclidean space with the Lebesgue measure or a discrete set with an atomic measure. For simplicity, we let $\mathbb{Y}=L^{2}(\Omega,\mu)$ throughout this study. Here $g[u]$ is a bounded (nonlinear) functional of $u$ and is assumed to be known. Note that the operator $R_{\phi}$ can be nonlinear in $u$, but it depends linearly on $\phi$. Such operators are widely seen in PDEs, matrix operators, and image processing. Examples include the Toeplitz matrix, integral and nonlocal operators; see Sect.2.1. In these examples, there is often limited prior knowledge about the kernel.

The inverse problem of learning the kernel $\phi$ is often ill-posed, due to the nonlocal dependence of the output data $f$ on the kernel. The ill-posedness, in terms of minimizing the negative log-likelihood of the data,

$$\mathcal{E}(\phi)=\frac{1}{N\sigma_{\eta}^{2}}\sum_{1\leq k\leq N}\|R_{\phi}[u^{k}]-f^{k}\|_{\mathbb{Y}}^{2}=\frac{1}{2\sigma_{\eta}^{2}}\left[\langle{{\mathcal{L}_{\overline{G}}}\phi,\phi}\rangle_{L^{2}_{\rho}}-2\langle{\phi^{\scriptscriptstyle\mathcal{D}},\phi}\rangle_{L^{2}_{\rho}}+C_{N}^{f}\right],$$

(1.4)

appears as the instability of the minimizer ${\mathcal{L}_{\overline{G}}}^{-1}\phi^{\scriptscriptstyle\mathcal{D}}$ when it exists. Here $L^{2}_{\rho}$ is space of square-integrable functions with measure $\rho$, the normal operator ${\mathcal{L}_{\overline{G}}}$ is a trace-class operator, and $\phi^{\scriptscriptstyle\mathcal{D}}$ comes from data; see Sect.2.3 for details. Thus, the ill-posedness is rooted in the unboundedness of ${\mathcal{L}_{\overline{G}}}^{-1}$ and the perturbation of $\phi^{\scriptscriptstyle\mathcal{D}}$.

There are two common strategies to overcome the ill-posedness: a prior in a Bayesian approach and regularization in a variational approach; see, e.g., [53, 23, 3] and a sample of the large literature in Sect.1.3).

This study focuses on selecting a prior for the Bayesian approach. The major challenge is the limited prior information about the kernel and the need to overcome the ill-posedness caused by a data-dependent normal operator.

Due to the lack of prior information about the kernel, one may use the default prior, a non-degenerate prior, e.g., a Gaussian distribution $\mathcal{N}(0,\mathcal{Q}_{0})$ with a nondegenerate covariance operator $\mathcal{Q}_{0}$, with the belief that it is a safe choice to ensure a well-defined posterior.

However, we show that the fixed non-degenerate prior has the risk of leading to a divergent posterior mean in the small noise limit. Specifically, the posterior mean,

$$\mu_{1}=({\mathcal{L}_{\overline{G}}}+\sigma_{\eta}^{2}\mathcal{Q}_{0})^{-1}\phi^{\scriptscriptstyle\mathcal{D}},$$

obtained from the prior $\mathcal{N}(0,\mathcal{Q}_{0})$ and the likelihood that yields (1.4), blows up when, the variance of the noise, $\sigma_{\eta}\to 0$ if $\phi^{\scriptscriptstyle\mathcal{D}}$ contains a perturbation in the null space of the normal operator ${\mathcal{L}_{\overline{G}}}$; see Proposition 3.2. Such perturbation can be caused by any of the four types of errors in data or computation: (a) discretization error, (b) model error, (c) partial observations, and (d) wrong noise assumption. Thus, a prior adaptive to ${\mathcal{L}_{\overline{G}}}$ and $\phi^{\scriptscriptstyle\mathcal{D}}$ is needed to remove the risk.

We propose a data-adaptive RKHS prior that ensures a stable posterior mean in the small noise regime. It is a Gaussian distribution $\mathcal{N}(0,\lambda_{*}^{-1}{\mathcal{L}_{\overline{G}}})$ with the parameter $\lambda_{*}$ selected adaptive to $\phi^{\scriptscriptstyle\mathcal{D}}$. We prove in Theorem 4.2 that it leads to a stable posterior whose mean

$$\mu_{1}^{\scriptscriptstyle\mathcal{D}}=({\mathcal{L}_{\overline{G}}}^{2}+\sigma_{\eta}^{2}\lambda_{*}I_{\mathrm{Null({\mathcal{L}_{\overline{G}}})}^{\perp}})^{-1}{\mathcal{L}_{\overline{G}}}\phi^{\scriptscriptstyle\mathcal{D}},$$

always has a small noise limit, and the small noise limit converges to the identifiable parts of the true kernel. Furthermore, we show that our prior outperforms the non-degenerate prior in producing a more accurate posterior mean and smaller posterior uncertainty in terms of the trace of the posterior covariance; see Sect.4.2. The prior is called an RKHS prior because its Cameron-Martin space ${\mathcal{L}_{\overline{G}}}^{1/2}L^{2}_{\rho}$ is the RKHS with a reproducing kernel ${\overline{G}}$ determined by the operator $R_{\phi}$ and the data. Importantly, the closure of this RKHS is the space in which the components of the true kernel can be identified from data.

We also study the computational practice of the data-adaptive prior and demonstrate it on the Toeplitz matrices and integral operators. We select the hyper-parameter by the L-curve method in [24]. Numerical tests show that while a fixed non-degenerate prior leads to divergent posterior means, the data-adaptive prior always attains stable posterior means with small noise limits; see Sect.5.

The outline of this study is as follows. Sect.2 presents the mathematical setup of this study, and shows the ill-posedness of the inverse problem through the variational approach. We introduce in Sect.3 the Bayesian approach and show the issue of a fixed non-degenerate prior. To solve the issue, we introduce a data-adaptive prior in Sect.4, and analyze its advantages. Sect.5 discusses the computational practice and demonstrates the advantage of the data-adaptive prior in numerical tests on Toeplitz matrices and integral operators. Finally, Sect.6 concludes our findings and provides future research directions.

We first present a few examples of learning kernels in operators. Note that in these examples, there is little prior information about the kernel.

To identify the kernel, the variational approach seeks a maximal likelihood estimator:

$$\widehat{\phi}=\operatorname*{arg\,min}_{\phi\in\mathcal{H}}\mathcal{E}(\phi),\quad\text{ where }\mathcal{E}(\phi)=\frac{1}{N\sigma_{\eta}^{2}}\sum_{1\leq k\leq N}\|R_{\phi}[u^{k}]-f^{k}\|_{\mathbb{Y}}^{2},$$

(2.3)

over a hypothesis space $\mathcal{H}$, where the loss function $\mathcal{E}(\phi)$ is the negative log-likelihood of the data (1.2) under the assumption that the noise $\eta$ in (1.1) is white.

The first step is to find a proper function space for $\phi$, in which one can further specify the hypothesis space $\mathcal{H}$. Clearly, given a data set, we can only identify the kernel where the data provides information. Examples 2.1– 2.4 show that the support of the kernel $\phi$ is yet to be extracted from data.

We set function space for learning $\phi$ to be $L^{2}(\mathcal{S},\rho)$ with $\mathcal{S}=\{x-y:x,y\in\Omega\}$, where $\rho$ is an empirical probability measure quantifying the exploration of data to the kernel:

$\displaystyle\rho(dr)=\frac{1}{ZN}\sum_{1\leq k\leq N}\int_{\Omega}\int_{\Omega}\delta(y-x-r)\left|g[u^{k}](x,y)\right|\mu(dx)\mu(dy),\quad r\in\mathcal{S}.$

(2.4)

Here $Z$ is the normalizing constant and $\delta$ is the Dirac delta function. We call $\rho$ an exploration measure. Its support is the region in which the learning process ought to work, and outside of which, we have limited information from the data to learn the function $\phi$. Thus, we can restrict $\mathcal{S}$ to be the support of $\rho$, and we denote $L^{2}_{\rho}:=L^{2}(\mathcal{S},\rho)$ for short. In other words, the exploration measure is a generalization of the measure $\rho_{X}$ in nonparametric regression that estimates $Y=f(X)$ from data $\{(x_{i},y_{i})\}$, where $\rho_{X}$ is the distribution of $X$ and the data are samples from the joint distribution of $(X,Y)$.

The next step is to find the minimizer of the loss function. Note that $\mathcal{E}(\phi)$ is quadratic in $\phi$ since the operator $R_{\phi}$ depends linearly on $\phi$. We get the minimizer by solving the zero of the (Fréchet) derivative of the loss function. To compute its derivative, we first introduce a bilinear form $\langle\langle{\cdot,\cdot}\rangle\rangle$: $\forall\phi,\psi\in L^{2}_{\rho}$,

	$\displaystyle\langle\langle{\phi,\psi}\rangle\rangle=$	$\displaystyle\frac{1}{N}\sum_{1\leq k\leq N}\langle{R_{\phi}[u^{k}],R_{\psi}[u^{k}]}\rangle_{\mathbb{Y}},$		(2.5)
	$\displaystyle=$	$\displaystyle\frac{1}{N}\sum_{1\leq k\leq N}\int\left[\int\int\phi(y-x)\psi(y-z)g[u^{k}](x,y)g[u^{k}](z,y)\mu(dx)\mu(dz)\right]\mu(dy)$
	$\displaystyle=$	$\displaystyle\int_{\mathcal{S}}\int_{\mathcal{S}}\phi(r)\psi(s){\overline{G}}(r,s)\rho(dr)\rho(ds),$

where the integral kernel ${\overline{G}}$ given by, for $r,s\in\mathrm{supp}(\rho)$,

$${\overline{G}}(r,s)=\frac{G(r,s)}{\rho(r)\rho(s)}\quad\text{ with }G(r,s)=\frac{1}{N}\sum_{1\leq k\leq N}\int g[u^{k}](x,r+x)g[u^{k}](x,s+x)\mu(dx),$$

(2.6)

in which by an abuse of notation, we also use $\rho(r)$ to denote either the probability of $r$ when $\rho$ defined in (2.4) is discrete or the probability density of $\rho$ when the density exists.

By definition, the bivariate function ${\overline{G}}$ is symmetric and positive semi-definite in the sense that $\sum_{i,j=1}^{n}c_{i}c_{j}G(r_{i},r_{j})\geq 0$ for any $\{c_{i}\}_{i=1}^{n}\subset\mathbb{R}$ and $\{r_{i}\}_{i=1}^{n}\subset\mathcal{S}$. In the following, we assume that the data is continuous and bounded so that ${\overline{G}}$ defines a self-adjoint compact operator which is fundamental for studying identifiability. This assumption holds under mild regularity conditions on the data $\{u^{k}\}$ and the operator $R_{\phi}$.

Under Assumption 2.5, the integral operator ${\mathcal{L}_{\overline{G}}}:L^{2}_{\rho}\to L^{2}_{\rho}$

$${\mathcal{L}_{\overline{G}}}\phi(r)=\int_{\mathcal{S}}\phi(s)\overline{G}(r,s)\rho(s)ds,$$

(2.7)

is a positive semi-definite trace-class operator (see Lemma A.1). Hereafter we denote $\{\lambda_{i},\psi_{i}\}$ the eigen-pairs of ${\mathcal{L}_{\overline{G}}}$ with the eigenvalues in descending order, and assume that the eigenfunctions are orthonormal, hence they provide an orthonormal basis of $L^{2}_{\rho}$. Furthermore, for any $\phi,\psi\in L^{2}_{\rho}$, the bilinear form in (2.5) can be written as

$$\langle\langle{\phi,\psi}\rangle\rangle=\langle{{\mathcal{L}_{\overline{G}}}\phi,\psi}\rangle_{L^{2}_{\rho}}=\langle{\phi,{\mathcal{L}_{\overline{G}}}\psi}\rangle_{L^{2}_{\rho}},$$

(2.8)

and we can write the loss functional in (2.3) as

	$\displaystyle\mathcal{E}(\phi)$	$\displaystyle=\langle\langle{\phi,\phi}\rangle\rangle-2\frac{1}{N}\sum_{1\leq k\leq N}\langle{R_{\phi}[u^{k}],f^{k}}\rangle_{\mathbb{Y}}+\frac{1}{N}\sum_{1\leq k\leq N}\|f^{k}\|_{\mathbb{Y}}^{2}$		(2.9)
		$\displaystyle=\langle{{\mathcal{L}_{\overline{G}}}\phi,\phi}\rangle_{L^{2}_{\rho}}-2\langle{\phi^{\scriptscriptstyle\mathcal{D}},\phi}\rangle_{L^{2}_{\rho}}+C_{N}^{f},$

where $\phi^{\scriptscriptstyle\mathcal{D}}\in L^{2}_{\rho}$ is the Riesz representation of the bounded linear functional:

$$\langle{\phi^{\scriptscriptstyle\mathcal{D}},\psi}\rangle_{L^{2}_{\rho}}=\frac{1}{N}\sum_{1\leq k\leq N}\langle{R_{\psi}[u^{k}],f^{k}}\rangle_{\mathbb{Y}},\quad\forall\psi\in L^{2}_{\rho}.$$

(2.10)

Then, by solving the zero of $\nabla\mathcal{E}(\phi)=2({\mathcal{L}_{\overline{G}}}\phi-\phi^{\scriptscriptstyle\mathcal{D}})$, one may obtain the minimizer ${\mathcal{L}_{\overline{G}}}^{-1}\phi^{\scriptscriptstyle\mathcal{D}}$, provided that it is well-defined. However, it is often ill-defined, e.g., when $\phi^{\scriptscriptstyle\mathcal{D}}$ is not in ${\mathcal{L}_{\overline{G}}}(L^{2}_{\rho})$ and ${\mathcal{L}_{\overline{G}}}$ is compact infinite-rank or rank-deficient. Thus, it is important to examine the inversion and the uniqueness of the minimizer, for which we introduce the function space of identifiability in the next section.

We define a function space of identifiability (FSOI) when one learns the kernel by minimizing the loss function.

The next theorem characterizes the FSOI. Its proof is deferred to Appendix A.1.

Theorem 2.7 enables us to analyze the ill-posedness of the inverse problems through the operator ${\mathcal{L}_{\overline{G}}}$ and $\phi^{\scriptscriptstyle\mathcal{D}}$. When $\phi^{\scriptscriptstyle\mathcal{D}}\in{\mathcal{L}_{\overline{G}}}(L^{2}_{\rho})$, the inverse problem has a unique solution in the FSOI $H$, even when it is underdetermined in $L^{2}_{\rho}$ due to $H$ being a proper subspace, which happens when the compact operator ${\mathcal{L}_{\overline{G}}}$ has a zero eigenvalue. However, when $\phi^{\scriptscriptstyle\mathcal{D}}\notin{\mathcal{L}_{\overline{G}}}(L^{2}_{\rho})$, the inverse problem $\nabla\mathcal{E}=0$ has no solution in $L^{2}_{\rho}$ because ${\mathcal{L}_{\overline{G}}}^{-1}\phi^{\scriptscriptstyle\mathcal{D}}$ is undefined. According to (2.11), this happens in one or more of the following scenarios:

• •

  when the model error leads to $\epsilon^{\xi}\notin{\mathcal{L}_{\overline{G}}}(L^{2}_{\rho})$.

• •

  when the observation noise leads to $\epsilon^{\eta}\notin{\mathcal{L}_{\overline{G}}}(L^{2}_{\rho})$. In particular, since $\epsilon^{\eta}$ is Gaussian $\mathcal{N}(0,{\mathcal{L}_{\overline{G}}})$, it has the Karhunen–Loève expansion $\epsilon^{\eta}=\sum_{i}\lambda_{i}^{1/2}\epsilon^{\eta}_{i}\psi_{i}$ with $\epsilon^{\eta}_{i}$ being i.i.d. $\mathcal{N}(0,1)$. Then, ${\mathcal{L}_{\overline{G}}}^{-1}\epsilon^{\eta}=\sum_{i}\lambda_{i}^{-1/2}\epsilon^{\eta}_{i}\psi_{i}$, which diverges a.s. if $\sum_{i:\lambda_{i}>0}\lambda_{i}^{-1}$ diverges. Thus, we have $\epsilon^{\eta}\notin{\mathcal{L}_{\overline{G}}}(L^{2}_{\rho})$ a.s. when $\sum_{i:\lambda_{i}>0}\lambda_{i}^{-1}$ diverges.

Additionally, $\phi^{\scriptscriptstyle\mathcal{D}}$ only encodes information of $\phi_{true}^{H}$, and it provides no information about $\phi_{true}^{\perp}$, where $\phi_{true}^{H}$ and $\phi_{true}^{\perp}$ form an orthogonal decomposition $\phi_{true}=\phi_{true}^{H}+\phi_{true}^{\perp}\in H\oplus H^{\perp}$. In other words, the data provides no information to recover $\phi_{true}^{\perp}$.

As a result, it is important to avoid absorbing the errors outside of the FSOI when using a Bayesian approach or a regularization method to mitigate the ill-posedness.

The RKHS with ${\overline{G}}$ as a reproducing kernel is a dense subspace of the FSOI. Hence, when using it as a prior, which we will detail in later sections, one can filter out the error outside the FSOI and ensure that the learning takes place inside the FSOI.

The next lemma is an operator characterization of this RKHS (see, e.g., [9, Section 4.4]). Its proof can be found in [38, Theorem 3.3].

The above RKHS has been used for regularization in [37, 38]. The regularizer, named DARTR, regularizes the loss by the norm of this RKHS,

$$\mathcal{E}_{\lambda}(\phi)=\mathcal{E}(\phi)+\lambda\|\phi\|_{H_{G}}^{2}=\langle{({\mathcal{L}_{\overline{G}}}+\lambda{\mathcal{L}_{\overline{G}}}^{-1})\phi,\phi}\rangle_{L^{2}_{\rho}}-2\langle{\phi^{\scriptscriptstyle\mathcal{D}},\phi}\rangle_{L^{2}_{\rho}}+C_{N}^{f}.$$

(2.14)

Selecting the optimal hyper-parameter $\lambda_{*}$ by the L-curve method, it leads to the estimator

$$\widehat{\phi}_{H_{G}}=({\mathcal{L}_{\overline{G}}}^{2}+\lambda_{*}I_{H})^{-1}{\mathcal{L}_{\overline{G}}}\phi^{\scriptscriptstyle\mathcal{D}},$$

(2.15)

where $I_{H}$ is the identity operator on $H$. Since the RKHS is a subset of the FSOI with more regular elements, DARTR ensures that the estimator is in the FSOI and is regularized.

We summarize the key terminologies and notations in this section in Table 1.

In this study, we focus on Gaussian priors, so the posterior is also a Gaussian measure in combination with a Gaussian likelihood. Also, without loss of generality, we assume that the prior is centered. Recall that the function space of learning is $L^{2}_{\rho}$ defined in (2.4). For illustration, we first specify the prior and posterior when the space $L^{2}_{\rho}$ is finite-dimensional, then discuss them in the infinite-dimensional case. We follow the notations in [53] and [14].

The prior distribution plays a crucial role in Bayesian inverse problems. To make the ill-posed inverse problem well-defined, it is often set to be a non-degenerate measure (i.e., its covariance operator $\mathcal{Q}_{0}$ has no zero eigenvalues). It is fixed in many cases and not adaptive to data. Such a non-degenerate prior works well for an inverse problem whose function space of identifiability does not change with data. However, in the learning of kernels in operators, a non-degenerate prior has a risk of leading to a catastrophic error: the posterior mean may diverge in the small observation noise limit, as we show in the Proposition 3.2.

Assumptions (A1-A2) are common in practice. Assumption (A1) holds because the operator ${\mathcal{L}_{\overline{G}}}$ is finite rank when the data is discrete, and it is not full rank for under-determined problems. It is natural to assume the prior has a full rank covariance $\mathcal{Q}_{0}$ as in (A2). We assume that $\mathcal{Q}_{0}$ commutes with ${\mathcal{L}_{\overline{G}}}$ for simplicity and one can extend it to the general case as in the proof of [12, Theorem 2.25, Feldman–Hajek theorem]. Assumption (A3), which requires $\phi^{\scriptscriptstyle\mathcal{D}}$ to be outside the range of ${\mathcal{L}_{\overline{G}}}$, holds when the regression vector ${\overline{b}}$ is outside the range of the regression matrix ${\overline{A}}$ in (5.3b), see Sect.5.2–5.3 for more discussions.

On the other hand, one may adjust the prior covariance by the standard deviation of the noise in the hope of removing the risk of divergence. However, the next proposition shows that such a noise-adaptive non-degenerate prior will have a biased small noise limit.