Sparse Reconstruction with Compressed Sensing

Author: Ngoc Ly (SysEn 5800 Fall 2021)

Compressed Sensing (CS)

^{Compressed Sensing summary here}

Compression is synonymous with sparsity. So when we talk about compression we are actually referring to the sparsity. We introduce Compressed Sensing and then focus on reconstruction.

Introduction

sub module goal

The goal of compressed sensing is to interact with an underdetermined linear system in which the number of variables is much greater than the number of observations, resulting in an infinite number of signal coefficient vectors $x$ for the same set of compressive measurements $y$. As a result, additional information is necessary for $x$ to recover from $y$. The objective is to reconstruct a vector $x$ in a given of measurements $y$ and a sensing matrix A. Instead of taking a large number of high-resolution measurements and discarding the majority of them, consider taking way fewer random measurements and reconstructing the original $x$ with high probability from its sparse representation.

sub modual

Begin with a linear equation $y=Ax+e$, where $A\in \mathbb {R} ^{M\times N}$ is a sensing matrix that must be obtained and will result in either exact or approximated optimum solution depending on how it is chosen, $x\in \mathbb {R} ^{N}$ is a signal vector with at most $k$-sparse entries, which means $x$ has $k$ non-zero entries, $[N]=\{1,\dots ,N\}$be an index set, $y\in \mathbb {R} ^{M}$ is a compressed measurement vector, $[M]=\{1,\dots ,M\}$, $e\in \mathbb {R} ^{M}$ is a noise vector and assumed to be bounded $\|e\|_{2}\leq \eta$ if it exists, and $M\ll N$.

The goal of compressed sensing is to being with the under determined linear system $y=\Phi x+e$, Where $\Phi \in \mathbb {R} ^{M\times N}$ for $M\ll N$ How can we reconstruct x from The goal is to reconstruct $x\in \mathbb {R} ^{N}$given $y$ and $\Phi$ Considerably fewer random measurements and reconstruct the original $x$ with high probability from its sparse representation instead of taking a large number of high-resolution measurements and discarding the majority of them. being a random matrix

let $[N]=\{1,\dots ,N\}$be an index set $[N]$ enumerates the columns of $\Phi$ and $x$. $\Phi$ is an under determined systems with infinite solutions since $M\ll N$. Why $\ell _{2}$ norm won't give sparse solutions, where asl $\ell _{1}$ norm will return a sparse solution.

Notation =

$x\in \mathbb {R} ^{N}$often not really sparse but approximately sparse

$\mathbf {x} \in \mathbb {R} ^{N}$

$\Phi \in \mathbb {R} ^{M\times N}$for $M\ll N$ Sensing matrix a Random Gaussian or Bernoulli matrix

$y\in \mathbb {R} ^{M}$are the observed y samples

$e\in \mathbb {R} ^{M}$noise vector $\|e\|_{2}\leq \eta$

put defn of p norm here

$x=\Psi \alpha$ where $\Psi$ is the sparsifying matrix and $\alpha$ are coeficients

sub module sparsity

A vector $x$ is said to be $k$-sparse in $\mathbb {R} ^{N}$ if it has at most $k$ nonzero coefficients. The support of $x$ is $supp(x)=\{i\in [N]:x_{i}\neq 0\}$, and $x$ is a $k$-sparse signal when $|supp(x)|\leq k$. The set of $k$-sparse vectors is denoted by $\Sigma _{k}=\{x\in \mathbb {R} ^{N}:\|x\|_{0}\leq k\}$. Consequently, there are ${\binom {N}{k}}$ different subsets of $k$-sparse vectors. If a random $k$-sparse $x$ is drawn uniformly from $\Sigma _{k}$, its entropy $\log {\binom {N}{k}}$ is approximately equivalent to $k\log {\frac {N}{k}}$ bits are required for compression of $\Sigma _{k}$ ~cite(Measurements vs Bits).

The idea is to search for the sparsest $x\in \Sigma _{k}$ from the measurement vector $y\in \mathbb {R} ^{M}$ and a sensing matrix $A\in \mathbb {R} ^{M\times N}$ with $M\ll N$. If the number of linear measurements is at least twice as its sparsity $x$, i.e., $M\geq 2k$, then there exists at most one signal $x\in \Sigma _{k}$ that satisfies the constraint $y=Ax$ and produce the correct result for any $x\in \Sigma _{k}$ [coluccia2015 book7]. Hence, the reconstruction problem can be formulated as an $l_{0}$"norm" program.

${\hat {x}}={\underset {x\in \Sigma _{k}}{argmin}}\|x\|_{0}\quad s.t.\quad y=Ax$

This optimization problem minimizes the number of nonzero entries of $x$ subject to the constraint $y=Ax$, that is to find the sparsest element in the affine space $\{x\in \mathbb {R} ^{N}:Ax=y\}$ [2019 book33]. It turns out to be a combinatorial optimization problem, which is NP-Hard because it includes all possible sets of $k$-sparse out of $[N]$. Furthermore, if noise is present, the recovery is unstable [Buraniuk "compressed sensing"].

In other words, only the smallest support set of $x$ is of interest, i.e., $min(supp(x))$.

The goal is to search for the sparsest $x\in \Sigma _{k}$ given the meassurment $y$ and the constraint matrix $\Phi$.

This is antiquated to find the $\min \|x\|_{0}$ in the set $\{x\in \mathbb {R} ^{N}:\Phi x=y\}$

This searching problem can be formulated to the following $\ell _{0}$ program

Another words we are interested in the smallest $supp(x)$ , i.e. $min(supp(x))$

sub module

Let $\Phi \in \mathbb {R} ^{M\times N}$ satisfy RIP, Let $[N]$ be an index set For $s$ is a restriction on $\mathbf {x}$ denoted by $x_{|s}$ $x\in \mathbb {R} ^{N}$ to $s$ k-sparse $\mathbf {x}$ s.t. RIP is satisfied the $s=supp(\mathbf {x} )$ i.e. $s\subseteq [N]$and $\Phi _{|s}\subseteq \Phi$ where the columns of $\Phi _{|s}$ is indexed by $i\in S$

In search for a unique solution we have the following $\ell _{0}=|supp(x)|$ optimization problem.

sub module zero norm program

$\mathbf {\hat {x}} ={\underset {\Sigma _{k}}{argmin}}\|\mathbf {x} \|_{0}\quad s.t.\quad \mathbf {y} =\Phi \mathbf {x}$, which is an combinatorial NP-Hard problem. Hence, if noise is presence the recovery is not stable. [Buraniuk "compressed sensing"]

Restricted Isometry Property (RIP)

A matrix $A$ is said to satisfy the RIP of order $k$ if for all $x\in \Sigma _{k}$ has a $\delta _{k}\in [0,1)$. A restricted isometry constant (RIC) of $A$ is the smallest $\delta _{k}$ satisfying this condition [2019 book38, coluccia2015 book7].

$(1-\delta _{k})\|x\|_{2}^{2}\leq \|Ax\|_{2}^{2}\leq (1+\delta _{k})\|x\|_{2}^{2}$

Under projections through matrix $A$, the restricted isometry property allows $k$-sparse vectors to have unique measurement vectors $y$. If $A$ meets RIP, then $A$ does not send two distinct $k$-sparse $x\in \Sigma _{k}$ to the same measurement vector $y$, indicating that $x$ is a unique solution under RIP.

If the matrix $A\in {\mathcal {R}}^{M\times N}$ satisfies the RIP condition of order $2k$ and the constant $\delta _{2k}\in [0,1)$, there are two distinct $k$-sparse vectors in $\Sigma _{2k}$. When they are equal, the restricted isometry property holds. If $A$ is a $2k$-order RIP matrix, it means that no two $k$-sparse vectors are mapped to the same measurement vector $y$ by $A$. In other words, when working with sparse vectors, the RIP ensures that the columns of $A$ are nearly orthonormal. Furthermore, $A$ is an approximately norm-preserving function, which means that it preserves its distance when mapping for $k$-sparse signals for all or more as $\delta _{k}$ approaches zero. [Candes, Romberg, Tao[4]] demonstrate that if $x$ is $k$-sparse, and $A$ satisfies the RIP of order $2k$ with RIP-constant $\delta _{2k}<{\sqrt {(}}2)-1$, then $l_{1}$ gives a unique sparse solution. The $l_{1}$ convex optimization problem is the same as the solution to the $l_{0}$ program and can solve by the Linear Program [2019 book38, coluccia2015 book7]. Hence, the $l_{1}$ reconstruction problem is as followed.

${\hat {x}}={\underset {x\in \Sigma _{k}}{argmin}}\|x\|_{1}\quad s.t.\quad y=Ax$

If $A$ satisfies RIP then $A$ doesn't send two distinct k-sparse $x\in \Sigma _{k}$ to the same measurment vector $y$. Anothers words $x$ is a unique solution under RIP.

RIP defined as

$\Phi$ satisfies RIP of order $k$ if for $\forall x\in \Sigma _{k}$ $\delta _{k}\in [0,1)$ satisfies for inequalaty

1. TODO switch s to k

sub module RIP matracies

If a sensing matrix $\Phi$must satisfy RIP, then the number of measurements $M={\mathcal {O}}(K/log(N/K))$ is required recover $x$ with high probability.

sub module RIC

Restricted Isometry Constant (RIC) is the smallest $\delta _{|s}$ in $\{\delta _{k}\in [0,1):(1-\delta _{s})\|x\|_{2}^{2}\leq \|\Phi x\|_{2}^{2}\leq (1+\delta _{s})\|x\|_{2}^{2}\}$

if $M\geq 2k$ i.e. twise the sparsity, then there exists an unique $x$ such that $y=\Phi x$.

for $x$ is k-sparse and and $\Phi$ satisfies RIP of order $2k$ RIC $\delta _{2k}<{\sqrt {2}}-1$ then the $l_{0}$ program can have to relaxed convex form $\ell _{1}$ program.

If $\Phi$ satisfies RIP and $\mathbf {y}$ is sparse the $\ell _{0}$ gives sparse solutions and is a unique. It is equivalent to the following $\ell _{1}$ convex optimization problem and can solve by Linear Program.

sub problem 1 norm program

From Results of Candes, Romberg, Tao, and Donoho

$\mathbf {\hat {x}} ={\underset {s}{argmin}}\|\mathbf {x} \|_{1}\quad s.t.\quad \mathbf {y} =\Phi \mathbf {x}$

Theory

Two things need to be considered when recovering $x$

• (1) The design of the sensing matrix

• (2) The recovery algorithm

Sensing Matrix

Although the sensing matrix $A$ satisfies RIP of order $k$ in some situations, confirming a given matrix $A$ meets RIP's criteria is NP-hard in general. As a result, designing an efficient sensing matrix is critical. These sensing matrices are responsible for signal compression at the encoder end and accurate or approximate reconstruction at the decoder end. For signal compression, different sensing matrices are utilized in compressed sensing. There are random, deterministic, structural, and optimized sensing matrices are used in compressed sensing [2020_book].

• Random Sensing Matrices

Some classes of random matrices satisfy RIP, specifically those matrices with independent and identically distributed (i.i.d.) entries drawn from a sub-Gaussian distribution. It requires the number of measurements, $M=O(klog(N/k))$, to recover $x$ with high probability. Other popular random sensing matrices are Gaussian, Bernoulli, or Rademacher distributions. However, those random dense matrices incur tremendous computational and memory costs, making them unsuitable for large-scale applications. Several researchers have turned to sparse measurement matrices such as binary or bi-adjacency matrices rather than random matrices to address this issue. Nonetheless, those sparse matrices are unstructured in the same way that acyclic networks or trees are. Fortunately, the new random Weibull matrices [9] are built with appropriate observations and provide exact sparse signal reconstruction with a greater probability [2020_book].

• Deterministic Sensing Matrices

Although deterministic sensing matrices are insufficient to satisfy the RIP condition, they can be used to provide instances for novel concepts. The Vandermond $k\times N$ matrices are one of the best deterministic matrices for recovering the $k$-sparse signal; however, the reconstruction technique gets unstable as $N$ rises. Despite the lack of RIP support in the deterministic sensing matrices, it successfully recovered the original sparse signal for the chirp function-based employing $M\times M^{2}$ complex-valued deterministic matrices. According to the researchers Amini and Marvasti [13], binary, bipolar, and ternary matrices are deterministic constructions of sensing matrices that satisfy the RIP of order $k$. Because of the cyclic feature, the reconstruction process can be sped up using a fast Fourier transform (FFT) technique. Another type of deterministic CS sensing matrix is one based on mutual coherence. The finite geometry-based sparse binary matrices were constructed with low coherence [18,19]. The sparseness property of matrices aids in reducing storage requirements and improving the reconstruction process [2020_book].

• Structural Sensing Matrices

Sensing technologies require structured measurement matrices to perform a variety of tasks. Those matrices can be easily created with a small number of parameters. Furthermore, structured matrices can be used to speed the recovery performance of algorithms, making these matrices suitable for big data challenges. Many researchers created Toeplitz and circulant random sensing matrices used in multipath sparse channel estimation and network systems. In terms of estimated accuracy, signal reconstruction speed, and coherence, these matrices perform similar to i.i.d. Gaussian satisfies RIP with high probability. The new sparse block circulant matrix (SBCM) structure greatly reduces computational complexity. Other structural sensing matrices include the Hadamard matrix, which provides near-optimal assurances of recovery while requiring less complexity and so allowing for simple hardware implementation [2020_book].

• Optimized Sensing Matrices

In order to accomplish high-quality signal reconstruction, the sensing matrices must satisfy RIP. Regardless of what might be expected, the RIP is difficult to verify. Another method for validating the RIP is to compute the mutual coherence between the sensing matrix and the sparse matrix or figure the Gram matrix as $G=D^{T}D$. The objective is to improve the sensing matrix to lower coherence using numerous strategies such as the random-detecting framework-based improvement strategy using the shrinkage technique and the irregular estimation of sensing matrix improvement employing a symmetrical strategy. Another technique is to optimize the sensing matrix using a block-sparsified dictionary approach, decreasing the total inter-block and sub-block coherence of the dictionary matrix and thereby significantly increasing the reconstruction [2020_book].

Check if $\Phi$ satisfies RIP Checking $\Phi$ satisfies RIP is combinatorial hard in general so it's unreasonable to ask a computer to verify a matrix satisfies RIP. In order to get around this problem, we need an understanding of what matrices satisfy RIP and recover $x$ with high probability.

• Random Sensing matrices: Gaussian, Bernoulli, Rademacher
• Deterministic Sensing Matrices: binary, bipolar, ternary, Vandermond
• Structural Sensing Matrices: Toeplitz, Circulant, Hadamard
• Optimized Sensing Matrices: (Parkale, Nalbalwar, Sensing Matrices in Compressed Sensing)

Are some examples. Different sensing matrices are more suited for different problems, but in general, we want to use an alternative to Gaussian because it reduces the computational complexity.

Verification of the Sensing matrix

Definition Mutual Coherence

Let $\Phi \in R^{M\times N}$, the mutual coherence $\mu _{\Phi }$ is defined by:</math>

$\mu _{\Phi }={\underset {i\neq j}{\frac {|\langle a_{i},a_{j}\rangle |}{\|a_{i}\|\|a_{j}\|}}}$

1. TODO switch s to k

$(1-\mu )\|x\|_{2}^{2}\leq \|\Phi x\|_{2}^{2}\leq (1+\mu )\|x\|_{2}^{2}$

Welch bound $\mu _{\Phi }\geq {\sqrt {\frac {n}{m(n-m)}}}$> $\mu \geq {\sqrt {\frac {N-M}{M(N-1)}}}$> is the coherence between $\Phi$ and $\Psi$ We want a small $\mu _{\Phi }$ because it will be close to the normal matrix, which satisfies RIP. Also, $\mu _{\Phi }$ will be needed for the step size for the following IHT.

Need to make the connection of Coherence to RIP and RIC.

1. TODO switch s to k

$(1-\mu )\|x\|_{2}^{2}\leq \|\Phi x\|_{2}^{2}\leq (1+\mu )\|x\|_{2}^{2}$

Algorithms

Three big groups of algorithms are:

• Optimization methods: includes $\ell _{1}$ minimization i.e. Basis Pursuit, and quadratically constraint $\ell _{1}$

minimization i.e. basis pursuit denoising.

• Greedy methods: include Orthogonal matching pursuit and Compressive Sampling Matching Pursuit (CoSaMP)

• thresholding-based methods: such as Iterative Hard Thresholding(IHT) and Iterative Soft Thresholding, Approximate IHT or AM-IHT, and many more.

More cutting-edge methods include dynamic programming.

We will cover one, i.e. IHT. WHY IHT THEN? Basis pursuit, matching pursuit type algorithms belong to a more general class of iterative thresholding algorithms. So IHT seems like the ideal place to start. If everything compliment with RIP, then IHT has fast convergence.

Algorithm IHT

The $\ell _{1}$ convex program mentioned in introduction has an equivalent nonconstraint optimization program.

${\underset {y}{min}}\|\mathbf {y} -\Phi \mathbf {x} \|_{2}^{2}+\lambda \|\mathbf {y} \|_{0}$ (cite IT for sparse approximations)  ??? ${\hat {\mathbf {x} }}=arg{\underset {s}{min}}{\frac {1}{n}}\|\mathbf {y} -\Phi \mathbf {x} \|_{2}^{2}+\lambda \|\mathbf {x} \|_{1}$ . In statistics we call the $l_{1}$ regularization LASSO with $\lambda$ as the regularization parameter. This is the closest convex relaxation to $l_{0}$ the first program menttioned in the introduction.[The Benefit of Group Sparsity]

$z_{v}^{(n)}=\nabla f_{v}(x^{(n)})=-\Phi _{v}^{T}(\mathbf {y} -\Phi \mathbf {x} )$ Then $x^{n+1}={\mathcal {H}}\left(\mathbf {x} ^{(n)}-\tau \sum _{j\in N}^{N}z_{v}^{(n)}\right)$

sub modual

Define the threashholding operators as: ${\mathcal {H}}_{s}[\mathbf {x} ]={\underset {z\in \sum _{s}}{argmin}}\|x-\Phi \mathbf {x} \|_{2}$ selects the best-k term approximation for some k

Stopping criterion is $\|y-\Phi \mathbf {x} ^{(n)}\|_{2}\leq \epsilon$ iff RIC $\delta _{3s}<{\frac {1}{\sqrt {32}}}$

• Input $\Phi ,\mathbf {y} ,\mathbf {e} \ {\mbox{with}}\ \mathbf {y} =\mathbf {\Phi } \mathbf {x} |\mathbf {e} and{\mathfrak {M}}$
• output $IHT(\mathbf {y} ,\mathbf {\Phi } ,{\mathcal {S}})$
• Set $x^{(0)}=\mathbf {0}$
• While Stopping criterion false do

• • $x^{(n+1)}\leftarrow {\mathcal {H}}_{|s}\left[x^{(n)}+\Phi ^{*}(\mathbf {y} -\mathbf {\Phi x} ^{(n)})\right]$

• • $n\leftarrow n+1$
  • end while
• return: $IHT(\mathbf {y} ,\mathbf {\Phi } ,{\mathfrak {M}})\leftarrow \mathbf {x} ^{(n)}$

$\Phi ^{*}$ is a Adjoint matrix i.e. the transpost of it's cofactor.

Numerical Example

$x\in \mathbb {R} ^{2}0$

$x=\left[0.1779,1.1649,-0.8262,-1.014,-0.7975,0.1883,0.2648,-0.9036,0.9409,0.4295\right]$

$A=\left[{\begin{matrix}1.24302&0.512088&0.727509&-0.621372&-0.167216&1.28429&1.80216&-0.860403&-0.190791&0.235666\\-0.199659&0.197401&-0.496226&1.1034&-1.06153&0.0130745&-0.402581&0.788863&0.591947&-0.355693\\0.0202145&-0.0043102&-1.39167&-0.873447&0.253648&1.59332&0.988392&-1.32271&-0.951241&-0.192874\\0.238206&1.43934&-1.08785&-0.807931&0.285027&-2.64908&-0.140972&-0.215606&-0.760412&1.77824\\-0.17931&-1.31525&-0.0348558&-0.550808&0.460217&-1.5188&-1.45967&0.508893&-0.942668&1.1059\\\end{matrix}}\right]$

${\hat {x}}=arg{\underset {s}{min}}{\frac {1}{a}}\|y-Ax\|_{2}^{2}+\lambda \|x\|_{1}$

$a=2$ $\lambda =10^{-1}$

sub module iteration 1

calculate $b=x+{\frac {1}{a}}*A^{T}(Ax-y)$

$b=[-1.1194851676858286,-0.5425433278394107,1.1296652382901562,3.3691735377460557,-1.7975322402864389,0.812781461654327,-1.5305641112291442,2.4886377108681788,2.7744632017898994,-2.450112036057253]$

$x^{(1)}={\begin{cases}b^{1}&if\ |b^{(1)}|>\lambda /2a\\0&if\ |b^{(1)}|\leq \lambda /2a\end{cases}}$

$x=[-1.1194851676858286,0.0,1.1296652382901562,3.3691735377460557,-1.7975322402864389,0.812781461654327,-1.5305641112291442,2.4886377108681788,2.7744632017898994,-2.450112036057253]$

Calculate the Error $\|A*x-y\|_{2}=24.020025581623976$

Calculate the new y

$y=Ax=[-7.325091157735414,10.392457780350403,-10.669635806747262,-13.670535847170495,-5.580526608842023]$

Check if error is less then $10^{-5}$

sub module iteration 2

calculate $b=x+{\frac {1}{a}}*A^{T}(Ax-y)$

$b=[-1.1194851676858286,0.0,1.1296652382901562,3.3691735377460557,-1.7975322402864389,0.812781461654327,-1.5305641112291442,2.4886377108681788,2.7744632017898994,-2.450112036057253]$

$x^{(2)}={\begin{cases}b^{2}&if\ |b^{(2)}|>\lambda /2a\\0&if\ |b^{(2)}|\leq \lambda /2a\end{cases}}$

$x=[-1.1194851676858286,0.0,1.1296652382901562,3.3691735377460557,-1.7975322402864389,0.812781461654327,-1.5305641112291442,2.4886377108681788,2.7744632017898994,-2.450112036057253]$

Calculate the Error $\|A*x-y\|_{2}=0.0$

Calculate the new y

$y=Ax=[-7.325091157735414,10.392457780350403,-10.669635806747262,-13.670535847170495,-5.580526608842023]$

Check if error is less then $10^{-5}$

Applications and Motivations

Low-Rank Matrices

The Netflix Prize was accompanied by low-rank matrix recovery or the matrix completion problem.  The approach then fills in the missing values in the user's ratings for movies that the user hasn't seen. These estimates are based on ratings from other users, who have similar ratings if a matrix is created with all the users as rows and the movie titles as columns. Because some users' interests will be similar and therefore overlap, it is possible to reduce the degrees of freedom significantly. This low-rank structure is frequently assumed for the problem domain of collaborative filtering ~cited citations.

Dictionary Learning

The goal in dictionary learning is to infer the original dictionary as possible. Instead of using a predefined dictionary, researchers have found that learning the dictionary by obtaining "dynamic features" from training data often yields representation. Biometric features can be taken from video clips of each subject in a dataset and used to populate the dictionary's columns. Using random projections and sparse representations for iris detection for noncontact biometrics-based authentication systems from video samples has been proposed ~cited citations.

Single-pixel cameras

Single-pixel cameras or single detector imaging are used in situations when detectors are either prohibitively expensive or difficult to miniaturize. A microarray is made up of a large number of miniature mirrors that can be individually turned on and off. The mechanism behind the random sampling, which results in low coherence between measurements, is the most important component of the single-pixel camera. This microarray reflects the light from the scene, and a lens combines all of the reflected beams into one sensor, which is the single detector of the camera used to capture the image ~cited citations.