Compressed Blind De-convolution

TL;DR

This paper introduces novel L1 minimization techniques for blind deconvolution and compressed sensing, enabling recovery of unknown filters and sparse signals from limited noisy measurements under certain conditions.

Contribution

It develops new L1 minimization methods and establishes recovery guarantees for blind deconvolution and compressed sensing with unknown stable systems.

Findings

Exact recovery of filter H and sparse signal u from noisy measurements.

Reconstruction of x, u, and H from O(k log n) measurements under specific conditions.

Proposed LP optimization algorithm reliably estimates unknown filter and sparse input.

Abstract

Suppose the signal x is realized by driving a k-sparse signal u through an arbitrary unknown stable discrete-linear time invariant system H. These types of processes arise naturally in Reflection Seismology. In this paper we are interested in several problems: (a) Blind-Deconvolution: Can we recover both the filter $H$ and the sparse signal $u$ from noisy measurements? (b) Compressive Sensing: Is x compressible in the conventional sense of compressed sensing? Namely, can x, u and H be reconstructed from a sparse set of measurements. We develop novel L1 minimization methods to solve both cases and establish sufficient conditions for exact recovery for the case when the unknown system H is auto-regressive (i.e. all pole) of a known order. In the compressed sensing/sampling setting it turns out that both H and x can be reconstructed from O(k log(n)) measurements under certain technical conditions on the support structure of u. Our main idea is to pass x through a linear time invariant system G and collect O(k log(n)) sequential measurements. The filter G is chosen suitably, namely, its associated Toeplitz matrix satisfies the RIP property. We develop a novel LP optimization algorithm and show that both the unknown filter H and the sparse input u can be reliably estimated.

Figures

Equations

$$x(t)=(Hu)(t)=\sum _{i}u({\tau }_i)h(t-{\tau }_i)$$

$$\min \Vert u{\Vert }_1\text{mbox}subjecttoy=GHu$$

$$y(t)=(g\ast x)(t)=(g\ast h\ast u)(t)=((g\ast h)\ast u)(t)=(h\ast g\ast u)(t)$$

$$z(t)\stackrel{\Delta }{=}y(t)-y(t-1)=g\ast (x(t)-x(t-1))=(g\ast u)(t)$$

$$x(n)+\sum _{i=1}^p{a}_{i}x(t-i)=z(t)$$

$$\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_{m-1}\end{array}\right]=\left[\begin{array}{ccccccc}{g}_{n-m}& g_{n-m-1}& \cdots & g_0& \cdots & 0& 0\\ g_{n-m+1}& g_{n-m}& \cdots & g_1& g_0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & \ddots & 0\\ g_{n-1}& g_{n-2}& \cdots & \cdots & \cdots & g_1& g_0\end{array}\right]\left[\begin{array}{c}{x}_0\\ x_1\\ ⋮\\ x_{n-1}\end{array}\right]$$

$$Y=\left[\begin{array}{ccccc}{y}_p& y_{p-1}& \cdots & y_1& y_0\\ y_{p+1}& y_p& \cdots & y_2& y_1\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ y_{m-1}& y_{m-2}& \cdots & y_{m-p}& y_{m-p-1}\end{array}\right]\left[\begin{array}{c}1\\ a_1\\ ⋮\\ a_p\end{array}\right]=G^{[m-p]}z$$

$$Ya+y^{[m-p]}=G^{[m-p]}z$$

$$y(t)=\sum _{\tau =0}^{n-1}{g}_{t,\tau }x(\tau ),t=0,1,\dots ,m$$

$$\underset{u\in {\mathbb{R}}^n,a\in {\mathbb{R}}^p}{\min }\Vert u{\Vert }_1\text{ subject to}Ya+y^{[m-p]}=G^{[m-p]}u$$

$$\underset{u\in {\mathbb{R}}^n,a\in {\mathbb{R}}^p}{\min }\frac{1}{2}\Vert Ya+y^{[m-p]}-G^{[m-p]}u{\Vert }_2^2+\lambda \Vert u{\Vert }_1$$

$$a=(Y^{T}Y{)}^{-1}{Y}^T\left(G^{[m-p]}u-y^{[m-p]}\right)$$

$$(I-Y(Y^{T}Y{)}^{-1}{Y}^T)y^{[m-p]}=(I-Y(Y^{T}Y{)}^{-1}{Y}^T)G^{[m-p]}u$$

$$\underset{u\in {\mathbb{R}}^n}{\min }\Vert u{\Vert }_1\text{ subject to}P{y}^{[m-p]}=P{G}^{[m-p]}u$$

$$(1-{\delta }_S)\Vert x{\Vert }_2^2\le \Vert G_T^{[m-p]}x{\Vert }_2^2\le (1+{\delta }_S)\Vert x{\Vert }_2^2$$

$$\underset{u\in {\mathbb{R}}^n,a\in {\mathbb{R}}^p}{\min }\Vert u{\Vert }_1\text{ subject to}Ya+y^{[m-p]}=G^{[m-p]}u$$

$$|h(i)|\le M{\rho }^i$$

$$\underset{u\in {\mathbb{R}}^n,a\in {\mathbb{R}}^p}{\min }\Vert u{\Vert }_1\text{ subject to}Ya+y^{[m-p]}=G^{[m-p]}u$$

$$y(n)=x(n)+w(n)$$

$$y(n)+\sum _{i=1}^p{a}_{i}y(n-i)=u(n)+e(n)$$

$$Y=\left[\begin{array}{ccc}0& \cdots & 0\\ y_0& \cdots & 0\\ ⋮& \ddots & ⋮\\ y_{p-1}& \cdots & y_0\\ ⋮& ⋮& ⋮\\ y_{n-2}& \cdots & y_{n-p}\end{array}\right]$$

$$y+Ya=u+e$$

$$\underset{u\in {\mathbb{R}}^n,a\in {\mathbb{R}}^p}{\min }\frac{1}{2}\Vert y+Ya-u{\Vert }_2^2+\lambda \Vert u{\Vert }_1$$

$$X=\left[\begin{array}{ccc}0& \cdots & 0\\ x_0& \cdots & 0\\ ⋮& \ddots & ⋮\\ x_{p-1}& \cdots & x_0\\ ⋮& ⋮& ⋮\\ x_{n-2}& \cdots & x_{n-p}\end{array}\right]$$

$${\hat{u}}_I$$

$${\hat{u}}_{I^c}$$

$$\hat{a}$$

$$x(n)+\sum _{i=1}^p{a}_{i}x(n-i)=u(n)+\sum _{i=1}^q{b}_{i}u(n-i)$$

$$\left[\begin{array}{cccc}{x}_0& 0& \cdots & 0\\ x_1& x_0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ x_p& x_{p-1}& \cdots & x_0\\ ⋮& ⋮& ⋮& ⋮\\ x_{n-1}& x_{n-2}& \cdots & x_{n-p}\end{array}\right]\left[\begin{array}{c}1\\ a_1\\ ⋮\\ a_p\end{array}\right]=\left[\begin{array}{ccccc}1& 0& 0& \cdots & 0\\ b_1& 1& 0& \cdots & 0\\ ⋮& ⋮& \ddots & \ddots & ⋮\\ b_q& \cdots & b_1& 1& \cdots \\ ⋮& ⋮& \ddots & \ddots & ⋮\\ 0& \cdots & b_q& \cdots b_1& 1\end{array}\right]\left[\begin{array}{c}{u}_0\\ u_1\\ ⋮\\ u_{n-2}\\ u_{n-1}\end{array}\right]$$

$$B:=\left[\begin{array}{ccccc}1& 0& 0& \cdots & 0\\ b_1& 1& 0& \cdots & 0\\ ⋮& ⋮& \ddots & \ddots & ⋮\\ b_q& \cdots & b_1& 1& \cdots \\ ⋮& ⋮& \ddots & \ddots & ⋮\\ 0& \cdots & b_q& \cdots b_1& 1\end{array}\right]$$

$$Ya+y^{[m-p]}=G^{[m-p]}Bu$$

Full text

1 Introduction

In this paper we focus on blind de-convolution problems for filtered sparse processes. These types of processes naturally arise in reflection seismology [1]. The LTI system $H$ is commonly referred to as the wavelet, which can be unknown, and serves as the input signal. This input signal passes through the different layers of earth and the reflected signal $z$ corresponds to the reflection coefficients from the different layers. The signal $z$ is typically sparse. The reflected output, which is referred to as the seismic trace, is recorded by a geophone. Other applications of filtered sparse processes include nuclear radiation [2], neuronal spike trains [3] and communications [4].

Specifically, a sparse input $u(t)$ is filtered by an unknown infinite impulse response (IIR) discrete time stable linear filter $H$ and the resulting output

[TABLE]

is measured in Gaussian noise, namely, $y(t)=x(t)+n(t)$ for $t=0,\,1,\,\ldots,N$. The goal is to detect $z(t)$, and estimate the filter $H$. The main approach heretofore proposed for blind de-convolution involves heuristic iterative block decomposition schemes (first proposed in [5]). Here the filter and sparse inputs are alternatively estimated by holding one of them constant. While these algorithms can work in some cases, no systematic performance guarantees currently exist. We explore a convex optimization framework for blind de-convolution.

In addition we consider the compressed sensing problem, namely, $x(t)$ is compressed by means of a random Gaussian filter ensemble, as described in Figure 1 and the resulting output is measured noisily. Analogously, we can consider a random excitation model as in Figure 2. Our task is to detect $z(t)$ and estimate $H$. Our goal is to characterize the minimum number of random samples required for accurate detection and estimation.

1.1 Comparison with Compressed Sensing

Note that this is significantly different from the standard Compressed sensing(CS) [6, 7] problem. In standard CS we have a signal or image, $x\in\mathbb{R}^{n}$, which is sparse in some transform domain. Specifically, there is a known orthonormal matrix $H$ such that the transformed signal $z=H^{T}x$ is k-sparse, namely, has fewer than $k$ non-zero components222This is often referred to as transform domain sparsity.. A matrix $G\in\mathbb{R}^{m\times n}$ then maps $x$ to measurements $y=Gx=GHu$. For suitable choices of matrices $G$, such as those satisfying the so called Restricted Isometry Property (RIP), the k-sparse signal $z$ can be recovered with $O(k\log(n))$ measurements as a solution to a convex optimization problem:

[TABLE]

This result holds for all sparsity levels $k\leq\alpha n,\,\alpha<1$, for sufficiently small $\alpha$. There has been significant effort in CS in recent years leading to various generalizations of this fundamental result. This includes the case when the signal $x$ is approximately sparse (see [8, 9]) and when the measurements are noisy, i.e., $y=GHu+e$ (see  [9]).

This paper is a significant extension of CS to cases where $H$ is not only not orthonormal but also arbitrary and unknown. Specifically, $H$, is a causal discrete linear time invariant system (LTI) with an unknown impulse response function $h(\cdot)$ as described above. A typical signal $x$ is neither sparse nor approximately sparse as we will see in Section 7.

1.2 Our Approach

Our CS problem (schematically shown in Figures 1 2) boils down to determining whether there is a sampling operator $G$ with $O(k\log(n))$ samples such that the signal $x$ can be recovered uniquely from the samples $y=Gx=GHu$ using a convex optimization algorithm. It turns out that this is indeed the case when $H$ is belongs to the set of stable finite dimensional AR processes of a known order.

At first glance the problem as posed appears difficult. For one there is no reason $GH$ satisfies isometry property when $H$ is not orthonormal. To build intuition we describe a practically relevant problem. A specific example is when $x$ is a one-dimensional piecewise constant signal. Such a signal is not sparse but does have a sparse derivative, namely, $u(t)=x(t)-x(t-1)$ is sparse. Clearly, the signal $x$ can represented as an output of an (integral) operator $H$ acting on a sparse input $u$, namely, $x=Hu$. However, $H$ is no longer orthonormal. To account for this scenario one usually minimizes the total variation (TV) of the signal. A compressed sensing perspective for this case has already been developed [10].

We develop an alternative approach here. Suppose we now filter $x$ through an LTI system $G$ whose impulse response is $g(t)$. Mathematically, we have,

[TABLE]

Since, the composite system $g*h$ is LTI we have that,

[TABLE]

Now we are in the familiar situation of $z=Gu$ of the standard CS problem, except that $G$ is a Toeplitz matrix. Consequently, if the Toeplitz matrix $G$ satisfies the RIP property we can recover $z$ using standard tools in CS. Indeed RIP properties of Toeplitz matrices have been studied [11]. Note that this idea generalizes to arbitrary but known finite dimensional stable LTI systems, $H$. The main idea being used here is the commutative property of convolutions.

However, the question arises as to how to deal with unknown system $H$? It turns out that corresponding to every finite dimensional LTI system there is an annihilating filter [12]. If $H$ is a pth order linear dynamical system it turns out that the annihilating filter, $H^{\perp}$ is parameterized by $p$ parameters. Now employing commutativity of convolution, namely, $g*h=h*g$, followed by filtering through the annihilator we are left with a linear characterization of the measurement equations. We are now in a position to pose a related $\ell_{1}$ optimization problem where the parameters are the sparse signal $z$ as well as the parameters governing the annihilating filter. Our proof techniques are based on duality theory.

Strictly speaking, for AR models commutativity is not necessary. Indeed, we could consider general random projections, but this comes at a cost of increasing the number of measurements as we will see later. On the other hand RIP properties for random projections is (provably) significantly stronger than Toeplitz matrices. Nevertheless note that in the random excitation scenario of Figure 2, the structure does not lend itself to a random projection interpretation. For these reasons we consider both constructions in the paper.

The paper is organized as follows. The mathematical formulation of the problem is presented in Section 2. Section 3 describes the new $\ell_{1}$ minimization algorithm. The result for recovery with AR filtered processes (Theorem 1) is stated in this section. The proof of Theorem 1 can be found in Section 5. To help the reader understand the main idea of the proof we first consider a very simple case and Section 5.2 provides the proof for the general case. Section 3.1 addresses the blind-deconvolution problem, which can be regarded as a noisy version of our problem. We use LASSO to solve this problem and the detailed proof is provided in Section 6. In Section 4, we extend the our techniques to two related problems, namely, decoding of ARMA process and decoding of a non-causal AR process. Finally, simulation results are shown in Section 7.

2 Problem Set-up

Our objective is to reconstruct an autoregressive (AR) process $x(t)$ from a number of linear and non-adaptive measurements. An autoregressive model is known as an “all-pole” model, and has the general form

[TABLE]

where $z(t)$ is a sparse driving process. We assume the vector $z=[z_{0},\cdots,z_{n-1}]^{T}$ is $k$-sparse, that is, there are only $k$ non-zero components in $z$. The task of compressed sensing is to find the AR model coefficients $a=[a_{1},\cdots,a_{p}]^{T}$ and the driving process $z=[z_{0},\cdots,z_{n-1}]^{T}$ from the measurement $y$. In this paper, we assume that the AR process $x(t)$ is stable, that is, the magnitude of all the poles of the system is strictly smaller than 1. In later discussion, we use $x_{t}$ or $x(t)$ interchangeably for convenience of exposition.

Note that in standard CS setup, the signal $x$ is assumed to be sparse in some known transform space. However, in our problem, the AR model is assumed to be unknown and the main contribution of this paper is to solve this new problem efficiently.

We consider two types of compressed sensing scenarios:

2.1 Toeplitz Matrices

Here we realize $m$ measurements by applying the sensing matrix $G$ to signal $x=[x_{0},\cdots,x_{n-1}]^{T}$.

[TABLE]

where each entry $g_{i}$ is independent Gaussian random variable $\mathcal{N}(0,1)$ or independent Bernoulli $\pm 1$ entries. Here the Toeplitz matrix $G$ preserves the shift-structure of the signal. Roughly speaking, assume $z^{\prime}$ is a shifted version of $z$ (disregarding the boundary effect), then $Gz^{\prime}$ is also just a shifted version of $Gz$. This is particularly suitable for the random excitation model of Figure 2.

For notational purposes we denote by $x^{[s]}$ (or $G^{[s]}$) to denote the subvector of $x$ (or submatrix of $G$) that is composed of the last $s$ components (or $s$ rows) of $x$ (or $G$). By rearranging the above Equation 2 and using the shift-property of $G$, we have the following equation.

[TABLE]

where we recall that $z=[z_{0},\cdots,z_{n-1}]^{T}$. Now Equation 3 is simplified to

[TABLE]

where $a=[a_{1},\cdots,a_{p}]^{T}\in\mathbb{R}^{p}$ and $z\in\mathbb{R}^{n}$ ($k$-sparse) need to be decoded from the model.

2.2 Random Projections

Here we consider randomly projecting the raw measurements $x(t)$, namely,

[TABLE]

where, each entry $g_{t,\tau}$ is an independent Gaussian random variable $\mathcal{N}(0,1)$ or independent Bernoulli $\pm 1$ entry. The reason for choosing random projections over random filters is that IID random Gaussian/Bernoulli matrix ensembles have superior RIP constants. The optimal RIP constants for toeplitz constructions has not been fully answered. Nevertheless, note that to form the matrix $Y$ with random projections requires significantly more projections. This is because we can no longer exploit the shift-invariant property of convolutions. For instance, consider again the matrix $Y$ of Equation 3 above: if random projections were employed instead of Toeplitz construction the entry $y_{1}$ on row 1 will not be equal to the entry $y_{1}$ in the second row. This means that for a $p$th order model we will require $m\times p$ measurements.

**Notation: ** To avoid any confusion, we use $u^{*}$ to denote the true spike train and $u$ refers to any possible solution in the decoding algorithm. Similarly, $a^{*}$ represents the true coefficients.

3 $\ell_{1}$-minimization Algorithm for AR Models

Since the AR model is unknown, standard decoding algorithms (e.g., Basis Pursuit [8], OMP [13], Lasso [14], etc.) can not be directly applied to this problem. However, we can regard the signal $(u,a)$ (the original signal $u$ together with the unknown coefficients $a$) as the new input to the model and $(u,a)$ is still sparse if $p$ (the length of $a$) is small.

With this in mind we solve the following $\ell_{1}$ minimization algorithm

[TABLE]

More generally, when the measurement $y$ is contaminated by noise, that is, the sensing model becomes $y=Gx+w$ where $w$ is Gaussian noise, the above LP algorithm will be replaced by Lasso,

[TABLE]

where $\lambda$ is a tuning parameter that adapts to the noise level.

Alternatively, the coefficient $a$ can be solved from Equation 4 by taking pseudo-inverse of $Y$,

[TABLE]

Then Equation 4 becomes

[TABLE]

and similar to Equation 5 we can apply the following $\ell_{1}$ minimization to find the solution for $u$.

[TABLE]

where $P$ denotes the projection matrix $I-Y(Y^{T}Y)^{-1}Y^{T}$ and $\|u\|_{1}$ denotes the $\ell_{1}$ norm of $u$. Suppose the solution of Equation 8 is $\hat{u}$. Then $a$ can be easily derived by $\hat{a}=(Y^{T}Y)^{-1}Y^{T}\left(G^{[m-p]}\hat{u}-y^{[m-p]}\right)$ and the signal $x(n)$ can be recovered through Equation 1.

We note that Equation 8 is equivalent to Equation 5 if $Y^{T}Y$ is invertible, which is always assumed to be true in this paper. To summarize the above discussion, our algorithm is summarized below.

(1) Inputs: Measurement $y$, sensing matrix $G$ and order of the system $p$.

(2) Compute $u$ and $a$: Solve the $\ell_{1}$ minimization (Equation 5 or 8) or Lasso (Equation 6).

(3) Reconstruction: Recover the signal $x(n)$ through forward propagation of the AR model of Equation 1.

Before stating the main result, we recall that for every integer $S$ the restricted isometry constant [9, 15], $\delta_{S}$ is defined to be the smallest quantity such that $G^{[m-p]}_{T}$ obeys

[TABLE]

for all subsets $T\subset\{0,1,\cdots,n-1\}$ of cardinality at most $S$ and all $(x_{j})_{j\in T}$.

Note that when the AR filter $a(n)$ is known the result is a direct application of standard compressed sensing results. We state this without proof below for the sake of completion. In other words, if the coefficients $a(n)$ are known, $u^{*}(\cdot)$ is the true driving process in Equation 1 then $u^{*}(\cdot)$ is the unique minimizer of

[TABLE]

A main result of our paper is the following where $a(n)$ is assumed to be unknown. We need the following assumptions before we state our the theorem.

Constant Order: We assume that $p$, the order of AR process $x(n)$, is a constant (i.e., $p$ does not scale with $n,m,S$). 2. 2.

Exponential Decay: Suppose the impulse response $|h(i)|$ of the AR model satisfies

[TABLE]

for some constant $M$ and $0<\rho<1$. 3. 3.

Distance between Spikes: We define the constant $l:=\left(\log(\frac{2}{1-\rho})+p\log(\frac{6\beta_{\max}M}{\beta_{\min}})\right)/\log(\rho^{-1})+p$ and impose the condition that any two spikes, $u_{i}^{*},u_{j}^{*}$ satisfy $|i-j|>l,\,\,i\not=j$. This implies that the sparsity $k:=|\text{Supp}(u^{*})|\leq\min\{S/l,S/3\}$. 4. 4.

Spike Amplitude: We also assume that any spike is bounded, $\beta_{\min}\leq|u_{k}|\leq\beta_{\max},\forall k\in\text{Supp}(u^{*})$.

Theorem 1**.**

Suppose assumptions $1$–$4$ above are satisfied. Let the integer $S$ satisfy $\frac{\delta_{S}}{1-3\delta_{S}}<1$. If $u^{*}(\cdot)$ is the true driving process in Equation 1 then it is the unique minimizer of

[TABLE]

Intuitively speaking, the condition in the theorem requires that the driving process $u(n)$ is sparse enough and any two spikes $(u_{i},u_{j})$ are reasonably far away from each other. This type of assumption is actually also necessary. In section 5.2, we give an example where two spikes are consecutive and show that in this case $x(n)$ can not be solved via equation 8. The proof of Theorem 1 is presented in Section 5.

Remark 3.1*.*

The reader might be curious as to whether a random convolution train provides benefits over random projection. Note that by using random convolutions we can naturally exploit shift-invariance property. Since $Y\in\mathbb{R}^{m-p\times p}$ as in Equation 3 is a partial Toeplitz matrix, we only need $m$ output measurements. In contrast for a random projection, since we can no longer exploit this property, we would require $O(mp)$ measurements.

3.1 Noisy Blind-deconvolution

We consider the noisy blind-deconvolution problem with IID Gaussian noise, $w_{i}\sim\mathcal{N}(0,\sigma^{2})$, and measurements

[TABLE]

where the process $x(n)$ is modeled by $x(n)+\sum_{i=1}^{p}a_{i}x(n-i)=u(n)$. In this section we consider the problem of reconstructing the sparse spike train $u(n)$ and coefficients $a$ from the observed signals $y(n)$. This problem is called “Blind deconvolution” [2, 16] and it is a simplified version of the Compressed Sensing problem where the sensing matrix $G$ is identity matrix. To the best of our knowledge, even this simplified problem is still not completely solved in literature. Therefore, we focus on the uncompressed noisy version here. The noisy compressed version is technically more involved and will be reported elsewhere.

Replacing $x(n)$ with $y(n)-w(n)$ in the AR model, we get

[TABLE]

where we denote $e(n):=w(n)+\sum_{i=1}^{p}a_{i}w(n-i)$.

Again by introducing

[TABLE]

we have the matrix-form system model

[TABLE]

Here Lasso is used to solve the problem:

[TABLE]

We can show that the solution of Lasso is very close to the true $a^{*}$ and $u^{*}$. Before stating the theorem, we first introduce some notation and technical conditions that will be used in the proof.

We denote the noiseless version of $Y$ as

[TABLE]

Denote the support of $u^{*}$ as $I$. We define $X_{1}$ as the matrix comprising of the rows of $X$ indexed by $I$ and $X_{2}$ as the matrix comprising of the rows of $X$ indexed by $I^{c}$. We also denote $x_{\max}=\max_{i}|x_{i}|,\,u_{\min}=\min_{i\in I}|u_{i}|$ and $a_{\max}=\max_{i}|a_{i}|$.

We assume that the AR process $x(n)$ satisfies the following set of conditions.

(1)

The smallest eigenvalue $\lambda_{min}(X_{2}^{T}X_{2})\geq\frac{\|x\|_{2}^{2}}{c}\geq\frac{4np\sigma^{2}}{(\sqrt{2}-1)^{2}}$ for some constant $c>1$.

(2)

$\|X_{1}^{T}\text{sgn}(z_{I}^{*})\|_{\infty}\leq\|x\|_{2}\sqrt{\log n}$,

(3)

$x_{\max}\geq 2\sigma\sqrt{\log n}$ and $\frac{x_{\max}^{2}}{\|x\|_{2}^{2}}\leq\min\{\frac{1}{4c\sqrt{2pn}},\left(\frac{1}{24cp\sqrt{\log n}}\right)^{2}\}$.

In practice, condition (1) is generally satisfied. For instance, if the signal $x$ is persistent, $\frac{1}{\|x\|_{2}}X_{2}^{T}X_{2}$ converges to a constant invertible matrix. Condition (3) is also standard in compressed sensing, which says we need $SNR\geq O(\log n)$. In addition, we also need the assumption that no components are dominantly large (compared with the total energy of $x$). The upper bound for $x_{\max}/\|x\|_{2}$ can be relaxed but the current setup simplifies the analysis.

Condition (2) is new. Let us consider two scenarios. In the first scenario, each spike in $u_{I}$ can be either positive or negative with equal probability (i.e. $\text{sgn}(u_{I})$ is Bernoulli $\pm 1$). In this case, $X_{1}^{T}\text{sgn}(u_{I}^{*})$ behaves like a sub-Gaussian sum and it is usually upper bounded by $\|x\|_{2}\sqrt{\log n}$ with high probability. On the other hand, let us also consider the case when all the spikes in $u_{I}$ are of the same sign, say positive. In this case each entry in $X_{1}^{T}$ and $\text{sgn}(u_{I}^{*})$ is positive and the inner product of these two aligned signals is typically much larger than the first scenario. This phenomena is also illustrated in the experiments shown in Figure 3. In the experiment, the AR model is $x_{t}-1.4x_{t-1}+0.45x_{t-2}=u(t)$. The blue curve corresponds to the scenario when $\text{sgn}(u_{i})$ ($u_{i}$ is a spike) is Bernoulli $\pm 1$. The red curve corresponds to the case when the sign of any spike $u_{i}$ is always $+1$. Each point on the curve is an average over 40 trials. We can see that in the first scenario (blue curve) we can tolerate many more spikes. To the best of our knowledge, this behavior does not exist in standard compressed sensing problem.

Theorem 2**.**

Denote $P:=I-Y(Y^{T}Y)^{-1}Y^{T}$ and assume condition (1),(2) and (3) stated above are satisfied. We also assume parameter $\lambda$ is chosen such that $\lambda\geq 6\sigma pa_{\max}\sqrt{\log n}$ and $u_{\min}\geq 2\lambda$, the solution to Lasso 15 is given by

[TABLE]

and we have $\text{sgn}(u^{*})=\text{sgn}(\hat{u})$ with probability at least $1-8p/n-(p+1)2^{-n/5}$.

**Remark: ** The assumption $u_{\min}\geq 2\lambda$ implicitly implies an SNR bound $O(\log n)$ for the smallest spike. The assumption $\lambda\geq 6\sigma pa_{\max}\sqrt{\log n}$ ensures $\lambda$ to be sufficiently large so that every non-spike element is shrunk to zero by the Lasso estimator. It is hard to analyze the case when parameter $\lambda$ is smaller because in this case it is not clear how to construct $\hat{u}_{I^{c}}$ which is critical for tractable KKT analysis. The choice of $\hat{u}$ in the Theorem 2 is motivated by the proof techniques used in [17]. The proof of Theorem 2 is presented in Section 6.

4 Extensions

In this section, we provide two interesting extensions to the AR model problem. First, we generalize AR model to the autoregressive moving average (ARMA) model, i.e., the process contains both poles and zeros in the transform function. Second, we develop an algorithm for the non-causal AR process, i.e., the current state not only depends on the past inputs but also depends on the future inputs.

4.1 ARMA model

The ARMA model takes the form

[TABLE]

Again we use Equation 2 to obtain the measurement $y=Gx$ where $G$ is a Toeplitz matrix as defined in Section 2. Similar to what we have done in Section 2, we write down the matrix representation of the ARMA model:

[TABLE]

We denote the lower triangular matrix $B$ as

[TABLE]

By multiplying $G^{[m-p]}$ to both sides of Equation 20, we get

[TABLE]

Note that for ARMA model we have an additional term $B$ compared to Equation 4. Generally, matrix $B$ is unknown. We first consider a simple situation when $B$ is assumed to be known to the decoder. Based on Theorem 1 we can derive the following result. .

Theorem 3** (Known Zero Locations).**

Given the same technical conditions as Theorem 1 and assume $u^{*}$ is the original sparse spike train that generates the ARMA process. Then $u^{*}$ is the unique minimizer of

[TABLE]

Proof.

Note that $B$ is also a Toeplitz matrix. From the commutativity of Toeplitz matrix, we have $G^{[m-p]}B=BG^{[m-p]}$. From Section 5, the KKT conditions claim that $u^{*}$ is the unique minimizer of Equation 23 if and only if there exists a vector $\pi$ such that:

$\left(\pi^{T}G^{[m-p]}B\right)_{i}=\text{sgn}(u^{*}_{i})$ for all $i\in\text{Supp}(u^{*})$, 2. 2.

$|\left(\pi^{T}G^{[m-p]}B\right)_{j}|<1$ for all $j\not\in\text{Supp}(u^{*})$, 3. 3.

$\pi^{T}Y=0$.

Applying the commutativity and define $\tilde{\pi}^{T}=\pi^{T}B$, the above three conditions are converted to

$\left(\tilde{\pi}^{T}G^{[m-p]}\right)_{i}=\text{sgn}(u^{*}_{i})$ for all $i\in\text{Supp}(u^{*})$, 2. 2.

$|\left(\tilde{\pi}^{T}G^{[m-p]}\right)_{j}|<1$ for all $j\not\in\text{Supp}(u^{*})$, 3. 3.

$\tilde{\pi}^{T}B^{-1}Y=0$.

Note that both the inverse $B^{-1}$ and the matrix $Y$ are Toeplitz. Therefore, from commutativity, the third equation is equivalent to $\tilde{\pi}^{T}YB^{-1}=0$. Finally, since $B^{-1}$ is invertible, the last equation can be further simplified to $\tilde{\pi}^{T}Y=0$. Now the KKT conditions look exactly the same as those in Section 5. Hence the corollary is proved by following the same argument as in Section 5. ∎

Now we consider the general situation when $B$ is unknown. The difficulty of decoding lies in the fact that we know neither $B$ nor the spike train $u(n)$. There might exist different combinations of $b_{i}$ and $u(n)$ that matches the measurements $y(n)$.

Here we propose an iterative algorithm for estimating $(u,a,b)$ in Equation 22. Each iteration comprises of two basic steps. First, if $B$ is known (from previous iteration), we can use the following $\ell_{1}$ minimization algorithm to solve $u$ and $a$ (Theorem 3).

[TABLE]

Here $\epsilon>0$ is required, even though there may not be any noise, to ensure that we do not get stuck in a local minima.

Now once $u$ is determined we switch from $u$ to $B$, as the optimization variable. This problem reduces to a standard regression problem. First we rewrite Equation 22 as follows:

[TABLE]

which can be simplified to $Ya+y^{[m-p]}=G^{[m-p]}u+G^{[m-p]}Ub$ where we denote

[TABLE]

Now we formulate the following least squares optimization problem:

[TABLE]

In summary our iterative algorithm consists of the following steps:

**Initialization: **

Set $b^{(0)}=0$, i.e., $B^{(0)}=I$.

**Iteration $k$: **

Compute $u^{(k)},a^{(k)},b^{(k)}$

1. 1.

Update $u^{(k)}$ and $a^{(k)}$ via solving Equation 24 with $B=B^{(k-1)}$; 2. 2.

Update $b^{(k)}$ via solving least-square (Equation 25) with $(u,a)=(u^{(k)},a^{(k)})$.

There is a subtlety in the choice of parameter $\epsilon$ in Equation 24. If $\epsilon$ is large, the iterative algorithm appears to have a faster convergence rate but at the cost of significant bias. On the other hand, if $\epsilon$ is small, the convergence rate is slow but the solution has small bias. Therefore, in practical implementation we choose $\epsilon$ to be reasonably large in the early stages of the iteration and then decrease it to $\epsilon/10$ at the later stages of the iteration.

Figure 4 illustrates a concrete example of solving the ARMA model $x(n)-1.9x(n-1)+1.06x(n-2)-0.144x(n-3)=u(n)+0.7u(n-1)+u(n-2)$ by using our iterative algorithm. We choose $\epsilon=3$ in the first 50 rounds of iteration and finally in the last 10 rounds of updates we set $\epsilon=0.3$. Figure 4(b) is a zoom-in version of Figure 4(a) which shows the final stage of the algorithm. We can see the effects of choosing different value of $\epsilon$ as well.

4.2 Non-causal AR model

Many real world signals are non-causal. For example, a 2D image is usually modeled by a Markov random field, where each pixel is dependent on all its neighboring pixels. In this subsection we consider this situation by modeling the signal to be a non-causal AR process.

A non-causal AR model is defined as

[TABLE]

A typical non-causal AR process is shown in Figure 5. Here the impulse response of each spike is two-sided as opposed to the one-sided impulse response of causal AR process. In this subsection, we discriminate between two boundary conditions for the non-causal AR process. As we will show later, there are subtle differences in dealing with these two boundary conditions.

Boundary is circulant, i.e., $x_{0}=x_{n},x_{1}=x_{n+1},\cdots$; 2. 2.

Boundary is not circulant.

4.2.1 Circulant Boundary

In this case we use the following circulant matrix in the sensing model $y=Gx$.

[TABLE]

where $g_{i}$ is i.i.d Gaussian random variable $\mathcal{N}(0,1)$ or Bernoulli $\pm 1$ random variable.

Since the boundary of $x$ is circulant ($x_{-i}=x_{n-i}$), we can write the matrix representation of Equation 26 as

[TABLE]

With an abuse of notation, we use $G^{[i:j]}$ to denote the submatrix of $G$ comprising rows $i$-th through $j$-th of $G$. Now we multiply $G^{[p+1:m-p]}$ to both sides of Equation 28 we get the following equation.

[TABLE]

We define matrix $\tilde{Y}$ to be

[TABLE]

and finally Equation 29 is simplified to

[TABLE]

where $a=[a_{1},\cdots,a_{p},a_{-1},\cdots,a_{-p}]^{T}\in\mathbb{R}^{2p}$.

As in Section 3 we can use either $\ell_{1}$-minimization or Lasso to solve this problem.

[TABLE]

4.2.2 Non-circulant Boundary

The case of non-circulant boundary is slightly more complicated. There are two ways of handling this situation. A simple approach is to view the problem as a perturbation of the circulant boundary case, namely,

[TABLE]

where

[TABLE]

Now one could use Lasso to solve this noisy model:

[TABLE]

Unfortunately, this approach will have a bias. To overcome this limitation, we consider the case where we can make an additional $2p$ set of measurements corresponding to the boundary conditions, namely,

[TABLE]

Then by the denoting

[TABLE]

the sensing model can be simplified to the noiseless version

[TABLE]

Again we can use either $\ell_{1}$-minimization or Lasso to solve this model:

[TABLE]

5 Proof of Theorem 1

We first write down the primal and dual formulation of algorithm 5.

[TABLE]

whose dual formualtion is:

[TABLE]

The proof is based on duality. $u^{*}$ is the unique minimizer of the primal problem 31 if we can find a dual vector $\pi$ with the following properties:

$\left(\pi^{T}G^{[m-p]}\right)_{i}=\text{sgn}(u^{*}_{i})$ for all $i\in\text{Supp}(u^{*})$, 2. 2.

$|\left(\pi^{T}G^{[m-p]}\right)_{j}|<1$ for all $j\not\in\text{Supp}(u^{*})$, 3. 3.

$\pi^{T}Y=0$.

where $\text{sgn}(u^{*}_{i})$ denotes the sign of $u^{*}_{i}$ ($\text{sgn}(u^{*}_{i})=0$ for $u^{*}_{i}=0$) and $\text{Supp}(u^{*})$ denotes the support of vector $u^{*}$. The above set of conditions ensure that the primal-dual pair $(u^{*},\pi)$ is not only feasible but also satisfy the complementary slackness condition, thus optimal. We call the above three conditions as the Dual Optimal Condition (DOC).

The rest of this section is to construct a $\pi$ that satisfies the DOC. Our construction relies on the following result (see [15]).

Lemma 4** ([15]).**

Let $S\geq 1$ be such that $\delta_{2S}\leq\frac{1}{3}$, and $c$ be a real vector supported on $T$ obeying $|T|\leq S$. Then there exists a vector $\pi\in\mathbb{R}^{m}$ such that $\left(\pi^{T}G^{[m-p]}\right)_{i}=c_{i}\,\forall i\in T$. Furthermore, $\pi$ obeys

[TABLE]

This lemma gives us the freedom to choose (arbitrarily) the value of $\pi^{T}G^{[m-p]}$ in the location of $T$ while the magnitude of the rest components is still bounded.

5.1 One Pole Case

In this section we provide a proof for the simple case when $x(n)$ is a first order AR process (i.e., $p=1$) and $u^{*}$ only contains one spike (i.e., every entry of $u^{*}$ is zero except one place). Though simple, it contains the main idea of proof techniques for the more general case. Note that in this simple case the assumptions in Theorem 1 are automatically satisfied.

For the $1$-sparse driving process $u^{*}$, without loss of any generality we assume $u^{*}_{0}=1$ and $u^{*}_{i}=0\,(\forall i\geq 1)$. We also denote $\alpha=-a$ as the root of the characteristic function of the first order AR process. Due to stability we have $|\alpha|<1$. Now in condition 3 of DOC, the term $\pi^{T}Y$ can be recast as

[TABLE]

In Lemma 4, we choose $c$ as $c_{0}=1$, $c_{1}=1/2$ and $c_{j}=0$ ($j=2,\cdots,S-1$). Then Lemma 4 tells us that there exists a $\pi_{1}$ such that $\left(\pi_{1}^{T}G^{[m-p]}\right)_{i}=c_{i}\,(\forall i=0,\cdots,S-1)$ and furthermore

[TABLE]

This implies

[TABLE]

where the summation $\left|\sum_{j=S}^{n-1}\left(\pi_{1}^{T}G^{[m-p]}\right)_{j}\alpha^{j-1}\right|\leq\frac{2|\alpha|^{S-1}}{\sqrt{S}(1-|\alpha|)}\ll\frac{1}{2}$. Therefore $\text{sgn}(\pi_{1}^{T}Y)=1$. To summarize the above discussion, we find $\pi_{1}$ such that:

$\left(\pi_{1}^{T}G^{[m-p]}\right)_{0}=1$ 2. 2.

$|\left(\pi_{1}^{T}G^{[m-p]}\right)_{j}|<1$ for all $j\geq 1$, 3. 3.

$\text{sgn}(\pi_{1}^{T}Y)=1$.

Similarly, by choosing $c_{0}=1$, $c_{1}=-1/2$ and $c_{j}=0$ ($j=2,\cdots,S-1$) in Lemma 4, there exists a $\pi_{2}$ such that condition 1 and 2 of DOC are also satisfied while $\text{sgn}(\pi_{2}^{T}Y)=-1$. Hence, by convexity there exists a $\lambda\in(0,1)$ such that for $\pi=\lambda\pi_{1}+(1-\lambda)\pi_{2}$, it satisfies $\pi^{T}Y=0$ and also condition 1 and 2, i.e., the whole DOC.

Finally we find a primal-dual pair $(u^{*},\pi)$ that satisfy all the feasible constraints and also the complementary slackness condition, which implies $u^{*}$ is the unique minimizer of the primal problem equation 31.

5.2 General Case

In this section we prove that in the general case the three conditions in Theorem 1 ensures the existence of a $\pi$ that satisfies the DOC. Before giving the proof, we point out that if some conditions in Theorem 1 are violated, there might not exist such a $\pi$. Let us consider the case of $p=1$ (first order AR process) and $k$=2 (only two entries of $u(n)$ are nonzero). Moreover, we choose $u^{*}_{0}=u^{*}_{1}=1$ and $u^{*}_{i}=0\,(\forall i>1)$, that is, the two spikes are next to each other.

In this case $[x_{0},x_{1},\cdots,x_{n-1}]^{T}=[1,1+\alpha,\alpha(1+\alpha),\cdots,\alpha^{n-2}(1+\alpha)]^{T}$. We pick $\alpha=-1/2$. Clearly the assumption $|i-i^{\prime}|>l,\forall i,i^{\prime}\in\text{Supp}(u^{*})$ in Theorem 1 is broken. On the other hand, we can also check that there does not exist a $\pi$ that satisfies the whole DOC condition. In fact, suppose $\pi$ is chosen such that condition 1 and 2 are satisfied, then in checking condition 3 we find

[TABLE]

which violates condition 3 in DOC. Hence there does not exist a $\pi$ that satisfies all the three conditions in DOC.

Before proving Theorem 1, we need the following lemma in constructing $\pi$.

Lemma 5**.**

Suppose the assumptions in Theorem 1 are satisfied. Denote $T=\{j+i:j\in\text{Supp}(u^{*}),0\leq i\leq l\}$. Then $|T|\leq S$ and we also have the following inequalities:

(1)

$\forall j\not\in T$* and $i=0,1,\cdots,p$, $|x_{j-i}|<\frac{\beta_{\max}M\rho^{l-p}}{1-\rho^{l}}$,*

(2)

$\forall k\in\cup_{i=1}^{p}\{j-i:j\in\text{Supp}(u^{*})\}$, $|x_{k}|<\frac{\beta_{\max}M\rho^{l-p}}{1-\rho^{l}}$,

(3)

$\forall i=0,1,\cdots,p-1$* and $\forall j\in\text{Supp}(u^{*}),|x_{j}/x_{j+i}|\geq r$ where $r:=\frac{\beta_{\min}(1-\rho^{l})}{\beta_{\max}M}-\rho^{l}$.*

Proof.

First, from the assumption of Theorem 1, $|T|\leq S$. Then we need to verify the three properties.

Suppose $u_{k}$ is a new spike and $k^{\prime}$ be the next spike. Given $i<(k^{\prime}-k)$, we clearly have

[TABLE]

Hence properties (1) and (2) follow.

We denote $\epsilon:=\frac{\beta_{\max}M}{1-\rho^{l}}$. Therefore, for any $j\in\text{Supp}(u^{*})$, $|x_{j}|>|u_{j}|-\epsilon\rho^{l}\geq\beta_{min}-\epsilon\rho^{l}$. Combining with the above argument, we have

[TABLE]

Note that when $\rho^{l}\leq\frac{\beta_{\min}}{3\beta_{\max}M}$ as given in the theorem assumption, we have $r\geq\frac{\beta_{\min}}{3\beta_{\max}M}$. ∎

Remark: Property (1) in Lemma 5 says that many components of $x(n)$ are small. Property (2) ensures that before a new ‘spike’ $u_{j}$ begins ($j\in\text{Supp}(u)$), the amplitude of $x_{j-p},\cdots,x_{j-1}$ is already negligible (i.e., very close to zero) such that the new impulse response caused by $u_{j}$ can be regarded as starting almost from zero level. Finally, property (3) says that when a new spike $u_{j}$ arrives, the corresponding output $x_{j}$ is reasonably large compared to its neighbors.

Now we are ready to prove Theorem 1. Similar to the last section’s argument, the objective is to find a sequence of vectors $\pi_{1},\cdots,\pi_{2^{p}}$ such that any of $\pi_{s}(s=1,\cdots,2^{p})$ satisfies the condition 1 and 2 of DOC while

[TABLE]

and this implies there exists a convex combination $\pi=\sum_{s=1}^{p}\lambda_{s}\pi_{s}$ which satisfies $\pi^{T}Y=0$ and also the condition 1 and 2 of DOC.

Based on Lemma 4, we construct $\pi_{1}$ via fixing the values of $\{(\pi_{1}^{T}G^{[m-p]})_{i}\}_{i\in T}:=\{c_{i}\}_{i\in T}$:

[TABLE]

This choice of $c$ gives the bound $\|c\|_{2}<\sqrt{k+k(1+2^{-1}+2^{-2}+\cdots)}\leq\sqrt{3k}$. Now by applying Lemma 4, we know there exists a $\pi_{1}$ such that $\left(\pi_{1}^{T}G^{[m-p]}\right)_{i}=c_{i}$ when $i\in T$ and

[TABLE]

where the last inequality follows from the assumption of the Theorem. Up to now we have shown that $\pi_{1}$ satisfies condition 1 and 2 of DOC. Next we will check the sign of $\pi_{1}^{T}Y$.

For $t=1,2,\cdots,p$,

[TABLE]

where the magnitude of $A_{t}$ can be lower bounded,

[TABLE]

based on property (3) of Lemma 5. And the magnitude of $B_{t}$ is upper bounded,

[TABLE]

When $l\geq\left(\log(\frac{2}{1-\rho})+p\log(\frac{6\beta_{\max}M}{\beta_{\min}})\right)/\log(\rho^{-1})+p$ as given by the assumption of the theorem, we have $|B_{t}|<|A_{t}|$, which implies that the sign of $(\pi_{1}^{T}Y)_{t}$ is determined by the sign of $A_{t}$.

Hence $\text{sgn}((\pi_{1}^{T}Y)_{t})=\text{sgn}(A_{t})=1$. This implies

[TABLE]

In general, for any sign pattern $[s_{1},\cdots,s_{p}]^{T}(s_{i}\in\{-1,1\})$, by choosing $\{c_{i}\}_{i\in T}$ (compare equation 34) in the following way

[TABLE]

and making similar arguments, we have

[TABLE]

6 Proof of Theorem 2

To prove Theorem 2, we only need to check that $(\hat{u},\hat{a})$ given in the theorem satisfy the KKT conditions. We denote the function $f(u,a)=\frac{1}{2}\|y+Ya-u\|_{2}^{2}+\lambda\|u\|_{1}$. Then the gradient of $f$ with respect to $a$ is

[TABLE]

and the subgradient of $f$ with respect to $u$ is

[TABLE]

where $v$ satisfies $v_{i}=\text{sgn}(u_{i})$ for $i\in I$ and $|v_{i}|<1$ for $i\in I^{c}$. Therefore, we only need to check the following set of (in)equalities

[TABLE]

We first check Equation 35.

Lemma 6**.**

Equation 35 is satisfied with $(\hat{u},\hat{a})$ given in Theorem 2.

Proof.

Actually,

[TABLE]

∎

Next we check Equation 36.

Lemma 7**.**

Equation 36 is satisfied with $(\hat{u},\hat{a})$ given in Theorem 2 with probability at least $1-8p/n-(p+1)2^{-n/5}$.

Proof.

Note that $P$ has the property that $P^{2}=P$ and $PY=0$. Therefore by multiplying $P$ to both sides of Equation 14, we have

[TABLE]

Now we can compute

[TABLE]

where (a) follows from Equation 38 and (b) follows from the fact that $\hat{u}_{I^{c}}=u_{I^{c}}^{*}=0$.

There is a small trick here. Since $y+Y\hat{a}-\hat{u}=P(y-\hat{u})$ as we have shown and $P^{2}=P$, we must have $P(y+Y\hat{a}-\hat{u})=y+Y\hat{a}-\hat{u}$. This implies Equation 36 is correct:

[TABLE]

where (a) follows from $P_{I}^{T}P=P_{I}^{T}$ (i.e., $P^{2}=P$) and the last equality holds true with probability at least $1-8p/n-(p+1)2^{-n/5}$. The proof of last equality is similar to the proof of Lemma 10 and is omitted here. ∎

Verifying inequality 37 requires more effort. We first simplifies the formula for $(y+Y\hat{a}-\hat{u})_{I^{c}}$.

Lemma 8**.**

With $(\hat{u},\hat{a})$ given in Theorem 2, we have

[TABLE]

where we denote $Y_{1}$ as the submatrix comprises of the rows of $Y$ indexed by $I$ and $Y_{2}$ as the submatrix comprises of the rows of $Y$ indexed by $I^{c}$.

Proof.

Following from the proof of Lemma 7, we have ,

[TABLE]

To simplify the above equation, we introduce $P_{11}\in\mathbb{R}^{k\times k}$ as the matrix comprises of the rows of $P$ indexed by $I$ and the columns of $P$ indexed by $I$. Similarly, $P_{12}\in\mathbb{R}^{(n-k)\times k}$ is the matrix comprises of the rows of $P$ indexed by $I$ and the columns of $P$ indexed by $I^{c}$; $P_{22}\in\mathbb{R}^{(n-k)\times(n-k)}$ is the matrix comprises of the rows of $P$ indexed by $I^{c}$ and the columns of $P$ indexed by $I^{c}$. By this definition, after some column and row permutations, $P$ can be rewritten as

[TABLE]

It is easy to check that $P_{11}=P_{I}^{T}P_{I}$ and $P_{12}^{T}=P^{T}_{I^{c}}P_{I}$ (since $P^{2}=P$). Furthermore,

[TABLE]

Hence, Equation 40 can be simplified to

[TABLE]

We note that $P_{11},P_{12},P_{22}$ can be expressed in terms of $Y,Y_{1}$ and $Y_{2}$.

[TABLE]

Moreover $P_{11}^{-1}$ can be derived via matrix inversion lemma:

[TABLE]

Finally, we get

[TABLE]

where (a) follows from the fact that $Y^{T}Y=Y_{2}^{T}Y_{2}+Y_{1}^{T}Y_{1}$. And similarly by repeatedly using this fact we can find the following simplification

[TABLE]

∎

In order to justify the condition 37, we also need the following lemma.

Lemma 9**.**

The following three claims hold true:

(i)

w.p. at least $1-p\cdot(4/n+2^{-n/5})$, $\|Y_{2}^{T}e_{I^{c}}\|_{\infty}\leq 2\sqrt{n\log n}\sigma a_{\max}x_{\max}\sqrt{2p}$.

(ii)

w.p. at least $1-4p/n$, $\|Y_{1}^{T}\lambda\text{sgn}(z_{I}^{*})\|_{\infty}\leq 2\lambda\|x\|_{2}\sqrt{\log n}$

(iii)

w.p. at least $1-2^{-n/5}$, $\lambda_{\max}\left((Y_{2}^{T}Y_{2})^{-1}\right)\leq 2\lambda_{\max}\left((X_{2}^{T}X_{2})^{-1}\right)\leq\frac{2c}{\|x\|_{2}^{2}}$

Proof.

To prove (i), we try to bound the first component $\left(Y_{2}^{T}e_{I^{c}}\right)_{1}$. By definition, the first column of $Y$ equals $[0,\,y_{0},\cdots,y_{n-2}]^{T}=[0,\,x_{0},\cdots,x_{n-2}]^{T}+[0,\,w_{0},\cdots,w_{n-2}]^{T}$. We also remember $e_{i}=w_{i}+\sum_{j=1}^{p}a_{j}w_{i-j}$ where $w_{i}$ are i.i.d. Gaussian $\mathcal{N}(0,\sigma^{2})$. Hence, we have

[TABLE]

It is easy to check that the first term of RHS is zero-mean Gaussian random variable with variance $\leq pa_{\max}^{2}x_{\max}^{2}n\sigma^{2}$. It is well known that for standard Gaussian random variable $t$, $\Pr(|t|\geq a)\leq 2e^{-a^{2}/2}$. So we conclude that with probability $\geq 1-2/n$

[TABLE]

It also can be proved that with probability $\geq 1-2/n-2^{-n/5}$

[TABLE]

We notice that $\sigma a_{\max}x_{\max}\sqrt{2pn\log n}\geq 2pa_{\max}\sigma^{2}\sqrt{n\log n}$ and hence claim (i) follows.

Next, we prove claim (ii). Again, $\left(Y_{1}^{T}\lambda\text{sgn}(z_{I}^{*})\right)_{1}$ can be decomposed into two terms;

[TABLE]

The first term is bounded from the assumption and the second term is Gaussian which is bounded by $\lambda\sigma\sqrt{2n\log n}\leq\lambda\|x\|_{2}\sqrt{\log n}$ (assumption (3) in Subsection 3.1) w.p. $\geq 1-2/n$.

For (iii), we only need to show that with high probability $\lambda_{\min}\left(Y_{2}^{T}Y_{2}\right)\geq\frac{1}{2}\lambda_{\min}\left(X_{2}^{T}X_{2}\right)$, or $\sigma_{\min}\left(Y_{2}\right)\geq\frac{1}{\sqrt{2}}\sigma_{\min}\left(X_{2}\right)$ where $\sigma_{\min}(A)$ denotes the smallest singular value of $A$.

We denote the Gaussian noise matrix

[TABLE]

and call $W_{2}$ as the submatrix that comprises of the rows of $W$ indexed by $I^{c}$. Then, we have

[TABLE]

So the remaining work is to upper bound $\sigma_{\max}(W_{2})$. A tight bound in this case is very difficult. However, the following bound is good enough for our proof. By denoting $W_{2,i}$ as the $i$-th column of $W_{2}$, we have

[TABLE]

where the second last inequality follows from Cauchy-Schwartz inequality. Then by the tail probability of $\chi^{2}$ distribution, we have with probability $1-2^{-n/5}$,

[TABLE]

Then by applying assumption (1) in Subsection 3.1 we have proved the claim (iii). ∎

Finally, we can show that $(\hat{u},\hat{a})$ satisfies the condition 37.

Lemma 10**.**

Equation 35 is satisfied with $(\hat{u},\hat{a})$ given in Theorem 2 with probability at least $1-8p/n-(p+1)2^{-n/5}$.

Proof.

From the tail probability of standard Gaussian $\Pr(|t|\geq a)\leq 2e^{-a^{2}/2}$, we know that with probability at least $1-2/n$, $\max_{i}|w_{i}|\leq 2\sigma\sqrt{\log n}$. Therefore the $\ell_{2}$ norm of all the rows of $Y_{2}$ is upper bounded $\sqrt{p}(x_{\max}+2\sigma\sqrt{\log n})$ with probability at least $1-2/n$. Combined with claim (iii) in Lemma 9, we know that the $\ell_{2}$ norm of all the rows of $Y_{2}(Y_{2}^{T}Y_{2})^{-1}$ is upper bounded $\frac{2c\sqrt{p}}{\|x\|_{2}^{2}}(x_{\max}+2\sigma\sqrt{\log n})\leq\frac{4c\sqrt{p}\,x_{\max}}{\|x\|_{2}^{2}}$ with probability at least $1-2/n-2^{-n/5}$.

Now we can verify that both $-Y_{2}(Y_{2}^{T}Y_{2})^{-1}Y^{T}_{1}\lambda\text{sgn}(u^{*}_{I})$ and $(I-Y_{2}(Y_{2}^{T}Y_{2})^{-1}Y^{T}_{2})e_{I^{c}}$ are small.

First, based on claim (ii) in Lemma 9, with probability at least $1-\frac{2+4p}{n}-2^{-n/5}$

[TABLE]

where the last inequality follows from condition (3) in Subsection 3.1.

Next, it is easy to bound $\|e_{I^{c}}\|_{\infty}\leq 2\sigma pa_{\max}\sqrt{\log n}\leq\lambda/3$ with probability at least $1-2/n$. Also, we have with probability at least $1-\frac{4p+2}{n}-(p+1)2^{-n/5}$

[TABLE]

where the last inequality follows from claim (i) of Lemma 9, condition (3) in Subsection 3.1 and the assumption $\lambda\geq 6\sigma pa_{\max}\sqrt{\log n}$. ∎

7 Numerical Experiments

We present simulations for some interesting cases. Theorem 1 asserts that as long as RIP is satisfied, stability assumptions on $H$ hold, and the spikes are well separated, our $\ell_{1}$-minimization algorithm reconstructs the AR process correctly. For general IID Gaussian or Bernoulli matrix ensemble (not Toeplitz), it is well known that [9] $m\geq O(S\log(n/S))$ ensures good RIP property. However, for our specific Toeplitz structured sensing matrix (Equation 2), this question (when RIP is satisfied) has not been fully answered.

We nevertheless experiment with Toeplitz constructions. First we simulate our algorithm for a third order process. The results are depicted in Figure 6. We see that the reconstruction reproduces both the spike train as well as the filtered process accurately. For the purpose of depiction we added a small amount of noise.

First, we fix the size of sensing matrix ($m=50,n=200$) and choose the entries of sensing matrix $G$ to be Gaussian. We also fix the order of the AR model ($p=2$) and let the sparsity $k$ vary from $1$ to $20$. For each fixed $k$, we run our $\ell_{1}$-minimization algorithm $50$ times to obtain the average performance. The result is shown in Figure 7(a). Similarly, we can choose the sensing matrix $G$ to be Bernoulli $\pm 1$ and do the same experiment again. The result is shown in Figure 7(b). We can see that in this example Toeplitz Bernoulli matrix is more preferable than Toeplitz Gaussian matrix.

Next, we run our algorithm on a case that does not satisfy our assumptions on stability. Specifically we consider the situation when the true process is governed by the equation $x(n)-x(n-1)=u(n)$. This type of model is closely associated with problems that arise when one is interested in minimizing total variations. Note that in this model $\alpha=1$ and it does not satisfy the assumptions of Theorem 1 where we assume $\alpha_{\max}<1$. We adopt the same sensing matrix as the last experiment (Gaussian or Bernoulli) and the empirical success rate of this experiment is shown in Figure 8.

Finally we test how the order of the AR process influences the performance of the algorithm. In this experiment, we fix the size of the sensing matrix as $80\times 200$ and also fix the sparsity $k=10$ (i.e., the $\#$ fraction of nonzero components in $z$ is $5\%$). We let $p$ (order of the AR process) vary from $1$ to $15$. Figure 9(a) shows that empirical success rate for the Gaussian sensing matrix and Figure 9(b) shows that success rate for the Bernoulli sensing matrix. We can see that again Bernoulli Toeplitz matrix outperforms the Gaussian Toeplitz matrix.