Recovery guarantees for compressed sensing with unknown errors

TL;DR

This paper establishes robust recovery guarantees for l1-minimization in compressed sensing without requiring prior noise estimates, enhancing practical applicability in high-dimensional and inverse problem contexts.

Contribution

It provides novel theoretical guarantees for l1-minimization robustness when noise estimates are unknown, applicable to various high-dimensional and inverse problem scenarios.

Findings

Robust recovery guarantees for basis pursuit without noise estimates.

Application to high-dimensional function approximation and inverse problems.

Relevance to MRI and other imaging techniques.

Abstract

From a numerical analysis perspective, assessing the robustness of l1-minimization is a fundamental issue in compressed sensing and sparse regularization. Yet, the recovery guarantees available in the literature usually depend on a priori estimates of the noise, which can be very hard to obtain in practice, especially when the noise term also includes unknown discrepancies between the finite model and data. In this work, we study the performance of l1-minimization when these estimates are not available, providing robust recovery guarantees for quadratically constrained basis pursuit and random sampling in bounded orthonormal systems. Several applications of this work are approximation of high-dimensional functions, infinite-dimensional sparse regularization for inverse problems, and fast algorithms for non-Cartesian Magnetic Resonance Imaging.

Equations

$$\mathbf{y}=\mathbf{Ax}+\mathbf{n},$$

$$\hat{\mathbf{x}}(\eta )\in \arg \underset{\mathbf{z}\in {\mathbb{C}}^N}{\min }\Vert \mathbf{z}{\Vert }_1,s.t.\Vert \mathbf{Az}-\mathbf{y}{\Vert }_2\le \eta ,$$

$$\Vert \mathbf{n}{\Vert }_2\le \eta .$$

$$\Vert \mathbf{x}-\hat{\mathbf{x}}(\eta ){\Vert }_2\lesssim \frac{{\sigma }_s(\mathbf{x}{)}_1}{\sqrt{s}}+\eta ,$$

$$\Vert \mathbf{n}{\Vert }_2\ge \eta .$$

$$\Vert \mathbf{x}-\hat{\mathbf{x}}(\eta ){\Vert }_2\lesssim \frac{{\sigma }_s(\mathbf{x}{)}_1}{\sqrt{s}}+\eta +{\mathcal{L}}^{\frac{1}{2}}\max \left\{\Vert \mathbf{n}{\Vert }_2-\eta ,0\right\},$$

$$\Vert {\mathbf{z}}_S{\Vert }_q\le \rho \Vert {\mathbf{z}}_{\overline{S}}{\Vert }_1+\tau \Vert \mathbf{Az}\Vert ,\forall \mathbf{z}\in {\mathbb{C}}^N.$$

$$(1-\delta )\Vert \mathbf{z}{\Vert }_2^2\le \Vert \mathbf{Az}{\Vert }_2^2\le (1+\delta )\Vert \mathbf{z}{\Vert }_2^2,\forall \mathbf{z}\in {\Sigma }_s^N.$$

$$A_{kj}:=m^{-\frac{1}{2}}{\varphi }_j({\tau }_k),\forall k\in [m],j\in [N].$$

$$m\gtrsim s\mathcal{L}(N,s,\delta ,\epsilon ,K),$$

$$\mathcal{L}=\frac{K^2}{{\delta }^2}\max \left\{{\ln }^2(s)\ln \left(\frac{K^2}{{\delta }^2}s\ln (N)\right)\ln (N),\ln \left(\frac{1}{\epsilon }\right)\right\}.$$

$$s_1(\mathbf{M})\ge s_2(\mathbf{M})\ge \cdots \ge s_m(\mathbf{M})=s_{\min }(\mathbf{M}).$$

$$\Vert \hat{\mathbf{x}}(\eta )-\mathbf{x}{\Vert }_2\lesssim \frac{{\sigma }_s(\mathbf{x}{)}_1}{\sqrt{s}}+\eta +\sqrt{\frac{m}{s}}\frac{\max \{\Vert \mathbf{n}{\Vert }_2-\eta ,0\}}{\min \{s_{\min }(\sqrt{\frac{m}{N}}{\mathbf{A}}^{\ast }),1\}},$$

$$\mathbf{M}=[{\mathbf{m}}_1|\cdots |{\mathbf{m}}_m]\in {\mathbb{C}}^{N\times m},$$

$$\mu :=\frac{1}{N^2}\mathbb{E}\underset{k\in [m]}{\max }\sum _{\ell \in [m]\setminus \{k\}}|⟨{\mathbf{m}}_k,{\mathbf{m}}_{\ell }⟩{|}^2.$$

$$\xi :=\mathbb{E}\underset{k\in [m]}{\max }\left|\frac{\Vert {\mathbf{m}}_k{\Vert }_2^2}{N}-1\right|.$$

$$\mathbb{E}\max_{j\in[m]}\bigg{|}s_{j}\bigg{(}\frac{1}{\sqrt{N}}\mathbf{M}\bigg{)}-1\bigg{|}\lesssim\xi+\sqrt{(1+\xi)\mu\ln m},$$

$$N\mu \le (Km{)}^2.$$

$$\xi\leq\min\bigg{\{}D_{1}\sqrt{\frac{m^{2}\ln(m)}{N}},D_{2}\bigg{\}},$$

$$s\le \frac{\epsilon \sqrt{N}}{c\mathcal{L}(N,s,\epsilon ,K){\ln }^{\frac{1}{2}}(N)},$$

$$m=\lceil ds\mathcal{L}(N,s,\epsilon ,K)\rceil .$$

$$\Vert \hat{\mathbf{x}}(\eta )-\mathbf{x}{\Vert }_2\le \frac{C}{\sqrt{s}}{\sigma }_s(\mathbf{x}{)}_1+D\eta +E{\mathcal{L}}^{\frac{1}{2}}\max \{\Vert \mathbf{n}{\Vert }_2-\eta ,0\},$$

$$\phi_{j}(\tau)=\exp\big{(}2\pi\textrm{i}j\tau/N),\quad\forall j,k\in[N]_{0},$$

$${\varphi }_0(\tau )\equiv 1,$$

$${\varphi }_j(\tau )=\sqrt{2}\cos (j\arccos (\tau )),\forall j\in [N-1].$$

$$f(x)=\ln\Big{(}d+1+\sum_{i=1}^{d}x_{i}\Big{)},\quad x\in[-1,1]^{d}.$$

$$y_i=f({\tau }_i)+\mathcal{N}(0,{\zeta }^2),\forall i\in [m].$$

$${\eta }_{opt}:=\Vert \mathbf{A}{\mathbf{x}}_{ref}-\mathbf{y}{\Vert }_2.$$

Full text

I Introduction

In Compressed Sensing (CS) and sparse representations we deal with underdetermined linear systems of equations

[TABLE]

where $\mathbf{A}\in\mathbb{C}^{m\times N}$, with $m\ll N$, is the sensing matrix, $\mathbf{x}\in\mathbb{C}^{N}$ is an unknown signal, and $\mathbf{y}\in\mathbb{C}^{m}$ is the vector of measurements perturbed by noise $\mathbf{n}\in\mathbb{C}^{m}$ [12, 15]. This corruption could be due to physical noise produced by the measuring device, to approximation errors in the model, or to numerical factors. Some examples are model error in inverse problems such as MRI [19, 22], the expansion error in infinite-dimensional CS when truncating the signal to its finite dimensional representation [1, 4], or the quadrature error involved in the evaluation of the bilinear form associated with a PDE [9, 10, 7].

A standard tool to regularize the inverse problem (25) and recover a good approximation $\hat{\mathbf{x}}(\eta)$ to the solution $\mathbf{x}$ (assumed to be sparse or compressible) is the Quadratically Constrained Basis Pursuit (QCBP) optimization program

[TABLE]

also called Basis Pursuit (BP) when $\eta=0$. Usually, in order to study the recovery guarantees of (2), the parameter $\eta$ is assumed to control the noise magnitude, i.e.,

[TABLE]

Indeed, under the regime (3) and with suitable hypotheses on the sensing matrix $\mathbf{A}$ (e.g., based on the restricted isometry property), the following type of recovery error estimate holds with high probability

[TABLE]

where $\sigma_{s}(\mathbf{x})_{1}$ is the best $s$-term approximation error of $\mathbf{x}$ with respect to the $\ell^{1}$-norm [14, 17]. Unfortunately, a priori estimates of the noise of the form (3) may not be available in real applications of CS. Moreover, since the recovery error estimate (4) is sensitive to $\eta$, the choice of this parameter is crucial (see Figure 1). In practice, one could resort to cross-validation in order to tune this parameter, but this technique could be time-consuming or inaccurate and it is not properly understood from a theoretical perspective [16].

The goal of this work is to establish robust recovery guarantees for QCBP (and BP) under the regime

[TABLE]

In this scenario, recovery estimates analogous to (4)–where $\eta$ is replaced by $\|\mathbf{n}\|_{2}$–hold for BP [17]. They are based on the so-called quotient property, which is known to be satisfied only by random Gaussian matrices [26] and by Weibull matrices [18], under suitable restrictions on the number of measurements $m$. Similar robust recovery estimates are also available for algorithms such as iterative hard thresholding, CoSaMP, and orthogonal matching pursuit [17]. Yet, all these techniques require an a priori knowledge of the sparsity level $s$ that is not necessary for QCBP.

Here, we prove robust recovery error estimates for QCBP (and BP) when the matrix $\mathbf{A}$ is built by random sampling from bounded orthonormal systems [20]. In particular, under suitable hypotheses involving the restricted isometry constants and the singular values of $\mathbf{A}$, we provide recovery error estimates in probability of the form (see Theorem 6)

[TABLE]

where the factor $\mathcal{L}$ is polylogarithmic in $N$ and $s$, and can be defined, for example, as in (10). The effect of the unknown error $\mathbf{n}$ is encapsulated in the third term on the right-hand side (compare with (4)). As is to be expected, this term approaches zero as the estimation of the model error $\eta$ improves.

This analysis has several potential extensions and applications. First, the case of weighted $\ell^{1}$-minimization. This is used notably in high-dimensional function approximation and interpolation [2, 21, 23], with applications in uncertainty quantification for parametric PDEs. Second, fast methods for non-Cartesian MRI, where model error arises from gridding non-uniform Fourier data to a uniform grid, and can seriously hamper reconstruction quality [3, 19]. Third, low-rank matrix recovery.

II Tools from compressed sensing

We first recall some concepts from CS that constitute the foundations that our analysis will be built upon: the null space and the restricted isometry properties (Section II-A), and random sampling in bounded orthonormal systems (Section II-B).

In the following, for every $k\in\mathbb{N}$, we define $[k]:=\{1,\ldots,k\}$ and $[k]_{0}:=\{0,\ldots,k-1\}$. Moreover, we denote the set of $s$-sparse vectors in $\mathbb{C}^{N}$ as $\Sigma_{s}^{N}$.

II-A Robust null space and restricted isometry properties

The first tool involved in our analysis is the so-called $\ell^{q}$-robust Null Space Property (NSP) [17, Chapter 4].

Definition 1** ($\ell^{q}$-robust null space property)**

Given $q\geq 1$, the matrix $\mathbf{A}\in\mathbb{C}^{m\times N}$ satisfies the $\ell^{q}$-robust null space property of order $s$ (with respect to the norm $\|\cdot\|$) with constants $0<\rho<1$ and $\tau>0$ if, for any set $S\subseteq[N]$ with $|S|\leq s$, it holds

[TABLE]

Moreover, we recall the well-known restricted isometry property (also known as “RIP”), where the sensing matrix is required to behave similarly to an isometry when its action is restricted to the set of sparse vectors [12].

Definition 2** (Restricted isometry property)**

The $s^{th}$ restricted isometry constant $\delta_{s}$ of a matrix $\mathbf{A}\in\mathbb{C}^{m\times N}$ is the smallest constant $\delta\geq 0$ such that

[TABLE]

The matrix $A$ has the Restricted Isometry Property (RIP) of order $s$ if $0<\delta_{s}<1$.

It is well-known that $\delta_{2s}<4/\sqrt{41}$ is a sufficient condition for the $\ell^{2}$-robust NSP to hold [17, Theorem 6.13] (we decide to use this condition and not that presented in [11] for ease of exposition).

II-B Bounded orthonormal systems

Our analysis focuses on the case of measurement matrices $\mathbf{A}$ arising from random sampling from a Bounded Orthonormal System (BOS) [20]. Some significant examples of random sampling from a BOS are the subsampled Fourier transform, nonharmonic Fourier measurements, and random evaluation of orthogonal polynomials. We will discuss these case studies in more detail in Section III-C.

Definition 3** (Random sampling from a BOS)**

Let $\mathcal{D}\subseteq\mathbb{R}^{d}$ be endowed with probability measure $\nu$. Then, a set $\Phi=\{\phi_{1},\ldots,\phi_{N}\}$ of complex-valued functions on $\mathcal{D}$ is called a Bounded Orthonormal System (BOS) with constant $K$ if, for every $j,k\in[N]$, it holds $\int_{\mathcal{D}}\phi_{j}(\tau)\overline{\phi_{k}(\tau)}d\nu(\tau)=\delta_{jk}$ and $\|\phi_{j}\|_{\infty}:=\sup_{\tau\in\mathcal{D}}|\phi_{j}(\tau)|\leq K$. Moreover, given $m$ independent random variables $\tau_{1},\ldots,\tau_{m}$, distributed according to $\nu$, we define the random sampling matrix $\mathbf{A}\in\mathbb{C}^{m\times N}$ associated with a the BOS $\Phi$ as

[TABLE]

A crucial property of this kind of matrices is the following. For every $\delta\in(0,1)$, assuming

[TABLE]

the $s^{th}$ restricted isometry constant of $\mathbf{A}$ satisfies $\delta_{s}\leq\delta$ with probability at least $1-\varepsilon$, where the factor $\mathcal{L}(N,s,\delta,\varepsilon,K)$ is polylogarithmic in $N$ and $s$ and can be chosen in different ways (see [20, 17, 21, 13]). For example, according to [17, Theorem 12.32], it can be defined as follows

[TABLE]

III Robustness analysis

In this section, we illustrate our robustness analysis. For proofs of the results stated here, we refer to [8].

III-A Singular values of tall random matrices

Given a “tall” matrix $\mathbf{M}\in\mathbb{C}^{N\times m}$, we denote and sort its singular values as follows

[TABLE]

First, we provide a robust error estimate for QCBP as defined in (2) assuming the $\ell^{2}$-robust NSP. This estimate depends on the minimum singular value of the Hermitian conjugate of the sensing matrix $\mathbf{A}$.

Proposition 4** (Robust error estimate based on the NSP)**

Let $\mathbf{A}\in\mathbb{C}^{m\times N}$ be a matrix of rank $m$ that satisfies the $\ell^{2}$-robust NSP of order $s$ with constants $0<\rho<1$ and $\tau>0$ with respect to $\|\cdot\|_{2}$ (see Definition 1). Then, the following error estimate holds for problem (2)

[TABLE]

where the hidden constant depends on $\rho$ and $\tau$.

Proposition 4 shows that, in order to get robust error estimates for QCBP, we need to understand the asymptotic behavior of $s_{\min}(\sqrt{\frac{m}{N}}\mathbf{A}^{*})$. In particular, our goal is to show that $s_{\min}(\sqrt{\frac{m}{N}}\mathbf{A}^{*})\approx 1$. With this aim, we employ principles and ideas from the theory of random matrices with isotropic heavy-tailed columns [25].

Recall that a random vector $\mathbf{z}\in\mathbb{C}^{N}$ is said to be isotropic if $\mathbb{E}[\mathbf{z}\mathbf{z}^{*}]=\mathbf{I}$. Then, consider random matrices of the form

[TABLE]

where the columns $\mathbf{m}_{j}$ are independent random isotropic vectors. We introduce the cross-coherence parameter, defined as

[TABLE]

This parameter controls the off-diagonal part of the Gram matrix of $\mathbf{M}$. We also define the distortion parameter as

[TABLE]

It measures how far the columns of $\mathbf{M}$ are from being normalized. We notice that if $\|\mathbf{m}_{k}\|_{2}=\sqrt{N}$ almost surely for every $k\in[m]$, then $\xi$ vanishes. Using this two parameters, we give a generalization of [25, Theorem 5.62].

Theorem 5** (Singular values of heavy-tailed matrices)**

Let $\mathbf{M}$ be an $N\times m$ matrix ($N\geq m$) whose columns are independent isotropic random vectors in $\mathbb{C}^{N}$. Then, its singular values satisfy the following asymptotic estimate

[TABLE]

where $\mu$ and $\xi$ are defined as in (13) and (14), respectively.

We observe that a necessary condition for estimate (15) to be nontrivial is that the quantity $\xi$ should be bounded uniformly in $N$ (but not necessarily in $m$). Regarding the incoherence parameter, exploiting the isotropy of the columns of $\mathbf{M}$, and assuming $\|\mathbf{m}_{j}\|_{2}\leq K\sqrt{N}$ for a suitable constant $K>0$ (this is always the case for $\mathbf{M}=\sqrt{m}\mathbf{A}^{*}$, where $\mathbf{A}$ is the random sampling matrix associated with a BOS), we can prove that

[TABLE]

We check the sharpness of this upper bound for moderate values of $m$ numerically in the case of the subsampled Fourier transform, where $K=1$ (see Figure 2).

III-B Recovery error estimate

The following theorem is the main result of the paper. It provides a robust error estimate in probability for QCBP, when the solution $\mathbf{x}$ is very sparse and $\mathbf{A}$ is a random sampling matrix associated with a BOS.

Theorem 6** (Robust recovery error estimate for QCBP)**

Consider a BOS $\Phi$ with constant $K\geq 1$ and with a distortion parameter that satisfies

[TABLE]

for suitable constants $D_{1},D_{2}$ independent of $m$ and $N$. Then, there exist constants $c,d,C,D,E>0$ and a function $\mathcal{L}=\mathcal{L}(N,s,\varepsilon,K)$ depending polylogarithmically on $N$ and $s$ such that the following holds. For every $N\in\mathbb{N}$ and $\varepsilon\in(0,1)$, assume that the sparsity level satisfies

[TABLE]

and let $\mathbf{A}\in\mathbb{C}^{m\times N}$ be the random sampling matrix associated with $\Phi$ and with a number of measurements

[TABLE]

Then, the following robust error estimate holds for QCBP

[TABLE]

with probability at least $1-\varepsilon$. The constant $c$ depends on $D_{1}$, $D_{2}$, and $K$, whereas the constants $d$, $C$, $D$, and $E$ are universal. The function $\mathcal{L}$ can be defined as in (10), with $\delta=1/2$.

The required relation between $s$ and $m$ is linear up to logarithmic factors, in accordance with the usual recovery error estimate in CS. However, there are three main limitations of Theorem 6 that are worth underlining. First, the result holds for a particular sparsity regime (18). Essentially, fixed the failure probability $\varepsilon$, we require $s\lesssim\sqrt{N}$ (up to logarithmic factors). Second, the error estimate (20) actually depends on $N$, but this dependence is only polylogarithmic (the factor $\mathcal{L}^{\frac{1}{2}}$ is due to the term $\sqrt{m/s}$ in the error estimate of Proposition 4). Third, there is a linear dependence between $s$ and $\varepsilon$ in (18). Thus, the failure probability of the estimate is not “overwhelmingly low”. These three issues are open problems currently under investigation.

III-C Applications

To conclude, we discuss some applications of Theorem 6 to concrete examples from signal processing and high-dimensional polynomial approximation.

III-C1 Fourier and Chebyshev BOSs

We discuss two examples of BOSs very popular in CS. First, we consider the subsampled Fourier transform. We have $\mathcal{D}=[N]_{0}$, the system is defined as

[TABLE]

and the sampling measure $\nu$ is the uniform discrete distribution on $\mathcal{D}$. It turns out that $K=1$ and, consequently, the distortion parameter is $\xi=0$. As a result, condition (17) of Theorem 6 is satisfied.

In the case of the Chebyshev system, we consider the Chebyshev orthogonal polynomials on $\mathcal{D}=[-1,1]$, defined as

[TABLE]

[TABLE]

They form a BOS with respect to the Chebyshev measure $d\nu(\tau)=\pi^{-1}(1-\tau^{2})^{-1/2}$ on $\mathcal{D}$ and with constant $K=\sqrt{2}$. By studying the normalized Christoffel function associated with this system, we estimate that the distortion parameter $\xi$ decays proportionally to $\sqrt{m/N}$. Therefore, hypothesis (17) holds true and we can apply Theorem 6.

III-C2 High-dimensional polynomial approximation

We assess the robustness of QCBP for polynomial approximation in high dimension [1]. Consider the multivariate function

[TABLE]

We fix $d=10$ and we employ the tensorized version of the Chebyshev polynomials defined in (22)-(23) over $[-1,1]^{10}$ as a sparsity basis. We set maximum degree $10$ on each variable and we restrict the multi-index space to the hyperbolic cross shape [5]. These choices leads to a total of $N=581$ degrees of freedom. We evaluate $f$ at $m=50$ random independent sampling points $\tau_{1},\ldots,\tau_{m}$ identically distributed according to the tensorized Chebyshev measure over $\mathcal{D}=[-1,1]^{10}$. In order to assess the robustness to unknown error, we artificially add centered gaussian noise with standard deviation $\zeta$ to the measurements, namely

[TABLE]

In the case of nonintrusive methods for the uncertainty quantification of PDEs with random parameters, we can interpret this noise as the numerical error associated with the black-box PDE solver used to produce point-wise samples of the quantity of interest [16].

In Figure 3, we plot the absolute error $\|\hat{\mathbf{x}}(\eta)-\mathbf{x}\|_{2}$ as a function of $\eta$ for different values of the standard deviation $\zeta$. The resulting curve always exhibits a global minimum. We observe that, for $\zeta=0.1$, underestimating $\eta$ is better than overestimating it. For $\zeta=1$, the minimum becomes more pronounced and underestimating $\eta$ becomes more penalizing. Recalling (20) in Theorem 6, this behavior could be justified as follows: for small values of $\zeta$, the recovery error is dominated by the term $s^{-\frac{1}{2}}\sigma_{s}(\mathbf{x})_{1}$, whereas, the more $\zeta$ gets larger, the more the term $\mathcal{L}^{\frac{1}{2}}(\|\mathbf{n}\|_{2}-\eta)$ becomes dominant.

In order to estimate the value of $\eta$ that minimizes the error, we evaluate the residual on a reference solution $\mathbf{x}_{ref}$ computed via least-square fitting over an oversampled random grid of size $40N=23240$, i.e.,

[TABLE]

We compare this value with the one computed by cross-validation [16]. Both approaches are able to approximate the minimum quite well for $\zeta=0.1$, whereas, in the case $\zeta=1$, cross-validation slightly underperforms.

III-C3 Non-Cartesian Magnetic Resonance Imaging

Finally, in Figure 4 we give an application of this analysis to non-Cartesian MRI. For fast reconstruction in sparse MRI, non-Cartesian data is often preprocessed by gridding it to a uniform integer grid, thus ensuring that the sampling matrix $\mathbf{A}$ can be expressed as a subsampled DFT matrix. This introduces $O(1)$ model errors, which, as seen in Figure 4, adversely affect the reconstruction. On the other hand, the gridding strategy introduced in [3] leads to model errors of order $1/n_{up}$, where $n_{up}$ is a user-controlled parameter. Theorem 6 theoretically establishes the advantage of this higher-fidelity gridding, as verified in Figure 4.