Scalable iterative data-adaptive RKHS regularization

Numerous regularization methods have been developed, and the literature on this topic is too vast to be surveyed here; we refer to [15, 11, 12] and references therein for an overview. In the following, we compare iDARR with the most closely related works.

The linear equation eq. 1.1 can arise from discrete or continuous inverse problems. In either case, we can present the inverse problem using the prototype of the Fredholm integral equation of the first kind. We consider only the 1D case for simplicity, and the extension to high-dimensional cases is straightforward. Specifically, let $\mathcal{S},\mathcal{T}\subset\mathbb{R}$ be two compact sets. We aim to recover the function $\phi:\mathcal{S}\to\mathbb{R}$ in the Fredholm equation

$\displaystyle y(t)=\int_{\mathcal{S}}K(t,s)\phi(s)\nu(ds)+\sigma\dot{W}(t)=:L_{K}\phi(t)+\sigma\dot{W}(t),\quad\forall t\in\mathcal{T}$

(2.1)

from discrete noisy data

$$b=(y(t_{1}),\ldots,y(t_{m}))^{\top}\in\mathbb{R}^{m},$$

where we assume the observation index $\mathcal{T}=\{t_{j}\}_{j=1}^{m}$ to be $0=t_{0}<t_{1}<\cdots<t_{m}$ for simplicity. Here the measurement noise $\sigma\dot{W}(t)$ is the white noise scaled by $\sigma$; that is, the noise at $t_{j}$ has a Gaussian distribution $\mathcal{N}(0,\sigma^{2}(t_{j+1}-t_{j}))$ for each $j$. Such noise is integrable when the observation mesh refines, i.e., $\max_{j}(t_{j+1}-t_{j})$ vanishes.

Here the finite measure $\nu$ can be either the Lebesgue measure with $\mathcal{S}$ being an interval or an atom measure with $\mathcal{S}$ having finitely many elements. Correspondingly, the Fredholm integral equation eq. 2.1 is either a continuous or a discrete model.

In either case, the goal is to solve for the function $\phi:\mathcal{S}\to\mathbb{R}$ in eq. 2.1. When seeking a solution in the form of $\phi=\sum_{i=1}^{n}x_{i}\phi_{i}$, where $\{\phi_{i}\}_{i=1}^{n}$ is a pre-selected set of basis functions, we obtain the linear equation eq. 1.1 with $x=(x_{1},\ldots,x_{n})^{\top}\in\mathbb{R}^{n}$ and the matrix $A$ with entries

$$A(j,i)=\int_{\mathcal{S}}K(t_{j},s)\phi_{i}(s)\nu(ds)=L_{K}\phi_{i}(t_{j}),\quad 1\leq j\leq m,1\leq i\leq n.$$

(2.2)

In particular, when $\{\phi_{i}\}$ are piece-wise constants, we obtain $A$ as follows.

• •

  Discrete model. Let $\nu$ be an atom measure on $\mathcal{S}=\{s_{i}\}_{i=1}^{n}$, a set with $n$ elements. Suppose that the basis functions are $\phi_{i}(s)=\mathbf{1}_{\{s_{i}\}}(s)$. Then, $\phi=x$ and the matrix $A$ has entries $A(j,i)=K(t_{j},s_{i})\nu(s_{i})$.

• •

  Continuous model. Let $\nu$ be the Lebesgue measure on $\mathcal{S}=[0,1]$, and $\phi_{i}(s)=\mathbf{1}_{[s_{i-1},s_{i}]}(s)$ be piecewise constant functions on a partition of $\mathcal{S}$ with $0=s_{0}<s_{1}<\ldots<s_{n}=1$. Then, $\phi=\sum_{i=1}^{n}x_{i}\phi_{i}$ and the matrix $A$ has entries $A(j,i)=K(t_{j},s_{i})(s_{i}-s_{i-1})$.

The default function spaces for $\phi$ and $y$ above are $L^{2}_{\nu}(\mathcal{S})$ and $L^{2}_{\mu}(\mathcal{T})$. The loss function $\mathcal{E}(x)=\|Ax-b\|_{2}^{2}$ over $L^{2}_{\nu}(\mathcal{S})$ becomes

$\displaystyle\mathcal{E}(\phi):=\left\|L_{K}\phi-y\right\|^{2}_{{L^{2}_{\mu}(\mathcal{T})}}$

$\displaystyle=\langle{L_{K}\phi,L_{K}\phi}\rangle_{{L^{2}_{\mu}(\mathcal{T})}}-2\langle{L_{K}\phi,y}\rangle_{{L^{2}_{\mu}(\mathcal{T})}}+\left\|y\right\|^{2}_{{L^{2}_{\mu}(\mathcal{T})}},$

(2.3)

Eq.eq. 2.1 is a prototype of ill-posed inverse problems, dating back from Hadamard [14], and it remains a testbed for new regularization methods [36, 28, 15, 23].

The $L^{2}_{\nu}(\mathcal{S})$ norm is often a default choice for regularization. However, it has a major drawback: it does not take into account the operator $L_{K}$, particularly when $L_{K}$ has zero eigenvalues, and it leads to unstable solutions that may blow up in the small noise limit [22]. To avoid such instability, particularly for iterative methods, we introduce a weighted function space and an RKHS that are adaptive to both the data and the model in the next sections.

We first introduce a weighted function space $L^{2}_{\rho}(\mathcal{S})$, where the measure $\rho$ is defined as

$$\frac{d\rho}{d\nu}(s):=\frac{1}{Z}\int_{\mathcal{T}}|K(t,s)|\mu(dt),\,\forall s\in\mathcal{S},$$

(2.4)

where $Z=\int_{\mathcal{S}}\int_{\mathcal{T}}|K(t,s)|\mu(dt)\nu(ds)$ is a normalizing constant. This measure quantifies the exploration of data to the unknown function through the integral kernel $K$ at the output set $\mathcal{T}$, i.e., $\{K(t_{j},\cdot)\}_{t_{j}\in\mathcal{T}}$, hence it is referred to as an exploration measure. In particular, when (1.1) is a discrete model, the exploration measure is the normalized column sum of the absolute values of the matrix $A$.

The major advantage of the space $L^{2}_{\rho}(\mathcal{S})$ over the original space $L^{2}_{\nu}(\mathcal{S})$ is that it is adaptive to the specific setting of the inverse problem. In particular, this weighted space takes into account the structure of the integral kernel and the data points in $\mathcal{T}$. Thus, while the following introduction of RKHS can be carried out in both $L^{2}_{\rho}(\mathcal{S})$ and $L^{2}_{\nu}(\mathcal{S})$, we will focus only on $L^{2}_{\rho}(\mathcal{S})$.

Next, we consider the variational inverse problem over $L^{2}_{\rho}$, and the goal is to find a minimizer of the loss function eq. 2.3 in $L^{2}_{\rho}$. Since the loss function is quadratic, its Hessian is a symmetric positive linear operator, and it has a unique minimizer in the linear subspace where the Hessian is strictly positive. We assign a name to this subspace in the next definition.

The next lemma specifies the FSOI for the loss function in eq. 2.3 (see [26] for its proof).

Lemma 2.2 reveals the cause of the ill-posedness, and provides insights on regularization:

• •

  The variational inverse problem is well-defined only in the FSOI $H$, which can be a proper subset of $L^{2}_{\rho}$. Its ill-posedness in $H$ depends on the smallest eigenvalue of the operator ${\mathcal{L}_{\overline{G}}}$ and the error in $\phi^{D}$.

• •

  When the data is noiseless, the loss function can only identify the $H$-projection of the true input function. When data is noisy, its minimizer ${\mathcal{L}_{\overline{G}}}^{-1}\phi^{D}$ is ill-defined in $L^{2}_{\rho}$ when $\phi^{D}\notin{\mathcal{L}_{\overline{G}}}(H)$.

As a result, when regularizing the ill-posed problem, an important task is to ensure the search takes place in the FSOI and to remove the noise-contaminated components making $\phi^{D}\notin{\mathcal{L}_{\overline{G}}}(H)$.

Our data-adaptive RKHS is the RKHS with ${\overline{G}}$ in eq. 2.5 as a reproducing kernel. Hence, it is adaptive to the integral kernel $K$ and the data in the model. When its norm is used for regularization, it ensures that the search takes place in the FSOI because its $L^{2}_{\rho}$ closure is the FSOI; also, it penalizes the components in $\phi^{D}$ corresponding to the small singular values.

The next lemma characterizes this RKHS, and we refer to [26] for its proof.

The next theorem shows the computation of the RKHS norm for the problem eq. 1.1 when it is written in the form eq. 2.1-eq. 2.2. A key component is solving the eigenvalues of ${\mathcal{L}_{\overline{G}}}$ through a generalized eigenvalue problem.

The orthonormality of $\{\psi_{i}\}$ implies that

$$I_{n}=\big{(}\langle{\psi_{k},\psi_{l}}\rangle_{L^{2}_{\rho}}\big{)}_{1\leq k,l\leq n}=\big{(}\langle{\sum_{i=1}^{n}V_{ik}\phi_{i},\sum_{j=1}^{n}V_{jl}\phi_{j}}\rangle_{L^{2}_{\rho}}\big{)}_{1\leq k,l\leq n}=V^{\top}BV.$$

Note that $[A^{\top}A](j,i)=\langle{L_{K}\phi_{i},L_{K}\phi_{j}}\rangle_{L^{2}_{\mu}(\mathcal{T})}=\langle{\phi_{j},{\mathcal{L}_{\overline{G}}}\phi_{i}}\rangle_{L^{2}_{\rho}}$ for $1\leq i,j\leq n$. Then,

$$\langle{\phi_{j},{\mathcal{L}_{\overline{G}}}\psi_{k}}\rangle_{L^{2}_{\rho}}=\sum_{i=1}^{n}[A^{\top}A](j,i)V_{ik}.$$

Meanwhile, the eigen-equation ${\mathcal{L}_{\overline{G}}}\psi_{k}=\lambda_{k}\psi_{k}$ implies that for each $\phi_{j}$,

$$\langle{\phi_{j},{\mathcal{L}_{\overline{G}}}\psi_{k}}\rangle_{L^{2}_{\rho}}=\lambda_{k}\langle{\phi_{j},\psi_{k}}\rangle_{L^{2}_{\rho}}=\lambda_{k}\langle{\phi_{j},\sum_{i=1}^{n}V_{ik}\phi_{i}}\rangle_{L^{2}_{\rho}}=\lambda_{k}\sum_{i=1}^{n}B_{ji}V_{ik},$$

i.e., $\big{(}\langle{\phi_{j},{\mathcal{L}_{\overline{G}}}\psi_{k}}\rangle_{L^{2}_{\rho}}\big{)}=BV\Lambda$. Hence, these two equations imply that $A^{\top}AV=BV\Lambda$.

Next, to compute the norm of $\phi=\sum_{i=1}^{n}x_{i}\phi_{i}\in H_{G}$, we write it as $\phi=x^{\top}\boldsymbol{\Phi}=x^{\top}V^{-1}\boldsymbol{\Psi}$. Then, its norm is

$$\|\phi\|_{H_{G}}^{2}=\sum_{k=1}^{n}\lambda_{k}^{-1}\big{(}x^{\top}V^{-1}\big{)}_{k}^{2}=x^{\top}V^{-1}\Lambda^{\dagger}V^{-\top}x=x^{\top}(V\Lambda V^{\top})^{\dagger}x.$$

Thus, $C_{rkhs}=(V\Lambda V^{\top})^{\dagger}=B(A^{\top}A)^{\dagger}B$ and $C_{rkhs}^{\dagger}=V\Lambda V^{\top}=B^{-1}(A^{\top}A)B^{-1}$.   

In particular, when eq. 1.1 is either a discrete model or a discretization of eq. 2.1 based on Riemann sum approximation of the integral, the exploration measure $\rho$ is the normalized column sum of the absolute values of the matrix $A$, and $B=\mathrm{diag}(\rho)$. See Section 5.1 for details.

We review DARTR, a data-adaptive RKHS Tikhonov regularization algorithm introduced in [26].

Specifically, it solves the problem eq. 1.1 with regularization:

$$(\widehat{x}_{\lambda_{*}},\lambda_{*})=\mathop{\mathrm{argmin}}_{x\in\mathbb{R}^{n},\lambda\in\mathbb{R}^{+}}\mathcal{E}_{\lambda}(x),\quad\text{ where }\mathcal{E}_{\lambda}(x):=\|Ax-b\|^{2}+\lambda\|x\|_{C_{rksh}}^{2},$$

where the norm $\|\cdot\|_{C_{rksh}}$ is the DA-RKHS norm introduced in Theorem 2.4. A direct solution minimizing $\mathcal{E}_{\lambda}(x)$ is to solve $(A^{\top}A+\lambda C_{rkhs})x_{\lambda}=A^{\top}b$. However, the computation of $C_{rkhs}$ requires a pseudo-inverse that may cause numerical instability.

DARTR introduces a transformation matrix $C_{*}:=V\Lambda^{1/2}$ to avoid using the pseudo-inverse. Note that $C_{*}^{\top}C_{rkhs}C_{*}=\begin{pmatrix}I_{r}&0\\ 0&0\end{pmatrix}:=\mathbf{I}_{r}$, where $I_{r}$ is the identity matrix with rank $r$, the number of positive eigenvalues in $\Lambda$. Then, the linear equation $(A^{\top}A+\lambda C_{rkhs})x_{\lambda}=A^{\top}b$ is equivalent to

$$(C_{*}A^{\top}AC_{*}+\lambda\mathbf{I}_{r})\widetilde{x}_{\lambda}=C_{*}A^{\top}b$$

(2.10)

with $\widetilde{x}_{\lambda}=C_{*}^{-1}x_{\lambda}$. Thus, DARTR computes $\widetilde{x}_{\lambda}$ in the above equation by least squares with minimal norm, and returns $x_{\lambda}=C_{*}\widetilde{x}_{\lambda}$.

DARTR is a direct method based on matrix decomposition, and it takes $O(mn^{2}+n^{3})$ flops. Hence, it is computationally infeasible when $n$ is large. The iterative method in the next section implements the RKHS regularization in a scalable fashion.

Our regularization method is based on subspace projection in the DA-RKHS. It iteratively constructs a sequence of linear subspaces $\mathcal{S}_{k}$ of the DA-RKHS $(\mathbb{R}^{n},\langle\cdot,\cdot\rangle_{C_{rkhs}})$, and recursively solves projected problems

$$x_{k}=\mathop{\mathrm{argmin}}_{x\in\mathcal{X}_{k}}\|x\|_{C_{rkhs}},\ \ \mathcal{X}_{k}=\{x:\min_{x\in\mathcal{S}_{k}}\|Ax-b\|_{2}\}.$$

(3.1)

This process yields a sequence of solutions $\{x_{k}\}$, each emerging from its corresponding subspace. We ensure the uniqueness of the solution within each iteration, as detailed in Theorem 4.7. The iteration proceeds until it meets an early stopping criterion, designed to exclude excessive noisy components and thereby achieve effective regularization. The spaces $\mathcal{S}_{k}$ are called the solution subspaces, and the iteration number $k$ plays the role of the regularization parameter.

Algorithm 1 outlines this procedure, which is a recursion of the following three parts.

1. P1

   Construct the solution subspaces. We introduce a new generalized Golub-Kahan bidiagonalization (gGKB) to construct the solution subspaces in the DA-RKHS iteratively. The procedure is presented in Algorithm 2.

2. P2

   Recursively update the solution to the projected problem. We solve the least squares problem in the solution subspaces in eq. 3.1 efficiently by a new LSQR-type algorithm, Algorithm 3, updating $\|x_{k}\|_{C_{rkhs}}$ and the residual norm $\|Ax-b\|_{2}$ without even computing the residual.

3. P3

   Regularize by early stopping. We select the optimal $k$ by either the discrepancy principle (DP) when we have an accurate estimate of $\|\mathbf{w}\|_{2}$ or the L-curve criterion otherwise.

In the next subsection, we present details for these key parts. Then, we analyze the computational complexity of the algorithm.

This subsection presents detailed derivations for the three parts P1-P3 in Algorithm 1.

Consider first the case where $C_{rkhs}$ is positive definite. In this scenario, $A$ has full column rank, and the $C_{rkhs}$-inner product Hilbert space $(\mathbb{R}^{n},\langle\cdot,\cdot\rangle_{C_{rkhs}})$ is a discrete representation of the RKHS $H_{G}$ with the given basis functions. Note that the true solution is mapped to the noisy $b$ by $A:(\mathbb{R}^{n},\langle\cdot,\cdot\rangle_{C_{rkhs}})\to(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{2})$. Let $A^{*}:(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{2})\to(\mathbb{R}^{n},\langle\cdot,\cdot\rangle_{C_{rkhs}})$ be the adjoint of $A$, i.e. $\langle Ax,b\rangle_{2}=\langle x,A^{*}b\rangle_{C_{rkhs}}$ for any $x\in\mathbb{R}^{n}$ and $b\in\mathbb{R}^{m}$. By definition, the matrix-form expression of $A^{*}$ is

$$A^{*}=C_{rkhs}^{-1}A^{\top},$$

(3.2)

since $\langle Ax,b\rangle_{2}=x^{\top}A^{\top}b$ and $\langle x,A^{*}b\rangle_{C_{rkhs}}=x^{\top}C_{rkhs}A^{*}b$ for any $x$ and $b$.

The Golub-Kahan bidiagonalization (GKB) process recursively constructs orthonormal bases for these two Hilbert spaces starting with the vector $b$ as follows:

	$\displaystyle\beta_{1}u_{1}=b,\ \ \alpha_{1}z_{1}=A^{*}u_{1},$		(3.3a)
	$\displaystyle\beta_{i+1}u_{i+1}=Az_{i}-\alpha_{i}u_{i},$		(3.3b)
	$\displaystyle\alpha_{i+1}z_{i+1}=A^{*}u_{i+1}-\beta_{i+1}z_{i},$		(3.3c)

where $u_{i}\in(\mathbb{R}^{m},\langle\cdot,\cdot\rangle_{2})$, $z_{i}\in(\mathbb{R}^{n},\langle\cdot,\cdot\rangle_{C_{rkhs}})$ with $z_{0}=\boldsymbol{0}$, and $\alpha_{i}$, $\beta_{i}$ are normalizing factor such that $\|u_{i}\|_{2}=\|z_{i}\|_{C_{rkhs}}=1$. The iteration starts with $u_{1}=b/\beta_{1}$ with $\beta_{1}=\|b\|_{2}$. Using $A^{*}$ in eq. 3.2, we write eq. 3.3c as

$\displaystyle\alpha_{i+1}z_{i+1}=C_{rkhs}^{-1}A^{\top}u_{i+1}-\beta_{i+1}z_{i}$

(3.4)

with $\alpha_{i+1}=\|C_{rkhs}^{-1}A^{\top}u_{i+1}-\beta_{i+1}z_{i}\|_{C_{rkhs}}$. To compute $\alpha_{i+1}$ without explicitly computing $C_{rkhs}$, define $\bar{z}_{i}=C_{rkhs}z_{i}$. Then we have

$\displaystyle\alpha_{i+1}\bar{z}_{i+1}=A^{\top}u_{i+1}-\beta_{i+1}\bar{z}_{i},$

(3.5)

where $\bar{z}_{0}:=\boldsymbol{0}$. Let $p=A^{\top}u_{i+1}-\beta_{i+1}\bar{z}_{i}$. Then we obtain $\alpha_{i+1}=\|C_{rkhs}^{-1}p\|_{C_{rkhs}}=(p^{\top}C_{rkhs}^{-1}p)^{1/2}$, which uses $C_{rkhs}^{-1}=B^{-1}A^{\top}AB^{-1}$ without computing $C_{rkhs}$.

Next, consider that $C_{rkhs}$ is positive semidefinite. The iterative process remains the same with $C_{rkhs}^{-1}$ replaced by the pseudo-inverse $C_{rkhs}^{{\dagger}}$, because $C_{rkhs}^{{\dagger}}=B^{-1}A^{\top}AB^{-1}$ has the same form as $C_{rkhs}^{-1}$ for the non-singular case. Specifically, the recursive relation eq. 3.4 becomes

$\displaystyle\alpha_{i+1}z_{i+1}=C_{rkhs}^{{\dagger}}A^{\top}u_{i+1}-\beta_{i+1}z_{i}.$

(3.6)

To compute $\alpha_{i+1}$, we use the property that $z_{i}\in\mathcal{R}(C_{rkhs})$, which will be proved in Proposition 4.3. Note that $C_{rkhs}^{{\dagger}}C_{rkhs}=P_{\mathcal{N}(C_{rkhs})^{\perp}}=P_{\mathcal{R}(C_{rkhs})}$ since $C_{rkhs}$ is symmetric, where $P_{\mathcal{X}}$ is the projection operator onto subspace $\mathcal{X}$. It follows that $C_{rkhs}^{{\dagger}}C_{rkhs}z_{i}=z_{i}$. Therefore, eq. 3.6 becomes

$$\alpha_{i+1}z_{i+1}=C_{rkhs}^{{\dagger}}(A^{\top}u_{i+1}-\beta_{i+1}C_{rkhs}z_{i}).$$

Letting $\bar{z}_{i}=C_{rkhs}z_{i}$ and $p=A^{\top}u_{i+1}-\beta_{i+1}\bar{z}_{i}$ again, we get $\alpha_{i+1}=\|C_{rkhs}^{{\dagger}}p\|_{C_{rkhs}}=(p^{\top}C_{rkhs}^{{\dagger}}p)^{1/2}$.

Thus, the two cases of $C_{rkhs}$ lead to the same iterative process. We summarize the iterative process in Algorithm 2, and call it generalized Golub-Kahan bidiagonalization (gGKB).