Monte Carlo wavelets: a randomized approach to frame discretization

TL;DR

This paper introduces Monte Carlo wavelets, a stochastic discretization method for continuous wavelets on general domains, leveraging spectral calculus and concentration of measure to ensure convergence and establish rates.

Contribution

It presents a novel randomized approach to discretize continuous wavelets called Monte Carlo wavelets, with theoretical guarantees of convergence.

Findings

Convergence of Monte Carlo wavelet discretization established

Derived convergence rates under regularity assumptions

Applicable to general domains using reproducing kernel Hilbert spaces

Abstract

In this paper we propose and study a family of continuous wavelets on general domains, and a corresponding stochastic discretization that we call Monte Carlo wavelets. First, using tools from the theory of reproducing kernel Hilbert spaces and associated integral operators, we define a family of continuous wavelets by spectral calculus. Then, we propose a stochastic discretization based on Monte Carlo estimates of integral operators. Using concentration of measure results, we establish the convergence of such a discretization and derive convergence rates under natural regularity assumptions.

Tables

Table 1. Notation

symbol	definition
$⟨\cdot ,\cdot ⟩$	inner product in a Hilbert space
$\parallel \cdot \parallel$	norm of a Hilbert space, or the operator norm
$\parallel \cdot {\parallel }_{\text{HS}}$	Hilbert-Schmidt norm
$\mathrm{span}\left\{S\right\}$	linear span of the set $S$
$\overline{S}$	topological closure of the set $S$
$P_S$	orthogonal projection onto the closed subspace $S$
$\mathrm{supp}\left(\rho \right)$	support of the measure $\rho$
${\delta }_x$	Dirac measure at $x$
$X\lesssim Y$	$X\le CY$ for some constant $C>0$
$𝐮\left[i\right]$	$i$-th component of the vector $𝐮\in {ℝ}^N$
$𝐌[i,j]$	$(i,j)$-th component of the matrix $𝐌\in {ℝ}^{N\times N}$

Table 2. Table 1:

method	$g_j\left(\lambda \right)$	qualification
Tikhonov regularization	$\frac{1}{\lambda +1/j}$	$1$
Iterated Tikhonov	$\frac{{\left(\lambda +1/j\right)}^m-{\left(1/j\right)}^m}{\lambda {\left(\lambda +1/j\right)}^m}$	$m$
Landweber iteration	$\frac{1}{\lambda }\left(1-{\left(1-\gamma \lambda \right)}^j\right)$	$\infty$
Asymptotic regularization	$\frac{1}{\lambda }\left(1-\exp \left(-j\lambda \right)\right)$	$\infty$

Table 3. Table 2:

method	rate
Tikhonov regularization	$N^{-\frac{\min \{\alpha ,1\}}{2\min \{\alpha ,1\}+2}}$
Iterated Tikhonov	$N^{-\frac{\min \{\alpha ,m\}}{2\min \{\alpha ,m\}+2}}$
Landweber iteration	$N^{-\frac{\alpha }{2\alpha +2}}$
Asymptotic regularization	$N^{-\frac{\alpha }{2\alpha +2}}$

Equations

$$K:\mathcal{X}\times \mathcal{X}\to \mathbb{C},$$

$$\mathcal{H}=\{f:\mathcal{X}\to \mathbb{C}:f(x)=⟨f,K_x⟩\}=\overline{\mathrm{span}}\{K_x{\}}_{x\in \mathcal{X}},$$

$$Tf(x)={\int }_{\mathcal{X}}K(x,y)f(y)d\rho (y).$$

$$T{v}_i=\{\begin{array}{ll}{\lambda }_i{v}_i& i\in I_{\rho }\\ 0v_i& i\in I_0\end{array}$$

$${\mathcal{H}}_{\rho }=\overline{\mathrm{span}}\{K_x{\}}_{x\in \mathrm{supp}(\rho )}=\overline{\mathrm{span}}\{v_i{\}}_{i\in I_{\rho }}$$

$$\sum _{j\in \mathbb{N}}\lambda G_j(\lambda {)}^2=1\lambda \in (0,{\kappa }^2].$$

$${\psi }_{j,x}=G_j(T)K_{x}j\in \mathbb{N},x\in \mathcal{X}.$$

$${\psi }_{j,x}=\sum _{i\in I_{\rho }}{G}_j({\lambda }_i)\overline{v_i(x)}{v}_i.$$

$$\Vert P_{{\mathcal{H}}_{\rho }}f{\Vert }^2=\sum _{j\in \mathbb{N}}{\int }_{\mathcal{X}}|⟨f,{\psi }_{j,x}⟩{|}^{2}d\rho (x).$$

$${\int }_{\mathcal{X}}|⟨f,{\psi }_{j,x}⟩{|}^{2}d\rho (x)$$

$$=⟨T{G}_j(T{)}^{2}f,f⟩$$

$$=\sum _{i\in I_{\rho }}{\lambda }_i{G}_j({\lambda }_i{)}^2|⟨f,v_i⟩{|}^2.$$

$$\sum _{j\in \mathbb{N}}⟨G_j(T{)}^{2}Tf,f⟩=\sum _{i\in I_{\rho }}|⟨f,v_i⟩{|}^2\sum _{j\in \mathbb{N}}{\lambda }_i{G}_j({\lambda }_i{)}^2=\sum _{i\in I_{\rho }}|⟨f,v_i⟩{|}^2=\Vert P_{{\mathcal{H}}_{\rho }}f{\Vert }^2.$$

$$0\le g_j\le g_{j+1},\underset{j\to \infty }{\lim }\lambda g_j(\lambda )=1.$$

$$G_0(\lambda )=\sqrt{g_0(\lambda )},G_{j+1}(\lambda )=\sqrt{g_{j+1}(\lambda )-g_j(\lambda )}.$$

$$\underset{\lambda \in (0,{\kappa }^2]}{\sup }{\lambda }^{\nu }(1-\lambda g_j(\lambda ))\lesssim j^{-\min \{\nu ,\overline{\nu }\}}.$$

$$x_1,\dots ,x_N\sim \rho .$$

$${\hat{\rho }}_N=\frac{1}{N}\sum _{i=1}^N{\delta }_{x_i},$$

$$\hat{T}f(x)=\frac{1}{N}\sum _{i=1}^{N}K(x,x_i)f(x_i),$$

$${\hat{\psi }}_{j,x_i}=G_j(\hat{T})K_{x_i}j\in \mathbb{N},i=1,\dots ,N.$$

$${\hat{\mathcal{H}}}_N={\mathcal{H}}_{{\hat{\rho }}_N}=\mathrm{span}\{K_{x_i}{\}}_{i=1}^N\simeq {\mathbb{R}}^N.$$

$${\hat{\mathcal{H}}}_N\subset {\hat{\mathcal{H}}}_{N+1}\subset \cdots \subset \mathcal{H}.$$

$$T_{j}f={\int }_{\mathcal{X}}⟨f,{\psi }_{j,x}⟩{\psi }_{j,x}d\rho (x),{\hat{T}}_{j}f=\frac{1}{N}\sum _{i=1}^N⟨f,{\hat{\psi }}_{j,x_i}⟩{\hat{\psi }}_{j,x_i}$$

$$I_{\mathcal{H}}=\sum _{j\in \mathbb{N}}{T}_j\approx \sum _{j\le \tau }{\hat{T}}_j.$$

$$\Vert f-\sum _{j\le \tau }{\hat{T}}_{j}f\Vert \le \Vert \sum _{j>\tau }{T}_{j}f\Vert +\Vert \sum _{j\le \tau }(T_j-{\hat{T}}_j)f\Vert ,$$

$$\Vert \sum _{j>\tau }{T}_{j}f\Vert \stackrel{\tau \to \infty }{\to }0$$

$$f=T^{\alpha }h\text{for some }h\in \mathcal{H},\alpha >0.$$

$$\Vert \sum _{j>\tau }{T}_{j}f\Vert \lesssim \Vert T^{-\alpha }f\Vert {\kappa }^{2(\alpha -\beta )}{\tau }^{-\beta },$$

$$\Vert \sum _{j>\tau }{T}_{j}f{\Vert }^2$$

$$=\sum _i{\lambda }_i^{2\beta }|1-{\lambda }_i{g}_{\tau }({\lambda }_i){|}^2|⟨T^{-\beta }f,v_i⟩{|}^2.$$

$$\Vert \sum _{j>\tau }{T}_{j}f\Vert$$

$$\Vert \sum _{j\le \tau }(T_j-{\hat{T}}_j)f\Vert \lesssim \Vert f\Vert {\kappa }^2\sqrt{t}\tau N^{-1/2}.$$

Full text

1 Introduction

To construct a discrete frame it is common practice to first construct a frame on a continuous parameter (measure) space, where the mathematical properties are nice and the structure is rich, and then obtain a discrete frame by carefully selecting a discrete subset of parameters. The discretization problem has been studied extensively and in several settings [5, 2, 7, 6]. Although theoretically relevant in establishing density theorems and complete characterizations, some general discretization procedures require assumptions which might be difficult to confirm, and they are often not constructive. On the other hand, discrete parameter selection in constructive discretization methods is usually sensitive, and it does not generalize trivially from the Euclidean space to more general geometries [1]. In the following, we propose a different approach based on random sampling.

First, we consider a general domain with an associated positive definite kernel, and exploit the theory of reproducing kernel Hilbert spaces (RKHS) and corresponding integral operators to define a family of continuous wavelets. Here, we borrow ideas from [8, 9, 2], and develop them using the properties of an RKHS. Then, we propose a stochastic discretization approach which replaces constructive discrete parameter selection with random sampling. We note that the corresponding discrete frame is actually only finite-dimensional. In particular, it does not solve the discretization problem in the classical sense. Indeed, each discrete frame is a frame only up to a certain sampling resolution. The frame dimensions increase with the number of samples, and our main technical contribution is proving convergence to the continuous frame as the number of samples goes to infinity. Further, we derive convergence rates under natural regularity assumptions [4]. Our analysis combines classical tools from approximation theory and spectral calculus with results of concentration of measure to deal with random sampling.

The paper is organized as follows. In Section 2 we introduce the general framework in which we will conduct our analysis. This will allow us to state the continuous and discrete settings within a unified formalism. The key ingredient is a reproducing kernel and the associated integral operator. Section 3 contains the general frame construction based on spectral filters and eigenfunctions of the integral operator. In Section 4 we consider a specific class of filters, stemming from the theory of inverse problems, for which we will derive explicit approximation rates. In Section 5 we specialize the construction of Section 3 in a continuous and a discrete frame, regarding the latter as a Monte Carlo estimate of the former. Our main result, Theorem 4, establishes quantitative consistency of the Monte Carlo estimate on a class of Sobolev-regular signals. In Section 6 we provide an implementable formula to compute our Monte Carlo wavelets from the eigendecomposition of a sample kernel matrix.

2 Setting

Let ${\mathcal{X}}$ be a locally compact, second countable topological space endowed with a Borel probability measure $\rho$. We fix a continuous positive semi-definite kernel

[TABLE]

and let ${\mathcal{H}}$ be the corresponding reproducing kernel Hilbert space

[TABLE]

with $K_{x}=K(\cdot,x)$ and inner product given by $\langle K_{x},K_{y}\rangle=K(x,y)$. Assuming $\sup_{x\in X}K(x,x)\leq{\kappa}^{2}$, we define the integral operator $T:{\mathcal{H}}\to{\mathcal{H}}$

[TABLE]

Under these assumptions, $T$ is a positive self-adjoint trace class operator with spectrum $\sigma(T)\subset[0,{\kappa}^{2}]$. Let

[TABLE]

be the spectral decomposition of $T$ over the orthonormal basis $\{v_{i}\}_{i\in I_{\rho}\cup I_{0}}$. The subspace

[TABLE]

splits ${\mathcal{H}}$ as ${\mathcal{H}}={\mathcal{H}}_{\rho}\oplus\ker T$.

As examples of this setting, we can think of ${\mathcal{X}}$ as ${\mathbb{R}}^{d}$, or non-Euclidean domains such as a compact connected Riemannian manifold or a weighted graph. In all these cases, we can take $K$ as the heat kernel associated with the proper notion of Laplacian, be it the Laplace-Beltrami operator or the graph Laplacian.

3 Frame construction

To construct our frame, we take inspiration from [8] and references therein. We start by defining filters on the spectrum of $T$. Let $(G_{j})_{{j}\in{\mathbb{N}}}$ be a family of bounded measurable functions $G_{j}:[0,\infty)\to[0,\infty)$ satisfying

[TABLE]

$G_{j}(T)$ is then a positive, self-adjoint operator on ${\mathcal{H}}$ with $\sigma(G_{j}(T))=G_{j}(\sigma(T))\\ \subset[0,G_{j}({\kappa}^{2})]$, for every ${j}$. We can thus define

[TABLE]

Note that, if $x\in\operatorname{supp}(\rho)$, then $\psi_{{j},x}\in{\mathcal{H}}_{\rho}$ since $G_{j}(T){\mathcal{H}}_{\rho}\subset{\mathcal{H}}_{\rho}$. Moreover, using self-adjointness of $G_{j}(T)$, the reproducing property (2), and expression (4), we can rewrite $\psi_{{j},x}$ as

[TABLE]

This reveals (6) as a generalization of classical wavelets, where $v_{i}$’s correspond to eigenfunctions of the Laplacian. We may therefore interpret $x$ and ${j}$ as location and scale parameters, respectively.

The following proposition shows that the family defined through (6) is a Parseval frame on ${\mathcal{H}}_{\rho}$.

Proposition 1**.**

For every $f\in{\mathcal{H}}$,

[TABLE]

Proof.

Let $f\in\operatorname{span}\{v_{i}\}_{i\in I_{\rho}\cup I_{0}}$. Since $G_{j}(T)$ is self-adjoint, we have $\langle f,\psi_{{j},x}\rangle=\langle G_{j}(T)f,K_{x}\rangle$. Hence, integrating over ${\mathcal{X}}$ we get

[TABLE]

Summation over ${\mathbb{N}}$ along with (5) thus gives

[TABLE]

The assertion follows since $\operatorname{span}\{v_{i}\}_{i\in I_{\rho}\cup I_{0}}$ is dense in ${\mathcal{H}}$. ∎

4 Examples of spectral filters

Let $(g_{j})_{{j}\in{\mathbb{N}}}$ be a family of functions $g_{j}:[0,\infty)\to[0,\infty)$ such that

[TABLE]

We can define filters obeying (5) by means of (8), as

[TABLE]

The qualification of a spectral function $g_{j}({\lambda})$ is a constant $\overline{\nu}\in(0,\infty]$ such that, for all $\nu>0$,

[TABLE]

As instances of $g_{j}({\lambda})$, we can take spectral functions in Table 1 (see also [3]).

These functions play a fundamental role in ill-posed inverse problems regularization, where the error of a given regularizer decays with a rate governed by its qualification.

5 Consistency

From now on, assume $\operatorname{supp}{\rho}={\mathcal{X}}$, whence ${\mathcal{H}}_{\rho}={\mathcal{H}}$ by (4). Therefore, (6) is a frame on the whole space ${\mathcal{H}}$, thanks to Proposition 1. Suppose we are given $N$ independent and identically distributed samples

[TABLE]

Then we can consider the empirical distribution

[TABLE]

define the empirical integral operator

[TABLE]

and repeat the construction of Section 3, to get

[TABLE]

In view of Proposition 1, (11) defines a discrete Parseval frame on

[TABLE]

Since (10) is a Monte Carlo estimate of (3), we call (11) a family of Monte Carlo wavelets. Note that, as we take more and more samples, we obtain a sequence of discrete Parseval frames on a sequence of nested subspaces of increasing dimension:

[TABLE]

Our goal is to establish consistency of the Monte Carlo estimate (11), that is, to see whether, and in what sense, (11) tends to (6) as $N\to\infty$.

Now, let

[TABLE]

be the frame operator at scale ${j}$ and its empirical counterpart. Proposition 1 gives a resolution of the identity, which can be approximated by a truncated empirical version:

[TABLE]

We split the resolution error of (13) on a signal $f\in{\mathcal{H}}$ into the approximation and the sampling error,

[TABLE]

and derive quantitative bounds for each term (Propositions 2 and 3). Tuning the resolution $\tau$ in terms of the sample size $N$ will then yield our consistency result in Theorem 4.

For what concerns the approximation error, note that Proposition 1 already implies

[TABLE]

by dominated convergence. However, in order to obtain a convergence rate, we need to assume some regularity of the signal $f$. Specifically, we will assume the following Sobolev smoothness condition [4]:

[TABLE]

Proposition 2** (Approximation error).**

Let $g_{j}({\lambda})$ as in (8). Suppose $g_{j}({\lambda})$ has qualification $\overline{\nu}$. Then, for every $f\in T^{\alpha}{\mathcal{H}}$,

[TABLE]

where ${\beta}=\min\{{\alpha},\overline{\nu}\}$.

Proof.

By Lemma 5 we have $\sum_{{j}\leq\tau}T_{j}=Tg_{\tau}(T)$, hence

[TABLE]

Thus,

[TABLE]

Bound on the sampling error follows by concentration of the empirical integral operator $\widehat{T}$ on the continuous operator $T$.

Proposition 3** (Sampling error).**

Let $g_{j}({\lambda})$ be as in (8). Suppose ${\lambda}\mapsto{\lambda}g_{j}({\lambda})$ is Lipschitz continuous on $[0,{\kappa}^{2}]$ with Lipschitz constant bounded by $j$. Then, for every $f\in{\mathcal{H}}$, with probability higher than $1-2e^{-t}$,

[TABLE]

Proof.

Using Lemma 6 we have

[TABLE]

Thanks to the reproducing property (2), we can write

[TABLE]

This allows for the concentration in Lemma 7, which yields

[TABLE]

with probability higher than $1-2e^{-t}$ . ∎

We are finally ready to state the consistency of our Monte Carlo discretization.

Theorem 4**.**

Let $g_{j}({\lambda})$ as in (8). Suppose $g_{j}({\lambda})$ has qualification $\overline{\nu}$, and ${\lambda}\mapsto{\lambda}g_{j}({\lambda})$ is Lipschitz continuous on $[0,{\kappa}^{2}]$ with Lipschitz constant bounded by $j$. Then, for every $f\in T^{\alpha}{\mathcal{H}}$, with probability higher than $1-2e^{-t}$,

[TABLE]

where ${\beta}=\min\{{\alpha},\overline{\nu}\}$.

Proof.

We split $\|f-\sum_{{j}\leq\tau}{\widehat{T}}_{j}f\|$ in (14), and combine the bounds in Propositions 2 and 3. In order to balance the two terms, we choose $\tau=N^{\frac{1}{2{\beta}+2}}$, which concludes the proof. ∎

In particular, in view of Lemma 8, we can use the spectral filters defined in Table 1, for which we obtain the following rates:

6 Numerical implementation

To implement the system of Monte Carlo wavelets in (11) we exploit equation (7), that is,

[TABLE]

where ${\widehat{T}}{\widehat{v}}_{i}={\widehat{{\lambda}}}_{i}{\widehat{v}}_{i}$ is the spectral decomposition of ${\widehat{T}}$. In particular, we need to compute eigenvalues and eigenvectors of ${\widehat{T}}$, translating the calculations from ${\widehat{{\mathcal{H}}}}$ to ${\mathbb{R}}^{N}$ and back [10]. To this end, we introduce the restriction map ${\widehat{S}}:{\widehat{{\mathcal{H}}}}\to{\mathbb{R}}^{N}$

[TABLE]

whose adjoint ${\widehat{S}}^{*}:{\mathbb{R}}^{N}\to{\widehat{{\mathcal{H}}}}$ defines the out-of-sample extension

[TABLE]

Let

[TABLE]

be the (normalized) sample kernel matrix. Now ${\widehat{T}}={\widehat{S}}^{*}{\widehat{S}}$ and $\widehat{\mathbf{K}}={\widehat{S}}{\widehat{S}}^{*}$, hence ${\widehat{T}}$ and $\widehat{\mathbf{K}}$ share the same eigenvalues. Moreover, the eigenvectors of ${\widehat{T}}$ are related to the ones of $\widehat{\mathbf{K}}$ by the singular value decomposition of ${\widehat{S}}$:

[TABLE]

Thus, we have

[TABLE]

which evaluated at $x_{k}$ gives

[TABLE]

We therefore obtain

[TABLE]

In summary, given a kernel (1), filters $\eqref{eq:G}$ and samples (9), a numerical implementation of the Monte Carlo wavelets (11) can be performed from the sample kernel matrix (16) alone following the steps below:

compute eigenvalues and eigenvectors $({\widehat{{\lambda}}}_{i},{\widehat{u}}_{i})_{i=1}^{N}$ of ${\widehat{K}}$; 2. 2.

apply equation (17).

Lemmata

Lemma 5**.**

Let $T$, $G_{j}$ and $T_{j}$ be defined as in (3), (5) and (12). Then $T_{j}=TG_{j}(T)^{2}$.

Proof.

For every $i\in I$ we have

[TABLE]

Thus,

[TABLE]

Since $\overline{\operatorname{span}}\{v_{i}\}_{i\in I}={\mathcal{H}}$, the claim follows by polarization. ∎

Lemma 6**.**

Let $A,B$ be self-adjoint operators on a Hilbert space ${\mathcal{H}}$, and let $F:{\mathbb{R}}\to{\mathbb{C}}$ be a Lipschitz continuous function with Lipschitz constant $L$. Then

[TABLE]

Proof.

Let $\{e_{i}\}_{i\in I}$ and $\{f_{j}\}_{j\in J}$ be orthonormal bases of ${\mathcal{H}}$ such that $Ae_{i}={\lambda}_{i}e_{i}$ and $Bf_{j}=\mu_{j}f_{j}$. Then

[TABLE]

Lemma 7** ([10, Theorem 7]).**

The operators $T$ and ${\widehat{T}}$ defined in (3) and (10) are Hilbert-Schmidt and

[TABLE]

with probability higher than $1-2e^{-t}$.

Lemma 8**.**

Let $g_{j}({\lambda})$ be as in Table 1. The function ${\lambda}\mapsto{\lambda}g_{j}({\lambda})$ is Lipschitz continuous on $[0,{\kappa}^{2}]$, with Lipschitz constant bounded by ${j}$.

Proof.

The claim follows by a direct computation of the derivative of ${\lambda}\mapsto{\lambda}g_{j}({\lambda})$. ∎