Compressed sensing with sparse corruptions: Fault-tolerant sparse collocation approximations

TL;DR

This paper develops a robust compressive sensing method to recover sparse coefficients in Polynomial Chaos Expansions despite measurement corruptions, addressing hardware/software failures in complex computational frameworks.

Contribution

It introduces a novel theoretical analysis and an iterative reweighted algorithm for fault-tolerant sparse recovery in polynomial chaos expansions.

Findings

The method effectively recovers sparse coefficients with corrupted measurements.

The iterative algorithm improves regularization parameter selection.

Numerical tests demonstrate robustness in high-dimensional differential equation solutions.

Abstract

The recovery of approximately sparse or compressible coefficients in a Polynomial Chaos Expansion is a common goal in modern parametric uncertainty quantification (UQ). However, relatively little effort in UQ has been directed toward theoretical and computational strategies for addressing the sparse corruptions problem, where a small number of measurements are highly corrupted. Such a situation has become pertinent today since modern computational frameworks are sufficiently complex with many interdependent components that may introduce hardware and software failures, some of which can be difficult to detect and result in a highly polluted simulation result. In this paper we present a novel compressive sampling-based theoretical analysis for a regularized $\ell^1$ minimization algorithm that aims to recover sparse expansion coefficients in the presence of measurement corruptions. Our recovery results are uniform, and prescribe algorithmic regularization parameters in terms of a user-defined a priori estimate on the ratio of measurements that are believed to be corrupted. We also propose an iteratively reweighted optimization algorithm that automatically refines the value of the regularization parameter, and empirically produces superior results. Our numerical results test our framework on several medium-to-high dimensional examples of solutions to parameterized differential equations, and demonstrate the effectiveness of our approach.

Tables

Table 1. Table 1: Notation used throughout this article.

$m$	number of measurements
$N$	length of sparse vector
$x$	sparse vector in ${ℂ}^N$
$c$	corruptions vector in ${ℂ}^m$
$A$	$m\times N$ measurement matrix
$n$	noise vector in ${ℂ}^m$
$ϵ$	noise bound
$\lambda$	non-negative weighting parameter for the corruptions vector
$\hat{x}$, $\hat{c}$	solutions of the optimization problem
$S$	subset of $\{1,\dots ,N\}$, indices corresponding to $x$
$T$	subset of $\{1,\dots ,m\}$, indices corresponding to $c$
$s$	sparsity of $x$
$k$	sparsity of $c$
${\Sigma }_s$	set of $s$-sparse vectors in ${ℂ}^N$
${\Sigma }_k$	set of $k$-sparse vectors in ${ℂ}^m$
${\sigma }_s{\left(x\right)}_1$	best $s$-term approximation error, measured in the ${\ell }^1$ norm
${\sigma }_k{\left(c\right)}_1$	best $k$-term approximation error, measured in the ${\ell }^1$ norm

Equations

$$f(\xi )=\sum _{j=1}^N{x}_j{\varphi }_j(\xi )+n(\xi ),$$

$${\sigma }_s(x{)}_p=\underset{\Vert \tilde{x}{\Vert }_0\le s}{\inf }{\Vert x-\tilde{x}\Vert }_p,$$

$$y_k$$

$$y=Ax+c+n\in {\mathbb{R}}^m.$$

$$\min_{z\in\mathbb{C}^{N},d\in\mathbb{C}^{m}}\|z\|_{1}+\lambda\|d\|_{1}\ \mbox{subject to $\|Az+d-y\|_{2}\leq\epsilon\sqrt{m}$}.$$

$$\Vert \hat{x}-x{\Vert }_2+\Vert c-\hat{c}{\Vert }_2$$

$${\delta }_{2s,2k}$$

$${\Vert x-\hat{x}\Vert }_1+\lambda {\Vert c-\hat{c}\Vert }_1$$

$${\Vert x-\hat{x}\Vert }_2+{\Vert c-\hat{c}\Vert }_2$$

$$K:=\underset{j=1,\dots ,N}{\max }\underset{\xi \in D}{\sup }\left|{\varphi }_j(\xi )\right|<\infty ,$$

$$\begin{array}{rl}m& \gtrsim {\delta }^{-2}\cdot K^2\cdot s\cdot \left({\log }^3(2s)\cdot \log (2N)+\log {ϵ}^{-1}\right),\\ m& \gtrsim {\delta }^{-2}\cdot K\cdot s\cdot k,\end{array}$$

$$m\gtrsim {\delta }^{-2}\cdot K\cdot k,$$

$$\frac{{\sigma }_k(c^{\prime }{)}_1}{\sqrt{k}}\le \frac{1}{\sqrt{km}}\left({\sigma }_k(c{)}_1+\Vert n{\Vert }_1\right)\le \frac{{\sigma }_k(c{)}_1}{\sqrt{km}}+\sqrt{\frac{m}{k}}ϵ.$$

$$\Vert x_S{\Vert }_1+\lambda \Vert c_T{\Vert }_1\le \rho \left(\Vert x_{S^c}{\Vert }_1+\lambda \Vert c_{T^c}{\Vert }_1\right)+\tau \Vert Ax+c{\Vert }_2,\forall x\in {\mathbb{C}}^N,c\in {\mathbb{C}}^m,$$

$$\sqrt{\Vert x_S{\Vert }_2^2+\Vert c_T{\Vert }_2^2}\le \frac{\rho }{\sqrt{s+{\lambda }^{2}k}}\left(\Vert x_{S^c}{\Vert }_1+\lambda \Vert c_{T^c}{\Vert }_1\right)+\tau \Vert Ax+c{\Vert }_2,\forall x\in {\mathbb{C}}^N,c\in {\mathbb{C}}^m,$$

$$\Vert x_S{\Vert }_1+\lambda \Vert c_T{\Vert }_1\le \sqrt{s}\Vert x_S{\Vert }_2+\lambda \sqrt{k}\Vert c_T{\Vert }_2\le \sqrt{s+{\lambda }^{2}k}\sqrt{\Vert x_S{\Vert }_2^2+\Vert c_T{\Vert }_2^2}.$$

$$\min_{z\in\mathbb{C}^{N},d\in\mathbb{C}^{m}}\|z\|_{1}+\lambda\|d\|_{1}\ \mbox{subject to $\|Az+d-y\|_{2}\leq\epsilon$}.$$

$$\Vert x-\hat{x}{\Vert }_1+\lambda \Vert c-\hat{c}{\Vert }_1\le C_1\left({\sigma }_s(x{)}_1+\lambda {\sigma }_k(c{)}_1\right)+C_2\sqrt{s+{\lambda }^{2}k}ϵ,$$

$$\Vert x-\hat{x}{\Vert }_2+\Vert c-\hat{c}{\Vert }_2\le C_3\left(1+{\eta }^{1/4}\right)\left(\frac{{\sigma }_s(x{)}_1}{\sqrt{s}}+\frac{{\sigma }_k(c{)}_1}{\sqrt{k}}\right)+C_4\left(1+{\eta }^{1/4}\right)ϵ,$$

$$\eta ={\eta }_{s,k}(\lambda )=\frac{s+{\lambda }^{2}k}{\min \{s,{\lambda }^{2}k\}}.$$

$$\Vert x{\Vert }_1+\lambda \Vert c{\Vert }_1+\Vert v_{S^c}{\Vert }_1+\lambda \Vert e_{T^c}{\Vert }_1$$

$$\le 2\Vert x_{S^c}{\Vert }_1+\Vert v_S{\Vert }_1+\Vert \hat{x}{\Vert }_1+\lambda \left(2\Vert c_{T^c}{\Vert }_1+\Vert e_T{\Vert }_1+\Vert \hat{c}{\Vert }_1\right).$$

$$\Vert v_{S^c}{\Vert }_1+\lambda \Vert e_{T^c}{\Vert }_1\le$$

$$+\left(\Vert \hat{x}{\Vert }_1+\lambda \Vert \hat{c}{\Vert }_1\right)-\left(\Vert x{\Vert }_1+\lambda \Vert c{\Vert }_1\right)$$

$$\le 2\left(\Vert x_{S^c}{\Vert }_1+\lambda \Vert c_{T^c}{\Vert }_1\right)+\left(\Vert v_S{\Vert }_1+\lambda \Vert e_T{\Vert }_1\right),$$

$$\Vert v_{S^c}{\Vert }_1+\lambda \Vert e_{T^c}{\Vert }_1\le \frac{2}{1-\rho }\left(\Vert x_{S^c}{\Vert }_1+\lambda \Vert c_{T^c}{\Vert }_1\right)+\frac{\tau \sqrt{s+{\lambda }^{2}k}}{1-\rho }\Vert Av+e{\Vert }_2,$$

$$\Vert Av+e{\Vert }_2\le \Vert A\hat{x}+\hat{c}-y{\Vert }_2+\Vert Ax+c-y{\Vert }_2\le 2ϵ,$$

$$\Vert v_{S^c}{\Vert }_1+\lambda \Vert e_{T^c}{\Vert }_1\le \frac{2}{1-\rho }\left({\sigma }_s(x{)}_1+\lambda {\sigma }_k(c{)}_1\right)+\frac{2\tau }{1-\rho }\sqrt{s+{\lambda }^{2}k}ϵ.$$

$$\Vert v{\Vert }_1+\lambda \Vert e{\Vert }_1$$

$$\le (1+\rho )\left(\Vert v_{S^c}{\Vert }_1+\lambda \Vert e_{T^c}{\Vert }_1\right)+\tau \sqrt{s+{\lambda }^{2}k}\Vert Av+e{\Vert }_2$$

$$\le 2\frac{1+\rho }{1-\rho }\left({\sigma }_s(x{)}_1+\lambda {\sigma }_k(c{)}_1\right)+\frac{4}{1-\rho }\tau \sqrt{s+{\lambda }^{2}k}ϵ.$$

$${\theta }_v=\underset{i\in S}{\min }|v_i|,{\theta }_e=\underset{j\in T}{\min }|e_j|,\theta =\max \{{\theta }_v,{\theta }_e/\lambda \}.$$

$$\Vert v_{S^c}{\Vert }_2^2+\Vert e_{T^c}{\Vert }_2^2=\sum _{i\notin S}|v_i{|}^2+\sum _{j\notin T}|e_j{|}^2\le {\theta }_v\sum _{i\notin S}|v_i|+{\theta }_e\sum _{j\notin T}|e_j|\le \theta \left(\Vert v_{S^c}{\Vert }_1+\lambda \Vert e_{T^c}{\Vert }_1\right).$$

Full text

In this paper we present a novel compressive sampling-based theoretical analysis for a regularized $\ell^{1}$ minimization algorithm that aims to recover sparse expansion coefficients in the presence of measurement corruptions. Our recovery results are uniform (the theoretical guarantees hold for all compressible signals and compressible corruptions vectors), and prescribe algorithmic regularization parameters in terms of a user-defined a priori estimate on the ratio of measurements that are believed to be corrupted. We also propose an iteratively reweighted optimization algorithm that automatically refines the value of the regularization parameter, and empirically produces superior results. Our numerical results test our framework on several medium-to-high dimensional examples of solutions to parameterized differential equations, and demonstrate the effectiveness of our approach.

1 Introduction

The approximation of function values using point evaluations or samples is necessary in a wide number of applications. Much attention has been focused recently on the approximation technique of compressive sampling (CS): The ability to recover sparse linear representations of a function from a given dictionary. This is a particularly important problem in parametric uncertainty quantification (UQ) where the number of parameters translates into the number of variables on which an unknown function depends (the “dimension” of the problem). It is common for dimension to be very large, and the number of degrees of freedom in classical approximation strategies generally grows exponentially with the dimension. This makes classical computational procedures for approximating functions infeasible for large dimensions.

In contrast, compressive sampling seeks a sparse representation of a function using only a small number of samples or measurements, regardless of the parametric dimension. In a non-intrusive UQ pipeline, each function sample corresponds to a potentially large-scale simulation, and so minimizing the requisite number of samples is desirable. When functions are sparse or compressible in a given basis or dictionary, this reconstruction procedure has the potential to mitigate the exponentially debilitating curse of dimensionality. Algorithms in UQ that utilize compressive sampling have enjoyed great success in recent years [43, 42, 44, 27, 18, 10, 17, 15, 22]. For related theoretical contributions, see [1, 2, 8, 16, 29, 28, 41].

Missing from the sparse recovery UQ contributions above is a concrete strategy for fault-tolerant or resilient algorithms. Ensuring modeling resilience for UQ in the presence of system failures is essential for credible prediction on new and emerging massively parallel systems. Fault-tolerant algorithms in general have become necessary in computational science since node failures on distributed architectures can yield corrupted data (the frequency of which increases as the number of processors increases), or algorithmic run-time software failures can result in polluted simulation results. These failures can generate polluted measurements in unpredictable and sometimes undetectable ways [6].

Faults can occur due to complex combination of internal and external conditions that are difficult to reproduce. For example, bits may suffer random corruption, or physical defects in hardware may cause data faults. Corruption errors during model simulation can be grouped into two main types, soft and hard. In this paper, we consider hard faults as errors that cause the simulation to terminate prematurely and/or return obvious, automatically detectable error values such as NaN or Inf. Hard faults by this definition are easy to identify and mark for discard, thus obviating or ameliorating the need for fault-tolerant algorithms.

In contrast, soft failures are essentially random systematic corruption of results that are not easily identifiable. These soft failures pose challenges in UQ: A soft failure will not cause obvious failure in fault-intolerant UQ methods; however, incorrect model values caused by soft failures can significantly degrade an approximation. It is in this case that we require the development of robust and resilient algorithms that can, ideally, deliver constant levels of performance when faced with a few highly corrupted data points.

To address this issue, fault-tolerant algorithms for UQ have been investigated in the context of multilevel Monte Carlo algorithms [24, 25, 26], and in overdetermined least-squares polynomial recovery problems [31]. To the best of our knowledge, there is no comprehensive research in the UQ literature on fault-tolerant sparse recovery algorithms, and in the compressive sampling literature only a handful of papers [19, 21, 23, 32, 33, 34, 35, 39] deal with the problem of corrupted measurements.

The operative distinction in the problem we consider in this paper is a hardware or software fault resulting in occasional large-magnitude errors; we call this the problem of corruptions. Existing CS algorithms are known to be stable with respect to small noise perturbations, but cannot handle sparse corruptions, i.e., situations when a small number of samples are highly corrupted with the corruption magnitudes much larger than typical noise. In this paper we present novel theory and application studies of a sparse corruptions algorithm for CS. The algorithm we use was considered in [21], but we present more general theoretical guarantees on recovery, including practical guidance for the choice of algorithmic regularization parameters.

For fault-tolerance in the context of the sparse recovery problem, the recovery properties of an ideal resilient algorithm would be agnostic to large-magnitude corruptions in a small number of function samples. As described above, these corruptions can arise due to unknown failure modes in computational models or because of large but intermittent measurement errors. Development of mathematical theory for the corrupted compressive sampling problem, and investigation of a corresponding resilient algorithm for sparse recovery of expansion coefficients are the central goals of this paper. The target applications we investigate are exemplars of a common task in UQ: recovery of approximately sparse expansion coefficients in an orthogonal polynomial (polynomial chaos) basis.

The theory and algorithms developed in this paper have the following features:

• •

The compressive sampling recovery theorems are uniform with respect to the function and the corruptions. That is, the recovery guarantees hold over all compressible functions having sparsely corrupted measurements for a single random sampling of measurements.

• •

The algorithm involves a tunable regularization parameter $\lambda$, and a theoretically optimal choice of this parameter is explicitly determined by our analytical results. This theoretically optimal value is defined only by the ratio of measurement corruptions to signal sparsity. Since signal sparsity is frequently comparable to the number of measurements, this optimal $\lambda$ loosely translates into the fraction of measurement samples that are corrupted. From a user’s point of view, our analysis thus suggests a value of $\lambda$ having knowledge only of the ratio of measurements believed to be corrupted.

• •

In experiments, we observe that optimal values of the regularization parameter are non-trivially dependent on the number of measurements, the signal sparsity, and the number of corruptions. We thus propose an iteratively reweighted algorithm for recovery that learns values of the regularization parameter. Our experiments suggest that these learned algorithmic parameters perform better than the value defined by our theoretical results, and thus this reweighted algorithm is more useful in practice.

• •

The location and magnitude of the corruptions amongst the collection of function samples can be unknown, but the algorithm recovers those locations and the corresponding corruption values.

• •

The algorithm is robust to small, but non-sparse measurement errors – e.g. due to noise, truncation of an infinite polynomial expansion or numerical error in computing function samples – and moreover is noise-blind. That is to say, it requires no a priori upper bound on such errors.

• •

The optimization problem we solve to compute solutions is from [21], but our work is both a theoretical and practical advancement over the results in that reference. In order to show the solution computed is indeed the original sparse solution, [21] uses conditions on the restricted isometry constant (RIC) of the measurement matrix. Our results are a significant relaxation of previously reported conditions on the RIC (compare conditions on $\delta_{2s,2k}$ in Lemma 2.3 of [21] versus our Theorem 3.7, equation (3.8), and the discussion in Section 3.4). The results for general sensing matrices in [21] are nonuniform with respect to the signal and corruptions support, and require certain models for the signal and corruptions; our results are uniform and require no model for the signal or corruptions, other than compressibility. Finally, our paper is also devoted to numerical investigation of the performance of the method, including practical guidance for choosing the regularization parameter $\lambda$; such thorough investigations are absent in [21].

We first introduce notation and summarize the main mathematical statements of this paper in Section 2. This is followed in Section 3 by our theoretical analysis. Section 4 presents numerical results to complement our theoretical analysis and verify the practical efficacy of the algorithm.

2 Model problem and main results

Let $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ denote an unknown function, and let $\left\{\phi_{j}\right\}_{j=1}^{N}$ be a given dictionary of functions, $\phi_{j}:\mathbb{R}^{d}\rightarrow\mathbb{R}$. For example, the functions $\phi_{j}$ are frequently multivariate polynomial chaos basis elements; our capstone numerical examples will show results from such a basis. In scenarios of interest, the size $N$ of the dictionary is very large.

The ultimate goal is to recover coefficients $x_{j}$ that determine the approximation

[TABLE]

using samples of $f$, where $n(\xi)$ is an assumed small discrepancy term between the exact function and its $N$-term linear approximation in $\phi_{j}$444Our notation suggests that $n=n(\xi)$ depends explicitly and deterministically on $\xi$; however, our theory encompasses the case when $n$ is a stochastic variable or process, e.g., independent Gaussian random variable additive perturbations of the measurements.. For the purposes of exposition we assume $|n(\xi)|\leq\epsilon$ for some known uniform noise bound $\epsilon$; we will show later that lack of a priori knowledge for this bound only affects theoretical results in benign ways. As described above, we assume the vector $x=\left(x_{1},\ldots,x_{N}\right)^{T}\in\mathbb{R}^{N}$ to be compressible. Sparsity or compressibility of a vector can be quantified via its best $s$-term approximation error,

[TABLE]

where $\|\cdot\|_{p}$ is the standard $\ell^{p}$ norm on vectors; for $p=0$, $\|x\|_{0}$ is the sparsity of $x$, i.e., the number of non-zero elements in the vector.

With $\left\{\xi_{1},\ldots,\xi_{m}\right\}\subset\mathbb{R}^{d}$ a collection of samples of $\xi$, we have the corresponding corrupted function measurements,

[TABLE]

where the corruption vector $c=\left(c_{1},\ldots c_{m}\right)^{T}\in\mathbb{R}^{m}$ is assumed to be $k$-sparse but can have large entries. To enforce an underdetermined system, we assume $m<N$. Defining the rectangular matrix $A$ with entries $(A)_{j,k}=\phi_{k}\left(\xi_{j}\right)$, then the unknown vectors $x$ and $c$ satisfy the underdetermined linear system

[TABLE]

In order to compute the solution $(x,c)$ having knowledge of only $A$ and $y$, we consider the following model problem (see also [21] and references therein):

[TABLE]

Let $(\hat{x},\hat{c})$ be a minimizer of this problem, where $\hat{x}\in\mathbb{R}^{N}$ and $\hat{c}\in\mathbb{R}^{m}$. Our objective is to obtain conditions on $A$ (in particular, on the number of measurements $m$) and $\lambda$ such that the error

[TABLE]

can be bounded by the best approximation numbers $\sigma_{s}(x)_{1}$ and $\sigma_{k}(c)_{1}$, and the noise magnitude $\epsilon$.

2.1 Main results

In all that follows, the statement $a\lesssim b$ means $a\leq Cb$ for some universal constant $C$. Our first main result shows that stable and robust recovery of $x$ and $c$ is implied by a certain modification of the classical Restricted Isometry Property (RIP) which incorporates the sparse corruptions term (Definition 3.5). Specifically, Theorem 3.7 establishes that if the matrix $A$ satisfies the RIP for the corruptions problem of order $(2s,2k)$ (see Definition 3.5) with constant $\delta_{2s,2k}$ satisfying

[TABLE]

then the following error bounds hold:

[TABLE]

Our second main result (Theorem 3.15) provides explicit conditions on $m$, $s$ and $k$ for (2.4) to hold for matrices of so-called bounded orthonormal system [12, Chpt. 12]. Specifically, suppose that $\{\phi_{j}\}^{N}_{j=1}$ is an $L^{2}_{\mathrm{d}\nu}(D)$-orthonormal system, where $\nu$ is a probability measure and $D\subset\mathbb{R}^{d}$ its support. Define

[TABLE]

and let $A=\left\{\phi_{j}(\xi_{i})\right\}^{m,N}_{i,j=1}$ where $\xi_{1},\ldots,\xi_{m}$ are drawn i.i.d. according to $\nu$. If

[TABLE]

then with probability at least $1-\epsilon$, the restricted isometry constant $\delta_{2s,2k}$ of the scaled matrix $\frac{1}{\sqrt{m}}A$ satisfies $\delta_{2s,2k}\leq\delta$.

One can see from these estimates that optimizing $\eta$ over values of $\lambda$ yields a minimum value of $\eta=2$ when $\lambda^{2}=s/k$. Assuming $s\sim m$, this provides a concrete determination of the parameter $\lambda$ for use in (2.3) having knowledge only of the ratio of corrupted measurements. We note in passing that we do not believe that the second condition in (2.6) is sharp in the dependence on the product $s\cdot k$. Improvement of this to a condition of the form

[TABLE]

is left as a topic for future work. Note that such a condition is known for Gaussian random matrices. Moreover, a nonuniform recovery result with the scaling (2.7) for exactly sparse coefficients $x$ and corruptions $c$ having random sign patterns was given in [21]. See Section 3.3 for further discussion.

It is common in compressed sensing to assume some a priori known noise bound $\epsilon$ based on the user’s knowledge of measurement noise or truncation error. Although there are some results that circumvent this assumption [1, 2], they typically yield somewhat weaker recovery guarantees. However, in the context of the sparse corruptions theory presented above, such prior knowledge of $\epsilon$ is not necessary for stable recovery: The error introduced by an unknown noise $\epsilon$ can be passed into theoretical estimates as a penalty of size $\epsilon$. To see this, note that if we define $c^{\prime}\coloneqq\frac{1}{\sqrt{m}}(c+n)$, then the system $y=Ax+c+n$ can be written as $\frac{1}{\sqrt{m}}y=\frac{1}{\sqrt{m}}Ax+c^{\prime}$. Solving (2.3) by setting $\epsilon=0$ results in the $\epsilon=0$ version of the estimate (2.5b) with $c^{\prime}$ replacing $c$. However, the normalized best $k$-term approximation error to $c^{\prime}$ appearing in (2.5b) is stable with respect to noise perturbations:

[TABLE]

Here $\epsilon\geq\|n\|_{\infty}$ is any bound for the perturbation $n$ in the uniform norm. Using (2.7), we see that $\sqrt{\frac{m}{k}}\epsilon\lesssim\epsilon$, which is on the same order as the estimate (2.5b) that uses a priori knowledge of $\epsilon$. A similar argument holds for the bound (2.5a).

While our theoretical results are thus insensitive to ignorance about small noise levels, we caution that it is always a good idea to use such information in practical recovery algorithms if available, e.g. as the result of cross validation. See, for example, [10, 43, 17].

2.2 Remarks on numerical results

We postpone presenting numerical results until the end of this paper in Section 4. However, some remarks on our findings are pertinent here in the context of the previous section’s theory. First, the optimal value of $\lambda^{2}=s/k$ that is suggested by (2.4) does not appear to be the computationally optimal value of $\lambda$. That this fixed value of $\lambda$ is not the best is not surprising since the bounds (2.5) are derived using some loose inequalities. However, such bounds can be useful in understanding qualitative trends. Results from our experimentation do suggest that large values of $\lambda$ more reliably recover corruptions when $s/k$ is large (see Figures 1 and 2). This general trend in numerical results is consistent with the behavior of $\eta$ in (2.4) as a function of $\lambda$ when $s/k$ is large.

We address this discrepancy between the theory and empirical results by propose an iteratively reweighted $\ell^{1}$ optimization scheme (see [7]) that learns and updates the value of $\lambda$. Our results show that this proposed algorithm performs much better in practice than algorithms that fix $\lambda$. However, we do not present any theory to support the observed superiority of reweighted $\ell^{1}$ optimization schemes for the corruptions problem.

Many of our capstone numerical examples are from applications using polynomial chaos expansions, where the compressible function has an expansion in a multivariate orthogonal polynomial basis. To simplify the presentation of our results, we focus on such examples where the basis is a tensor-product Legendre polynomial or Chebyshev polynomial system. Much recent work has shown that randomly generating measurements using samples from standard distributions (e.g., the uniform distribution) can accurately and near-optimally recover orthogonal polynomial expansions from such basis sets [28, 17, 41]. Recovery in more general polynomial spaces has been investigated [16, 18, 15], but these methods usually rely on sophisticated sampling strategies and optimal sampling schemes are still an active area of research.

3 Theory for the sparse corruptions problem

We recall and summarize our notation for the sparse corruptions problem in Table 1. Our previous discussion was framed for real-valued signals $x$ and measurements $y$, but we now generalize to the complex-valued setting. This adds generality with no additional mathematical difficulty.

We follow a familiar path for deriving conditions on $m$ such that $\ell^{1}$ optimization problems recover sparse solutions (see, for example, [12]). Section 3.1 defines an appropriate robust Null Space Property (NSP) for the matrix $A$ in the sparse corruptions setting. Under this property, we show that the recovery estimates (2.5) hold. In order to construct matrices $A$ that satisfy the robust NSP, Section 3.2 generalizes the concept of the Restricted Isometry Property (RIP) for matrices to the sparse corruptions setting. That section shows that matrices satisfying the RIP for the sparse corruptions problem also satisfy the robust NSP. Sections 3.3 and 3.3.2 show that if the dictionary elements $\phi_{j}$ form a bounded orthonormal system, then under the condition (2.6), the matrix $A$ satisfies the RIP with high probability. Finally, using these various results, we discuss a theoretically-optimal choice for $\lambda$ in Section 3.4.

3.1 The Robust Null Space Property for the sparse corruptions problem

The following two definitions are generalizations of robust null space properties (cf. [12, Definition 4.17] and [12, Definition 4.21], respectively), and prescribe classes of matrices whose kernels do not contain sparse vectors.

Definition 3.1**.**

Let $1\leq s\leq N$, $1\leq k\leq m$ and $\lambda>0$. A matrix $A\in\mathbb{C}^{m\times N}$ satisfies the $\ell^{1}$-robust null space property of order $(s,k)$ with weight $\lambda$ if there exist constants $0<\rho<1$ and $\tau>0$ such that

[TABLE]

for all sets $S\subseteq\{1,\ldots,N\}$ and $T\subseteq\{1,\ldots,m\}$ with $|S|\leq s$ and $|T|\leq k$. Above, $S^{c}$ is the complement of $S$ in $\{1,\ldots,N\}$, and similarly for $T^{c}$.

Definition 3.2**.**

Let $1\leq s\leq N$, $1\leq k\leq m$ and $\lambda>0$. A matrix $A\in\mathbb{C}^{m\times N}$ satisfies the $\ell^{2}$-robust null space property of order $(s,k)$ with weight $\lambda$ if there exist constants $0<\rho<1$ and $\tau>0$ such that

[TABLE]

for all sets $S\subseteq\{1,\ldots,N\}$ and $T\subseteq\{1,\ldots,m\}$ with $|S|\leq s$ and $|T|\leq k$.

These definitions yield the following two results:

Lemma 3.3**.**

If $A\in\mathbb{C}^{m\times N}$ satisfies the $\ell^{2}$-robust null space property of order $(s,k)$ with weight $\lambda>0$ and constants $0<\rho<1$, $\tau>0$ then it satisfies the $\ell^{1}$-robust null space property of order $(s,k)$ with weight $\lambda>0$ and constants $\rho$, $\tau\sqrt{s+\lambda^{2}k}$.

Proof.

Observe that

[TABLE]

We now use the definition of the $\ell^{2}$-robust null space property. ∎

Theorem 3.4**.**

Let $1\leq s\leq N$, $1\leq k\leq m$ and $\lambda>0$ and suppose that $A\in\mathbb{C}^{m\times N}$ satisfies the $\ell^{2}$-robust null space property of order $(s,k)$ with weight $\lambda$. Let $x\in\mathbb{C}^{N}$, $c\in\mathbb{C}^{m}$, $y\in\mathbb{C}^{m}$ and $\epsilon>0$ be such that $\|Ax+c-y\|_{2}\leq\epsilon$, and suppose that $(\hat{x},\hat{c})$ is a minimizer of

[TABLE]

*Then *

[TABLE]

and

[TABLE]

where the constants $C_{1},C_{2},C_{3},C_{4}$ depend on $\rho$ and $\tau$ only and $\eta$ is given by

[TABLE]

Proof.

We first prove (3.2). Lemma 3.3 implies that $A$ satisfies the $\ell^{1}$-robust null space property. Let $S\subseteq\{1,\ldots,N\}$, $|S|\leq s$ and $T\subseteq\{1,\ldots,m\}$, $|T|\leq k$ be such that $\|x_{S^{c}}\|_{1}=\sigma_{s}(x)_{1}$ and $\|c_{T^{c}}\|_{1}=\sigma_{k}(c)_{1}$. Then, if $v=x-\hat{x}$ and $e=c-\hat{c}$ we have

[TABLE]

Rearranging now gives

[TABLE]

where in the second inequality we note that $\|x\|_{1}+\lambda\|c\|_{1}\geq\|\hat{x}\|_{1}+\lambda\|\hat{c}\|_{1}$ since $(x,c)$ is feasible and $(\hat{x},\hat{c})$ is a minimizer. The $\ell^{1}$-robust null space property now implies that

[TABLE]

and since $\|x_{S^{c}}\|_{1}=\sigma_{s}(x)_{1}$, $\|c_{T^{c}}\|_{1}=\sigma_{k}(c)_{1}$ and

[TABLE]

we deduce that

[TABLE]

Finally, to complete the proof of (3.2) we argue as follows:

[TABLE]

Here, we use the $\ell^{1}$-robust null space property in the second step, and (3.5) and (3.6) in the third step.

We now consider (3.3). Writing $v=x-\hat{x}$ and $e=c-\hat{c}$ as before, let $S$ be the index of the largest $s$ elements of $v$ in absolute value and $T$ be the index set of the largest $k$ elements of $e$ in absolute value. Define

[TABLE]

Then

[TABLE]

Now observe that $\theta_{v}\leq\|v_{S}\|_{2}/\sqrt{s}$ and $\theta_{e}\leq\|e_{T}\|_{2}/\sqrt{k}$, and therefore

[TABLE]

where in the second step we use the $\ell^{2}$-robust null space property and (3.5). Combining this with the previous estimate and using the definition of $\eta$ gives

[TABLE]

where we have defined the non-negative scalar $w$ as

[TABLE]

Completing the square with respect to $w$ under the brackets yields

[TABLE]

Using the $\ell^{2}$-robust NSP on the pair $(v,e)$ along with the above estimate, we have

[TABLE]

We note that

[TABLE]

Combining the above with (3.7) proves (3.3). ∎

3.2 The Restricted Isometry Property for the sparse corruptions problem

The robust NSP is typically difficult to prove directly. Hence we now introduce the Restricted Isometry Property (RIP) for the sparse corruptions problem, and show that it implies the robust NSP. Note that this has been defined previously in [21, Defn. 2.1].

Definition 3.5**.**

Let $1\leq s\leq N$, $1\leq k\leq m\leq N$ and $A\in\mathbb{C}^{m\times N}$. The $(s,k)^{\mathrm{th}}$ Restricted Isometry Constant (RIC) $\delta=\delta_{s,k}$ of the matrix $A$ is the smallest constant such that

[TABLE]

for all $x\in\Sigma_{s}$ and $c\in\Sigma_{k}$. If $0<\delta_{s,k}<1$ then we say that $A$ has the Restricted Isometry Property (RIP) of order $(s,k)$.

Our first result is the following:

Lemma 3.6**.**

Let $1\leq s\leq N$, $1\leq k\leq m\leq N$, $\lambda>0$ and $A\in\mathbb{C}^{m\times N}$. If $A$ satisfies the RIP of order $(2s,2k)$ with constant

[TABLE]

where $\eta$ is as in (3.4), then $A$ satisfies the $\ell^{2}$-robust NSP of order $(s,k)$ with weight $\lambda$ and constants $0<\rho<1$ and $\tau>0$ depending only on $\delta_{2s,2k}$.

The proof of this result is given next. Combining this lemma with Theorem 3.4 now yields our main result:

Theorem 3.7**.**

Let $1\leq s\leq N$, $1\leq k\leq m$ and $\lambda>0$ and suppose that $A\in\mathbb{C}^{m\times N}$ satisfies the RIP of order $(2s,2k)$ with constant $\delta_{2s,2k}$ satisfying (3.8) and $\eta$ as in (3.4). Let $x\in\mathbb{C}^{N}$, $c\in\mathbb{C}^{m}$, $y\in\mathbb{C}^{m}$ and $\epsilon>0$ be such that $\|Ax+c-y\|_{2}\leq\epsilon$, and suppose that $(\hat{x},\hat{c})$ is a minimizer of

[TABLE]

Then

[TABLE]

*where the constants $C_{1},C_{2},C_{3},C_{4}$ depend on $\delta_{2s,2k}$ only. *

We now prove Lemma 3.6. We first require the following:

Lemma 3.8**.**

Let $1\leq s\leq N$, $1\leq k\leq m\leq N$, and let $A\in\mathbb{C}^{m\times N}$ satisfy the RIP of order $(2s,2k)$ with constant $\delta_{2s,2k}$. Suppose that $x\in\Sigma_{s}$ and $c\in\Sigma_{k}$ are such that

[TABLE]

for some $t$ with $0\leq|t|\leq\delta_{2s,2k}$. If $z\in\Sigma_{s}$ and $d\in\Sigma_{k}$ are orthogonal to $x$ and $c$, respectively, then

[TABLE]

Proof.

Assume that $\|x\|^{2}_{2}+\|c\|^{2}_{2}=\|z\|^{2}_{2}+\|d\|^{2}_{2}=1$ without loss of generality. Let $\alpha,\beta\in\mathbb{R}$ and $\gamma\in\mathbb{C}$ and notice that $\alpha x+\gamma z,\beta x-\gamma z\in\Sigma_{2s}$ and $\alpha c+\gamma d,\beta c-\gamma d\in\Sigma_{2k}$. Therefore

[TABLE]

Note that in the second step we use orthogonality of the vectors $x$ and $z$ and $c$ and $d$. Similarly,

[TABLE]

Subtracting the second equation from the first gives

[TABLE]

On the other hand

[TABLE]

Combining this with (3.9) gives

[TABLE]

Now let $\gamma$ be such that $|\gamma|=1$ and $\mathrm{Re}\,\left(\bar{\gamma}\langle Ax+c,Az+d\rangle\right)=|\langle Ax+c,Az+d\rangle|$. Then, after rearranging, we get

[TABLE]

We now seek values $\alpha$ and $\beta$ which minimize the right-hand side of this expression. If $t=\delta_{2s,2k}$ then the minimal value [math] is attained by setting $\beta=0$ and letting $\alpha\rightarrow\infty$. Conversely, if $t<\delta_{2s,2k}$ the minimal value is attained when $\alpha=\sqrt{\frac{\delta_{2s,2k}+t}{\delta_{2s,2k}-t}}$ and $\beta=\frac{1}{\alpha}$. This gives

[TABLE]

which completes the proof. ∎

Proof of Lemma 3.6.

Let $x\in\mathbb{C}^{N}$ and $c\in\mathbb{C}^{m}$. To prove the $\ell^{2}$-robust NSP for $A$ it is enough to show that (3.1) holds when $S=S_{0}$ is the index set of the $s$ largest coefficients of $x$ in absolute value and $T=T_{0}$ is the set of the $k$ largest values of $c$ in absolute value. Given $S_{0}$, let $S_{1}$ be the index set of the next $s$ largest coefficients of $x$ in absolute value, $S_{2}$ be the index set of the next $s$ largest coefficients and so on. Define $T_{1},T_{2},\ldots$ in a similar way. We now have the following:

[TABLE]

Let $0\leq|t|\leq\delta_{2s,2k}$ be such that

[TABLE]

and note that this gives

[TABLE]

For the second term of (3.10), we use the disjointness of $S_{0}$ and $S_{j}$ and $T_{0}$ and $T_{j}$ for $j\geq 1$ in combination with Lemma 3.8 to get

[TABLE]

Let $x^{+}_{j}$ and $x^{-}_{j}$ be the largest entries of $x_{S_{j}}$ in absolute value. Then, by [12, Lem. 6.14], we have

[TABLE]

Similarly,

[TABLE]

which gives

[TABLE]

Therefore, combining this with (3.10), (3.11), (3.12) and (3.13) yields

[TABLE]

Consider the function $g(t)=\frac{\delta^{2}_{2s,2k}-t^{2}}{(1+t)^{2}}$, where $0\leq t\leq\delta_{2s,2k}$. This function attains its maximum value at $t=-\delta^{2}_{2s,2k}$ and takes value $\frac{\delta^{2}_{2s,2k}}{1-\delta^{2}_{2s,2k}}$ there. Additionally $\frac{1}{\sqrt{1+t}}\leq\frac{1}{\sqrt{1-\delta_{2s,2k}}}$. Hence we get

[TABLE]

After noting that $\|x_{S_{0}}\|_{2}+{\left\|c_{T_{0}}\right\|}_{2}\leq\sqrt{2}\sqrt{{\left\|x_{S_{0}}\right\|}^{2}_{2}+{\left\|c_{T_{0}}\right\|}^{2}_{2}}$ and rearranging, we obtain

[TABLE]

where

[TABLE]

To complete the proof we note that $\rho,\tau>0$ provided $\delta_{2s,2k}<\sqrt{8/9}$. This holds by assumption, since $\eta\geq 2$ and therefore the condition (3.8) implies that $\delta_{2s,2k}<\sqrt{8/33}<\sqrt{8/9}$. Also, after rearranging we see that $\rho<1$ if

[TABLE]

which again holds by assumption. ∎

• Remark 3.9

The RIP for the sparse corruptions problem is a special case of the RIP in levels (RIPL), introduced in [5]. The RIPL applies to vectors that are sparse in levels; namely, having different amounts of sparsity in different (but fixed) sections of the vector. In the context of the sparse corruptions problem, this corresponds to the concatenated vector $z=[x;c]$, which is $s$-sparse in its first $N$ entries and $k$-sparse in its remaining $m$ entries. As a general tool, sparsity in levels has been used in the context of compressive imaging [3, 4, 30], radar [11] and multi-sensor acquisition [9]. It is interesting that the same model also occurs naturally in the, seemingly unrelated, sparse corruptions problem. We note in passing that Theorems 3.4 and 3.7 follow a similar approach to that of [5] with some changes made to incorporate the weighted optimization problem.

3.3 Matrices that satisfy the RIP for sparse corruptions

We first recall the classical RIP for sparse vectors:

Definition 3.10**.**

Let $1\leq s\leq N$ and $A\in\mathbb{C}^{m\times N}$. The $s^{\mathrm{th}}$ Restricted Isometry Constant (RIC) $\delta=\delta_{s}$ of the matrix $A$ is the smallest constant such that

[TABLE]

for all $x\in\Sigma_{s}$. If $0<\delta_{s}<1$ then we say that $A$ has the Restricted Isometry Property (RIP) of order $s$.

To distinguish it from the RIP for the sparse corruptions problem (Definition 3.5), we shall refer to this as the RIP for sparse vectors.

Lemma 3.11**.**

Let $1\leq s\leq N$, $1\leq k\leq m$, $A\in\mathbb{C}^{m\times N}$ and define

[TABLE]

where $A_{S,T}\in\mathbb{C}^{|T|\times|S|}$ is the submatrix of $A$ with entries $\{A_{ij}\}_{i\in T,j\in S}$. Suppose that $A$ has the RIP for sparse vectors with constant $\delta_{s}$ and that $\sigma_{s,k}<\sqrt{1-\delta_{s}}$. Then $A$ has the RIP of order $(s,k)$ for the sparse corruptions problem with constant

[TABLE]

In other words,

[TABLE]

for all $x\in\Sigma_{s}$ and $c\in\Sigma_{k}$.

Proof.

Let $x\in\Sigma_{s}$ and $c\in\Sigma_{k}$ and write $S=\mathrm{supp}(x)$ and $T=\mathrm{supp}(c)$. Then

[TABLE]

By Young’s inequality

[TABLE]

for any $\epsilon>0$. Hence

[TABLE]

Solving the equation $\delta_{s}+\sigma_{s,k}/\epsilon=\sigma_{s,k}\epsilon$ yields the value $\epsilon=\frac{\delta_{s}+\sqrt{\delta^{2}_{s}+4\sigma^{2}_{s,k}}}{2\sigma}$, and substituting this value of $\epsilon$ into the previous expression yields the proof. ∎

This result shows that any matrix satisfying the RIP for sparse vectors also satisfies the RIP for the sparse corruptions problem, provided the all $k\times s$ submatrices have small spectral norm.

3.3.1 Gaussian random matrices

Gaussian random matrices in the context of the sparse corruptions problem were considered in [21]. The following result essentially recaps the main result for this case given therein. We include a short proof for completeness:

Theorem 3.12**.**

Let $0<\delta,\epsilon<1$, $1\leq s\leq$, $1\leq k\leq m$ and suppose that

[TABLE]

Let $A\in\mathbb{C}^{m\times N}$ be a matrix whose entries are independent Gaussian random variables with mean zero and variance $1$. Then with probability at least $1-\epsilon$, the matrix $\frac{1}{\sqrt{m}}A$ has the RIP for the sparse corruptions problem of order $(s,k)$ with constant $\delta_{s,k}\leq\delta$.

Proof.

Lemma 3.11 asserts that $A$ has the RIP of order $(s,k)$ for the sparse corruptions problem with constant $\delta_{s,k}\leq\delta$ provided (i) $A$ has the RIP of order $s$ with $\delta_{s}\leq\delta/\sqrt{2}$ and (ii) the constant $\sigma_{s,k}$ defined in (3.15) satisfies $\sigma_{s,k}\leq\delta/(2\sqrt{2})$. Hence, by the union bound it suffices to show that (3.16) and (3.17) imply both (i) and (ii) separately with probabilities at least $1-\epsilon/2$. Due to a standard result in compressed sensing (see, for example, [12, Thm. 9.2]), property (i) holds with probability at least $1-\epsilon/2$ whenever the condition (3.16) is satisfied. We now consider property (ii). First, notice that $\sigma_{s,k}$ is increasing in $k$. Therefore, we may assume that $k\asymp\delta^{2}\cdot m$, i.e. $k\gtrsim\delta^{2}\cdot m$ and $k\lesssim\delta^{2}\cdot m$. Fix subsets $S\subseteq\{1,\ldots,N\}$ and $T\subseteq\{1,\ldots,m\}$ with $|S|=s$ and $|T|=k$. Then, due to a known result for singular values of random Gaussian matrices (see, for example, [38, Cor. 5.35]), we have

[TABLE]

The conditions (3.16) and (3.17) imply that $\sqrt{s/m}\leq\delta/(6\sqrt{2})$ and $\sqrt{k/m}\leq\delta/(6\sqrt{2})$. Hence, by the union bound

[TABLE]

In particular, $\mathbb{P}\left(\sigma_{s,k}>\delta/(2\sqrt{2})\right)\leq\epsilon/2$ provided

[TABLE]

Since $k\asymp\delta^{2}\cdot m$, we have $\log(\mathrm{e}m/k)\lesssim\log(2\delta^{-1})$. Hence this condition is implied by (3.16) and (3.17). This establishes property (ii) and completes the proof. ∎

This result asserts that Gaussian random matrices can recover a fixed fraction $k/m\leq c$ of corruptions (see (3.17)) and (up to constants) the same level of sparsity $s$ as in the uncorrupted case (see (3.16)).

3.3.2 Bounded orthonormal systems

Gaussian random matrices, while mathematically appealing, are of little relevance to multivariate approximation using Polynomial Chaos expansions. In this case, a more suitable framework is that of bounded orthonormal systems (see, for example, [12, Chpt. 12]):

Let $D$ be a domain with a probability measure $\nu$ and $\phi_{1},\ldots,\phi_{N}$ be an orthonormal system of complex-value functions in $L^{2}(D)$. Recall that this system is bounded if

[TABLE]

Given such a system, we construct the measurement matrix $A$ as

[TABLE]

where $t_{i}$ are drawn independently at random according to the probability measure $\nu$.

Theorem 3.13**.**

Let $A\in\mathbb{C}^{m\times N}$ be the matrix of a bounded orthonormal system, $1\leq s\leq N$ and $0<\delta,\epsilon<1$. If

[TABLE]

then $A$ satisfies the RIP for sparse vectors with probability at least $1-\epsilon$.

We remark in passing that the logarithmic dependence in $s$ can be improved by one power, at the expense of a larger factor in $\delta^{-1}$ [8]. However, this may not be best for the purposes of this paper, since in view of Theorem 3.7, $\delta^{-2}$ scales linearly in the parameter $\eta$ (see next).

The following lemma estimates the constant $\sigma_{s,k}$ for matrices of the form (3.18):

Lemma 3.14**.**

Let $A\in\mathbb{C}^{m\times N}$ be the matrix of a bounded orthonormal system, $1\leq s,k\leq N$ and $\sigma_{s,k}$ be as in (3.15). Then

[TABLE]

Proof.

Fix subsets $S\subseteq\{1,\ldots,N\}$, $|S|=s$ and $T\subseteq\{1,\ldots,m\}$, $|T|=k$ and let $x\in\mathbb{C}^{N}$ and $c\in\mathbb{C}^{m}$ with $\mathrm{supp}(x)=S$ and $\mathrm{supp}(c)=T$. Then

[TABLE]

Hence $\|P_{T}AP_{S}\|_{2}\leq\sqrt{\frac{K^{2}sk}{m}}$. This now gives the result. ∎

With this in hand, we now deduce the following result:

Theorem 3.15**.**

Let $1\leq s\leq N$, $1\leq k\leq m$, $0<\delta,\epsilon<1$ and suppose that

[TABLE]

and

[TABLE]

Then, with probability at least $1-\epsilon$, $A$ has the RIP of order $(s,k)$ for the sparse corruptions problem with constant $\delta_{s,k}\leq\delta$.

Proof.

Theorem 3.13 and (3.19) imply that $A$ has the RIP of order $s$ with $\delta_{s}\leq\delta/\sqrt{2}$ with probability at least $1-\epsilon$. Moreover, Lemma 3.14 and (3.15) imply that $\sigma_{s,k}\leq\delta/(2\sqrt{2})$. We now apply Lemma 3.11. ∎

• Remark 3.16

This result asserts that the number of corruptions that can be tolerated is a fraction of $m/s$. This is inferior to the case of Gaussian random measurements, where Theorem 3.12 gives that a fraction of $m$ corruptions are permitted. We conjecture, however, that a similar estimate can be proved for the bounded orthonormal systems case – indeed, a nonuniform recovery result of this form was proved in [21] for the case of exactly sparse coefficients $x$ and corruptions $c$ with random sign sequences – albeit with a substantially more sophisticated argument than the proof of Theorem 3.12. In particular, while estimates for the singular values of matrices of bounded orthonormal systems are known [38], they are more stringent than those for Gaussian random matrices. Using these estimates and arguing via the union bound (as in the proof of Theorem 3.12) unfortunately results in an estimate similar to (3.15). We also note in passing that while there exist RIP estimates for quite general matrices under the sparsity in levels model [20] (see Remark 3.2), these unfortunately do not apply to the setup of the sparse corruptions problem. We therefore leave the problem of improving (3.15) for future work.

3.4 Strategy for choosing $\lambda$

Regardless of the matrix $A$, our main theorems (Theorems 3.7 and 3.13) suggest an optimal strategy for choosing the parameter $\lambda$. Notice that the restricted isometry constant $\delta$ enters into the measurement condition in Theorem 3.13 as $\delta^{-2}$. Since Theorem 3.7 requires that (3.8) holds, the measurement condition contains a factor that is at least as large as

[TABLE]

We wish to minimize this factor so as to reduce the measurement condition as much as possible. This can be done by minimizing $\eta$, which in turn yields the theoretically-optimal scaling

[TABLE]

Notice that this gives the value $\eta=2$. In particular, the condition (3.8) becomes

[TABLE]

with right-hand side independent of $s$ and $k$. We remark in passing that the choice (3.20) is implicitly made in [21]. However, the condition given in [21, Lem. 2.3] is $\delta_{2s,2k}<1/18\approx 0.056$ which is significantly more stringent than (3.20). Moreover, [21] only considers exact sparsity, whereas Theorem 3.7 also treats the case of stable recovery of inexactly sparse coefficients and corruptions.

4 Numerical experiments

We divide our numerical results into two main sections. The goal of Section 4.1 is to study the behavior of numerical algorithms in the context of the theoretical estimates presented earlier. In particular, we investigate the influence that the regularization parameter $\lambda$ has on recovery properties. We confine these investigations to problems with manufactured sparsity so that systematic studies may be carried out. The lessons learned from these studies allow us to formulate and propose an iteratively reweighted alternative to the one-time optimization (2.3). Note that none of our theoretical error estimates apply to algorithms with weighted norms. However, weighted $\ell^{1}$ schemes can provide empirically superior results, e.g., [43]. Thus, we explore weighted algorithms because their use is natural from a practical point of view, but is not in the scope of our theoretical analysis. Our simulations in this section use the SPGL1 package [36, 37].

The second collection of results, Section 4.2, focuses on more practical scenarios in scientific computing, dealing with recovery of sparse or compressible polynomial Chaos expansions of solutions to parameterized differential equations. Here we use the algorithmic lessons learned from Section 4.1 to illustrate the efficacy and fault-tolerance of our approaches on realistic problems in the presence of measurement corruptions.

4.1 Recovery of manufactured solutions with sparse corruptions

This section is primarily concerned with the generation of phase recovery diagrams for the sparse corruptions problem. In particular, our tests here are not necessarily motivated by sensing matrices and corruptions from function approximation, but instead are designed to understand behavior of the algorithms. The following standard experiment for accomplishing this is carried out: We fix the number of measurements $m$ and the dictionary size $N$, and we vary the signal sparsity $s$ and the number of measurement corruptions $k$. For each $s$ and $k$ we generate an $s$-sparse signal $x$, and for a given model of a measurement matrix $A$, we generate $m$ measurements $y$ from the signal $x$, and subsequently corrupt (highly pollute) these measurements with a $k$-sparse vector $c$, whose non-zero entries are $CZ$, where $Z$ is a random draw from a certain probability distribution and $C>0$ is a scaling constant. In this test, $Z$ is a standard normal random variable and $C=1$.

We then run the recovery algorithm (2.3) for a given value of $\lambda$, producing a recovered signal $\widehat{x}$ and measurement corruption vector $\widehat{c}$. We define the recovery as successful if $\left\|x-\widehat{x}\right\|^{2}+\left\|c-\widehat{c}\right\|^{2}<\epsilon_{\mathrm{tol}}$. In this test, we set the success tolerance to be $\epsilon_{\mathrm{tol}}=10^{-4}$.

In the test above, the generation of $x$, and of $y$, and of $c$, are statistically independent555Measurement corruptions are generated as iid standard normal random variables, and support indices in a sparse vector are generated using the uniform probability law (draws without replacement) on the index set.. For each $s$ and $k$, the above procedure is run $T\in\mathbb{N}$ times with independent draws, and an empirical estimate of the probability of “success” is computed. In the phase transitions plots below, we use $T=10$ simulations.

The phase transitions color each pixel, corresponding to a particular value of $s$ and $k$, according to the empirical success probability. The phase transition axes are $s/m$ and $k/m$, and thus each ranges in the interval $[0,1]$, but we truncate to $[0,0.5]$ in our plots because this region is sufficient to illustrate behavior. We consider the following two models of measurement matrix $A$:

• $\bullet$  

Model 1: a Gaussian random matrix

• $\bullet$  

Model 2: a randomly-subsampled Discrete Fourier Transform (DFT) matrix

Note that Model 2 is an example of a bounded orthonormal system. We compare several different choices of $\lambda$ for each model.

4.1.1 Phase transition plots for fixed $\lambda$

Figures 1 and 2 display the results for models 1 and 2 described above, respectively. Each figure shows an array of plots; the columns correspond to differing values of $m$, increasing from left to right; the rows correspond to differing values of $\lambda$, increasing from top to bottom, except the last two rows, which show the “optimal” value of $\lambda=\sqrt{s/k}$ suggested by the theory, and the iterative reweighting procedure described in the next section.

Comparing the results for $\lambda=\sqrt{s/k}$ (row 5 in the plots) with the other plots with $\lambda$ fixed, we see that $\lambda=\sqrt{s/k}$ does not behave optimally in practice, even though this is suggested by our theory. Indeed, further experimentation reveals that the behavior of these transition plots changes notably when $m$ is varied. However, the following observations are consistent across all our runs:

• •

When there are few corruptions relative to the signal sparsity ($k\ll s$), larger values of $\lambda$ tend to perform better. This general trend is consistent with the theory from previous sections: Our recovery results are stated in terms of a quantity $\eta$ defined in (3.4), and when $k\ll s$, we require large $\lambda$ to make $\eta$ small.

• •

When there are many corruptions relative to signal sparsity ($k\sim s$), smaller values of $\lambda$ tend to perform better. Again, this is consistent with the theory in terms of the parameter $\eta$.

4.1.2 Iteratively reweighted $\ell^{1}$ minimization

The results from the previous section show that our a priori postulated optimal values of $\lambda$ are not optimal in practice; this suggests that an adaptive learning of $\lambda$ may produce better results. See, for example, [7]. This section introduces an iteratively reweighted $\ell^{1}$ optimization procedure that effects this learning of $\lambda$.

We compute minimizers $\hat{x}$ and $\hat{c}$ using an initial value of $\lambda$. We then update $\lambda$ based on $\hat{x}$ and $\hat{c}$, and then recompute minimizers with the new $\lambda$. Such an approach not only allows for a single parameter $\lambda$ to be updated, it also permits individual (i.e. non-equal) weights to be used for term in the regularization functional. This aims to enhance recovery performance by both iteratively estimating an optimal weighting $\lambda$ between the coefficients and corruptions term, and iteratively estimating the support sets of $x$ and $c$.

We outline the procedure below:

• •

Step 1. Set $r=1$, $\mu_{i}=1$ for $i=1,\ldots,N$, and $\lambda_{j}=1$ for $j=1,\ldots,m$. Prescribe noise tolerance $\epsilon$ and a small positive number $\eta>0$.

• •

Step 2. Compute the solution $(\widehat{x},\widehat{c})$ to

[TABLE]

where ${\left\|z\right\|}_{1,\mu}=\sum^{N}_{i=1}\mu_{i}|z_{i}|$ and ${\left\|d\right\|}_{1,\lambda}=\sum^{m}_{j=1}\lambda_{j}|d_{j}|$.

• •

Step 3. Update $\mu$ and $\lambda$ as follows:

[TABLE]

• •

Step 4. If $r<r_{\max}$, set $r=r+1$ and go back to step 2, otherwise stop.

Numerical results in the bottom row of plots in Figures 1 and 2 show this approach (implemented with $r_{\max}=10$ iterations) significantly improves the recovery over a fixed choice of $\lambda$. We therefore use this iteratively reweighted $\ell^{1}$ approach for optimization for all our simulations in the next section.

4.1.3 Large corruption values

This section is devoted to understanding the behavior of our algorithm with respect to the magnitude of the corruptions.

We run the same experiment as outlined at the beginning of Section 4.1 on Model 2 (the measurement matrix is a subsampled DFT matrix) using the iteratively reweighted algorithm outlined in Section 4.1.2. For this test, we vary $C$ between $1$ and $10^{6}$, and choose the random variable $Z$ defining the corruptions as a standard Cauchy random variable.666The point of generating from a Cauchy distribution is to show that measurement corruption by heavy-tailed distributions does not adversely affect the algorithm’s results.

A straightforward application of the iteratively reweighted algorithm in Section 4.1.2 when $C$ is very large produces suboptimal results. The reason for this is the scale differential between $x$ and $c$, so that the algorithm heavily favors recovery of the corruptions and devotes little effort to recovering the signal. To overcome this limitation, we leverage a significant advantage of our algorithm: Corruption indices and values in the measurement vector are identified. This allows us to formulate a slight modification of the algorithm in Section 4.1.2:

Run the algorithm from Section 4.1.2, generating computed solutions $\widehat{x}$ and $\widehat{c}$. 2. 2.

If $\|\widehat{c}\|<C_{\mathrm{max}}\|y-\widehat{c}\|$, then return the solutions $\widehat{x}$ and $\widehat{c}$. 3. 3.

If instead $\|\widehat{c}\|\geq C_{\mathrm{max}}\|\widehat{y}-\widehat{c}\|$, then define a support set for the vector $\widehat{c}$ as

[TABLE]

and let $\widehat{c}_{S}$ equal to $\widehat{c}$ on $S$ and zero otherwise. 4. 4.

Remove the large corruptions from the measurements and resolve with the measurements $\widetilde{y}\leftarrow y-\widehat{c}_{S}$. This yields a new solution pair $\widetilde{x}$ and $\widetilde{c}$. Return $x=\widetilde{x}$ and $c=\widehat{c}_{S}+\widetilde{c}$.

This procedure uses the algorithm to identify and remove highly corrupted measurements, and then uses another instance of the algorithm to accurately compute the signal. We use the procedure above with the choices $C_{\mathrm{max}}=10$ and $\tau=\frac{1}{5\sqrt{m}}$.

We can now generate a phase transition plot for a fixed value of $C$. Figure 3 shows the transition plots for values $C=1,10^{3}$, and $10^{6}$. We see that the algorithm detects and removes corruptions just as well when $C=1$ as when $C=10^{6}$.

• Remark 4.1

The iteratively reweighted procedure in (4.13) updates weights for both the corruptions ($\lambda_{i}$) and the signal ($\mu_{i}$). Since our focus here is recovery of the corruptions, one may wonder which set of weights is more influential. We have conducted tests in this direction by performing an experiment parallel to the results in Figure 3, where we iteratively update $\lambda_{i}$ according to (4.13), but keep $\mu_{i}$ fixed at unity for all $i$. Our results, shown in Figure 4, indicate that fixing the weights $\mu_{i}$ results in significant deterioration of the algorithm’s performance when $C=1$. However, it results in notable improvement of the algorithm when $C=10^{3}$ or $C=10^{6}$. In the context of soft faults, the $C=1$ behavior of the algorithm is more relevant since when $C\geq 10^{3}$ it is likely that the corruptions can be easily identified and removed by other means. In this small-$C$ context, allowing both sets of weights $\mu_{i}$ and $\lambda_{i}$ to vary appears to be beneficial. On the other hand, the deterioration of the algorithm for very large $C$ is an interesting phenomenon whose investigation we leave for future work.

4.2 Recovery of compressible polynomial Chaos expansions

In this section we test our algorithm on more realistic problems in UQ: sparse recovery of multivariate polynomial Chaos expansion coefficients with corrupted measurements. Polynomial chaos expansions (PCE) [40, 14] have become a popular means of quantifying parametric uncertainty in expensive computer simulations. To formulate our problem using our earlier notation, let $f(\xi)$ denote a scalar-valued response of a model (e.g., a differential equation) where $\xi\in\mathbb{R}^{d}$ is a random parameter appearing in the model. The dependence of $y$ on $\xi$ thus encodes uncertainty in the response. We are interested in building the approximation $\xi\mapsto\sum_{n=1}^{N}x_{n}\phi_{n}(\xi)$, where $\left\{\phi_{n}\right\}_{n=1}^{N}$ are computable orthonormal polynomials constructed from the probability density of the random vector $\xi$, and we wish to compute the unknown coefficients $x_{n}$. In a CS recovery procedure, we construct $m$ samples $\left\{\xi_{j}\right\}_{j=1}^{m}$ of the random vector $\xi$, collect the measurements $y_{j}=f(\xi_{j})$, and then attempt to find a sparse coefficient vector $x$ minimizing $\|y-Ax\|$, where $A$ is the measurement matrix with entries $(A)_{j,n}=\phi_{n}(\xi_{j})$. The underlying assumption is that $\xi\mapsto f(\xi)$ is expensive to evaluate so that $m$ should be as small as possible. To focus our study on the corruptions problem, we consider the case where the vector $y$ can have a sparse number of entries that are polluted by large-magnitude errors.

The models $f(\xi)$ we consider here reflect the types of large scale models that are susceptible to soft failures. However, these test models can be evaluated repeatedly with almost zero probability of corruptions. Therefore, to simulate the effect of soft failures we randomly generate soft faults according to the corruptions model from Section 4.1. After constructing components of $y$ as $f(\xi)$, we pollute $k$ of these entries as described at the beginning of Section 4.1. In our tests below we fix a value $r\coloneqq k/m$, the ratio of corrupted measurements.

4.2.1 Genz test functions

We compare the algorithm presented in this paper against a classical $\ell^{1}$ minimization approach in the presence of measurement corruptions for the purposes of computing compressible PCE expansion coefficients of a function. A classical $\ell^{1}$ minimization algorithm sets the corruptions vector $d=0$ in (2.3) and minimizes over all $x\in\mathbb{R}^{N}$.

Our function $f(\xi)$ will be one of the multidimensional test functions used by Genz [13]. For $\xi\in\mathbb{R}^{d}$, $d\in\mathbb{N}$, we investigate computing expansion coefficients for the following two functions on the hypercube $[-1,1]^{d}$:

[TABLE]

We use $d=4$ and $d=10$ in our tests, with the dictionary elements $\phi_{n}$ given by tensor-product Chebyshev polynomials of total degree $10$ and $4$, respectively, over $[-1,1]^{d}$. We set the corruptions ratio to the value $r=0.1$ uniformly over all tests, and vary the corruptions magnitude $C$. After computing a coefficient vector $x$ solving either a classical $\ell^{1}$ problem or (2.3), we compute a discrete $\ell^{2}$ error metric defined by

[TABLE]

where $Q=10^{3}$ for each test, and $\tau_{q}$ are iid samples drawn from the product Chebyshev distribution over $[-1,1]^{d}$.

Figure 5 shows the result of this test. (See the figure caption for additional details of the test.) The results indicate that when corruptions are present, a standard $\ell^{1}$ minimization algorithm suffers severe degradation of the quality of the computed expansion coefficients. However, the corruptions algorithm of this paper is able to compute accurate coefficients in the presence of corruptions, whether they have large or small magnitude.

This example shows that there may be a penalty for using our algorithm when no corruptions are present. This is mostly easily noticed in the product peak example with no corruptions ($C=0$): The corruptions algorithm of this paper computes a PCE that is less accurate than the result using a standard $\ell^{1}$ minimization approach. (Compare the black lines in row 3 of Figure 5.)

4.2.2 Damped Harmonic Oscillator

In this section we investigate the fault-tolerance of our algorithm for recovery of PCE coefficients in a damped linear oscillator subject to external forcing with six unknown parameters. The model is

[TABLE]

where we assume the damping coefficient $\gamma$, spring constant $k$, forcing amplitude $g$ and frequency $\omega$, and the initial conditions $u_{0}$ and $u_{1}$ are all uncertain, defining components of a 6-dimensional random vector $\xi$. We solve (4.14) analytically to circumvent the impact of discretization errors in our study.

Defining $\xi=(\gamma,k,g,\omega,u_{0},u_{1})$, we restrict the components $\xi^{(j)}$ of $\xi$ to the following ranges:

[TABLE]

We define $I_{\xi}\in\mathbb{R}^{6}$ to be the range of $\xi$ defined by the product of these intervals. For any parameter realization in $I_{\xi}$ the harmonic oscillator is underdamped. In the following, we choose our quantity of interest as $f(\xi)=u(20,\xi)$. We set the corruptions magnitude $C$ as the mean of the function, i.e. $C=\mathbb{E}_{\xi}[f]$.

Figure 6 compares, as a function of the number of measurements, the error in classical $\ell^{1}$ recovery for uncorrupted sparse recovery versus the iteratively reweighted version of the sparse corruptions $\ell^{1}$ optimization proposed in Section 4.1.2. The results show that the sparse corruptions optimization notably outperforms standard $\ell^{1}$ minimization when corruptions are present, and is competitive without corruptions.

In Figure 7 we run the iteratively reweighted sparse corruptions optimization but vary the corruptions rate $r$, and the corruptions magnitude $C$. The left-hand plot shows predictable behavior: increasing corruptions has deleterious effects on the error in recovery, but notably the algorithm is reasonably stable for increasing $r$. The right-hand plot shows that the algorithm is relatively insensitive to the magnitude of the corruptions.

5 Summary and conclusion

We have developed novel theoretical guarantees and algorithms for recovery of sparse or compressible signals where measurements have been polluted by high-magnitude corruptions. Our results are uniform theoretical recovery estimates for general linear systems where the measurement matrix satisfies a corruptions-based RIP-like condition.

We have refined an existing regularized $\ell^{1}$ minimization algorithm into an iteratively reweighted $\ell^{1}$ minimization algorithm that shows superior performance for the examples that we have investigated. An application of these examples to recovery of polynomial Chaos expansions from model UQ problems illustrates that our algorithms are resistant to highly-corrupted measurement data that may result from hardware or software faults in modern large-scale parallel computing paradigms.

Empirical tests suggest that refinements of our algorithm is relatively stable with respect to the magnitude of the corruptions, but our theory is not applicable to these algorithmic refinements and some observed behavior (e.g., Remark 4.1.3) remains theoretically unexplained, which can be the subject of future explorations.

Acknowledgments

B. Adcock thanks Simone Brugiapaglia and Xiaodong Li for helpful discussions. The authors acknowledge an anonymous referee whose report led to the investigations outlined in Remark 4.1.3.

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA-0003525. The views expressed in the article do not necessarily represent the views of the U.S. Department of Energy or the United States Government.