Sparse Solutions of an Undetermined Linear System

TL;DR

This paper explores methods for finding sparse solutions to underdetermined linear systems, focusing on compressive sensing techniques using l_1 and l_q approaches, and proposes an algorithm based on restricted Isometry to ensure solution uniqueness.

Contribution

It introduces a new algorithm leveraging restricted Isometry for solving underdetermined systems to achieve sparse, unique solutions.

Findings

The l_1 and l_q approaches effectively recover sparse solutions.

The proposed algorithm improves solution uniqueness under restricted Isometry conditions.

Application potential in compressive sensing and signal processing.

Abstract

This work proposes a research problem of finding sparse solution of undetermined Linear system with some applications. Two approaches how to solve the compressive sensing problem: using l_1 approach , the l_q approach with 0 < q < 1. Compressive sensing algorithms are designed to cope with ambiguities introduced by undersampling. We propose an algorithm for restricted Isometry and how it can be used to constrain the undetermined linear system to eventually get a unique solution.

Equations

$$\begin{array}{rl}{\text{min}}_{\mathbf{x}}& \Vert \mathbf{x}{\Vert }_0\\ \text{subject to}& y_i={\mathbf{b}}_i^H\mathbf{x},\mathbf{y}\in {\mathbb{C}}^N{\mathbf{b}}_i\in {\mathbb{C}}^n\mathbf{x}\in {\mathbb{C}}^n\\ & i=1,\dots ,N\end{array}$$

$$\begin{array}{rl}{\text{min}}_{\mathbf{x}}& \Vert \mathbf{x}{\Vert }_1\\ \text{subject to}& y_i={\mathbf{b}}_i^H\mathbf{x},\mathbf{y}\in {\mathbb{C}}^N{\mathbf{b}}_i\in {\mathbb{C}}^n\mathbf{x}\in {\mathbb{C}}^n\\ & i=1,\dots ,N.\end{array}$$

$$\begin{array}{rl}{\text{min}}_{\mathbf{x}}& \Vert \mathbf{x}{\Vert }_0\\ \text{subject to}& y_i=a_i+{\mathbf{b}}_i^H\mathbf{x}+{\mathbf{x}}^H{\mathbf{c}}_i+{\mathbf{x}}^H{\mathbf{Q}}_i\mathbf{x}\\ & i=1,\dots ,N\end{array}$$

$$\begin{array}{rl}{\text{min}}_{\mathbf{x}}& \Vert \mathbf{x}{\Vert }_0\\ \text{subject to}& y_i={\mathbf{b}}_i^H\mathbf{x},\mathbf{y}\in {\mathbb{C}}^N{\mathbf{b}}_i\in {\mathbb{C}}^n\mathbf{x}\in {\mathbb{C}}^n\\ & i=1,\dots ,N\end{array}$$

$$\begin{array}{rl}{\text{min}}_{\mathbf{x}}& \Vert \mathbf{x}{\Vert }_1\\ \text{subject to}& y_i={\mathbf{b}}_i^H\mathbf{x},\mathbf{y}\in {\mathbb{C}}^N{\mathbf{b}}_i\in {\mathbb{C}}^n\mathbf{x}\in {\mathbb{C}}^n\\ & i=1,\dots ,N.\end{array}$$

$$\begin{array}{rl}{\text{min}}_{\mathbf{x}}& \Vert \mathbf{x}{\Vert }_0\\ \text{subject to}& y_i=a_i+{\mathbf{b}}_i^H\mathbf{x}+{\mathbf{x}}^H{\mathbf{c}}_i+{\mathbf{x}}^H{\mathbf{Q}}_i\mathbf{x}\\ & i=1,\dots ,N\end{array}$$

$$y_i=\left[\begin{array}{c}1{\mathbf{x}}^H\end{array}\right]\left[\begin{array}{c}{a}_i{\mathbf{b}}_i^H\\ {\mathbf{c}}_i{\mathbf{Q}}_i\end{array}\right]\left[\begin{array}{c}1\\ \mathbf{x}\end{array}\right]\in \mathbb{C},i=1,\dots ,N.$$

$$\begin{array}{rl}{y}_i=& \text{Tr}(\left[\begin{array}{c}1{\mathbf{x}}^H\end{array}\right]\left[\begin{array}{c}{a}_i{\mathbf{b}}_i^H\\ {\mathbf{c}}_i{\mathbf{Q}}_i\end{array}\right]\left[\begin{array}{c}1\\ \mathbf{x}\end{array}\right])\\ =& \text{Tr}(\left[\begin{array}{c}{a}_i{\mathbf{b}}_i^H\\ {\mathbf{c}}_i{\mathbf{Q}}_i\end{array}\right]\left[\begin{array}{c}1\\ \mathbf{x}\end{array}\right]\left[\begin{array}{c}1{\mathbf{x}}^H\end{array}\right])\\ =& \text{Tr}({\mathbf{\Phi }}_i\mathbf{X})\end{array}$$

$$\begin{array}{rl}{\text{min}}_{\mathbf{X}}& \Vert \mathbf{X}{\Vert }_0\\ \text{subject to}& y_i=\text{Tr}({\mathbf{\Phi }}_i\mathbf{X}),i=1,\dots ,N\\ & \text{rank}(\mathbf{X})=1,{\mathbf{X}}_{1,1}=1,\mathbf{X}⪰0.\end{array}$$

$$\begin{array}{rl}{\text{min}}_{\mathbf{X}}& \text{Tr}(\mathbf{X})+\lambda \Vert \mathbf{X}{\Vert }_1\\ \text{subject to}& y_i=\text{Tr}({\mathbf{\Phi }}_i\mathbf{X}),i=1,\dots ,N\\ & {\mathbf{X}}_{1,1}=1,\mathbf{X}⪰0.\end{array}$$

$$\begin{array}{rl}{\text{min}}_{\mathbf{x}}& \Vert \mathbf{x}{\Vert }_0\\ \text{subject to}& \mathbf{y}=\Phi \mathbf{x}\end{array}$$

$$\begin{array}{rl}{\text{min}}_{\mathbf{x}}& \Vert \mathbf{x}{\Vert }_1\\ \text{subject to}& \mathbf{y}=\Phi \mathbf{x}\end{array}$$

$$\begin{array}{rl}{\text{min}}_{\mathbf{x}}& \Vert \mathbf{x}{\Vert }_q\\ \text{subject to}& \mathbf{y}=\Phi \mathbf{x}\end{array}$$

$$\begin{array}{rl}{\text{min}}_{\mathbf{x}}& \Vert \mathbf{x}{\Vert }_1\\ \text{subject to}& \mathbf{y}=\Phi \mathbf{x}\end{array}$$

$$B:X\in {\mathbb{C}}^{n\times n}\mapsto {\text{Tr}(\mathbf{\Phi X})}_{1\le i\le N}\in {\mathbb{C}}^N$$

$$\left|\frac{\Vert B(\mathbf{X}){\Vert }_2^2}{\Vert \mathbf{X}{\Vert }_2^2}\right|<ϵ$$

$$\mu (\mathbf{A})={\text{max}}_{i,j}\frac{|{\mathbf{A}}_i^H{\mathbf{A}}_j|}{\Vert {\mathbf{A}}_i{\Vert }_2\Vert {\mathbf{A}}_j{\Vert }_2},i\ne j.$$

$$B:X\in {\mathbb{C}}^{n\times n}\mapsto {\text{Tr}(\mathbf{\Phi X})}_{1\le i\le N}\in {\mathbb{C}}^N$$

$$\left|\frac{\Vert B(\mathbf{X}){\Vert }_2^2}{\Vert \mathbf{X}{\Vert }_2^2}\right|<ϵ$$

$$\mu (\mathbf{A})={\text{max}}_{i,j}\frac{|{\mathbf{A}}_i^H{\mathbf{A}}_j|}{\Vert {\mathbf{A}}_i{\Vert }_2\Vert {\mathbf{A}}_j{\Vert }_2},i\ne j.$$

$$\text{min}\Vert \mathbf{y}-\mathbf{Bx}{\Vert }_2^2+\lambda \Vert \mathbf{x}{\Vert }_1.$$

Full text

I System Model and Background

Compressive sensing (CS) aims to solve the following optimization problem to estimate the sparsest vector $\mathbf{x}$

[TABLE]

where the $\|\mathbf{x}\|_{0}$ counts the number of nonzeros in $\mathbf{x}$. This problem is known to be NP-hard, and the theory of compressive sensing has found many approximations to solve this problem. One such approximation is known as basis pursuit (BP) and solves the following convex approximation to (14)

[TABLE]

This approach allows one to solve for the sparsest vector $\mathbf{x}$ for a linear model under certain conditions, which have been studied extensively [3]. Recently, there has been an interest to generalize CS to nonlinear models. In particular in this report we will focus on solving

[TABLE]

where, $\mathbf{y}\in\mathbb{C}^{N},\;\mathbf{a}_{i}\in\mathbb{C}^{n},\;\mathbf{b}_{i}\in\mathbb{C}^{n},\;\mathbf{c}_{i}\in\mathbb{C}^{n},\;\mathbf{Q}\in\mathbb{C}^{n\times n}\;\text{and }\mathbf{x}\in\mathbb{R}^{n}.$ Solving (6) will allow us to solve CS problems where the model takes a quadratic form rather than a linear form.

II Motivation

In this section we give reasons why we want to solve the sparse solution of undetermined systems of linear equations.

II-A Signal and image Compression

This is the most direct application. Lets consider a equation, $y=Hx$ in which $y$ is an under sampled version of image $x$ in a fourier domain. if we no before hand that $x$ is sparse Image, we can recover $x$ from $y$ using equation that provides an approximation of $x$ from $y$.

II-B Compressed Sensing

Compressive sensing (CS) aims to solve the following optimization problem to estimate the sparsest vector $\mathbf{x}$

[TABLE]

where the $\|\mathbf{x}\|_{0}$ counts the number of nonzeros in $\mathbf{x}$. This problem is known to be NP-hard, and the theory of compressive sensing has found many approximations to solve this problem. One such approximation is known as basis pursuit (BP) and solves the following convex approximation to (14)

[TABLE]

This approach allows one to solve for the sparsest vector $\mathbf{x}$ for a linear model under certain conditions, which have been studied extensively [3]. Recently, there has been an interest to generalize CS to nonlinear models. In particular in this report we will focus on solving

[TABLE]

where, $\mathbf{y}\in\mathbb{C}^{N},\;\mathbf{a}_{i}\in\mathbb{C}^{n},\;\mathbf{b}_{i}\in\mathbb{C}^{n},\;\mathbf{c}_{i}\in\mathbb{C}^{n},\;\mathbf{Q}\in\mathbb{C}^{n\times n}\;\text{and }\mathbf{x}\in\mathbb{R}^{n}.$ Solving (6) will allow us to solve CS problems where the model takes a quadratic form rather than a linear form.

II-C Error Correction Code

Let $z$ be a vector which has encoded $x$ by a linear system. A of size $m\times n$ where $m>n$. That is $z=Ax$. Now if the the $z$ obtained is corrupted by noise such that $z=Ax+e$. In such a case we need to obtain the error corrupting $z$, lets consider a matrix $B$ which is $n\times m$ such that $BA=0$ hence we obtain $Bz=Be$, we consider $y=Bz$ therefore we get a linear equation $y=Be$ we can use (1) to solve for $e$ if $y$ is constant.

III Problem Formulation

Similarly to (14) the optimization problem in (6) is not convex due to the $\|\mathbf{x}\|_{0}$ objective function. Therefore first we introduce a convex relaxation to (6). First note that we can rewrite the (6) into a general quadratic form

[TABLE]

Since $y_{i}$ is a scalar using the trace operator we have

[TABLE]

where we define $\mathbf{\Phi}_{i}=\begin{bmatrix}a_{i}\;\;\;\mathbf{b}^{H}_{i}\\ \mathbf{c}_{i}\;\;\;\mathbf{Q}_{i}\end{bmatrix}$ and $\mathbf{X}=\begin{bmatrix}1\\ \mathbf{x}\end{bmatrix}\begin{bmatrix}1\;\;\;\mathbf{x}^{H}\end{bmatrix}$. By definition $\mathbf{X}$ is a Hermitian matrix and it satisfies the constraints $\mathbf{X}_{1,1}=1$ and $\text{rank}(\mathbf{X})=1$. The optimization problem (6) can then be rewritten as

[TABLE]

Recasting the optimization problem to solve for the matrix $\mathbf{X}$ in this form is known as matrix lifting [5], and it was shown that solving for $\mathbf{X}$ one can obtain $\mathbf{x}$ by the rank 1 decomposition of $\mathbf{X}$ by the singular value decomposition (SVD). The problem (9) is still not convex, and the convex approximation is given by

[TABLE]

Here, the $\|\mathbf{X}\|_{1}$ denotes the element-wise $\ell_{1}$-norm, or the sum of magnitudes of all elements in $\mathbf{X}$, and $\lambda$ is a design parameter. The trace of $\mathbf{X}$ is known to be the convex surrogate of the low rank constraint, and the $\|\mathbf{X}\|_{1}$ is the convex surrogate of $\|\mathbf{X}\|_{0}$. The problem (10) is referred to as quadratic basis pursuit (QBP) and the remainder of this report will be based on solving this optimization problem.

IV Restricted Isometry Property

We provide a solution of an undetermined Linear System. consider a linear problem $y=\Phi x$ in which $\Phi$ is a a $m\times N$ where $m$ is the number of measurement and $m<N$, in which. Here $x$ will have many solution but we require the sparsest solution. The formulation is

[TABLE]

The solution thus obtained will be the sparse solution for the linear equation $y=\Phi x$ but the above optimization problem will be a NP hard problem.

The method described above is very computationally expensive. Eq(1) is a NP-hard problem and can not be solved optimization algorithm hence we can construct other ways of obtaining sparse solution like following problems:

[TABLE]

or

[TABLE]

IV-A Restricted Isometry Property

Consider a Linear equation, $\mathbf{y}=\Phi\mathbf{x}=\Phi\mathbf{x}_{0}$ now using equation (1) to obtain $\mathbf{x}$ from $\mathbf{y}$ that will be equal to or close to $x_{0}$. The condition for the equality is derived using the restricted Isometry property. Further restricted Isometry property of $\Phi$ prevents any Amplification of Noise. A vector is said to be s-sparse if it has at most $s$ nonzero entries. Consider now equation (1), we will explain how the restricted Isometry is applied to (1) to construct a constrain on $\Phi$ in order to maintain the uniqueness of the solution. From (6) we know that $\mathbf{y}=\Phi\mathbf{x}=\Phi\mathbf{x}_{0}$, and we know that $\mathbf{x}$ can be derived from $\mathbf{y}$ using equation (1), lets consider that $x*$ is the solution for (1) along with $x_{0}$.

V The $l_{1}$ approach

Although the problem in (1) needs a non polynomial time to solve, in general it can be much more effectively solved using $l_{1}$ minimisation approach. Let us review this approach in this approach in the section. The $l_{1}$ minimization problem is the following,

[TABLE]

Since the $l_{1}$ minimization is a convex problem and can be converted easily into linear programming problem it can be easily solved. Now we have to see under what condition does the solution of (1) is same as the solution of (11). we can define the condition by too Concept: 1) Mutual Coherence, 2) Restricted Isometric property. In this paper we will provide solution using Restricted Isometric property.

V-A The Exact recovery condition for $l_{1}$ norm

In this section some, theoretical results are presented on uniqueness and recoverability of QBP in the noiseless case. Proofs are omitted in this report but can be found in [1]. For convenience we first introduce a linear operator $B:$

[TABLE]

using (18) we can define the restricted isometry property.

Definition 1** (RIP).**

A linear operator $B(\cdot)$ is ($\epsilon,k$)-RIP if

[TABLE]

for all $\|\mathbf{X}\|_{0}\leq k$ and $\mathbf{X}\neq 0$.

From [1] it we have the following theorem for uniqueness,

Theorem 1** (Uniqueness).**

$\overline{\mathbf{x}}$* is a solution to (6) and if $\mathbf{X}^{*}$ satisfies $y=B(\mathbf{X}^{*})$, $\mathbf{X}^{*}\succeq 0$, $\text{rank}(\mathbf{X}^{*})=1$, $\mathbf{X}^{*}_{1,1}=1$ and if $B(\cdot)$ is a $(\epsilon,2\|\mathbf{X}^{*}\|_{0})$-RIP with $\epsilon<1$ then the rank-1 decomposition of $\mathbf{X}^{*}$ is equal to $[1\;\;\;\overline{x}]^{T}$.*

The next concern is on recoverability of the QBP algorithm, first we define the mutual coherence of a dictionary matrix $\mathbf{A}$

Definition 2** (Mutual Coherence).**

For a matrix $\mathbf{A}$ the mutual coherence $\mu(\mathbf{A})$ is given by

[TABLE]

Also define a matrix $\mathbf{D}$ be a matrix satisfying $\mathbf{y}=\mathbf{D}\mathbf{X}^{s}=B(\mathbf{X})$ where $\mathbf{X}^{s}=\text{vec}(\mathbf{X})$. Then the rank-1 decomposition of the solution of QBP is equal to $[1\;\;\;\overline{\mathbf{x}}]^{T}$ if $\|\mathbf{X}\|_{0}<0.5(1+1/\mu(D)).$

V-B Stability

In the previous section we proved condition for exact recovery for $l_{1}$ norm now we move further and prove the condition for stability. When we apply an undetermined Linear transform $H$ to a signal $x_{0}$ there could be some error induced in the undetermined linear system or x0 could have noise in it . In such cases we need to guarantee that the x calculated from measurement y will not blow up in comparison to x0 For this there are two set of theorem defining the upper limit of error between original signal and the recovered signal. we consider that y is contaminated and incomplete observation $y=Hx_{0}+e$, e is the error term .

VI The $l_{q}$ approach

In this section some, theoretical results are presented on uniqueness and recoverability of QBP in the noiseless case. Proofs are omitted in this report but can be found in [1]. For convenience we first introduce a linear operator $B:$

[TABLE]

using (18) we can define the restricted isometry property.

[TABLE]

for all $\|\mathbf{X}\|_{0}\leq k$ and $\mathbf{X}\neq 0$.

A linear operator $B(\cdot)$ is ($\epsilon,k$)-RIP if

[TABLE]

for all $\|\mathbf{X}\|_{0}\leq k$ and $\mathbf{X}\neq 0$.

From [1] it we have the following theorem for uniqueness,

Theorem 2** (Uniqueness).**

$\overline{\mathbf{x}}$* is a solution to (6) and if $\mathbf{X}^{*}$ satisfies $y=B(\mathbf{X}^{*})$, $\mathbf{X}^{*}\succeq 0$, $\text{rank}(\mathbf{X}^{*})=1$, $\mathbf{X}^{*}_{1,1}=1$ and if $B(\cdot)$ is a $(\epsilon,2\|\mathbf{X}^{*}\|_{0})$-RIP with $\epsilon<1$ then the rank-1 decomposition of $\mathbf{X}^{*}$ is equal to $[1\;\;\;\overline{x}]^{T}$.*

The next concern is on recoverability of the QBP algorithm, first we define the mutual coherence of a dictionary matrix $\mathbf{A}$

Definition 3** (Mutual Coherence).**

For a matrix $\mathbf{A}$ the mutual coherence $\mu(\mathbf{A})$ is given by

[TABLE]

Also define a matrix $\mathbf{D}$ be a matrix satisfying $\mathbf{y}=\mathbf{D}\mathbf{X}^{s}=B(\mathbf{X})$ where $\mathbf{X}^{s}=\text{vec}(\mathbf{X})$. Then the rank-1 decomposition of the solution of QBP is equal to $[1\;\;\;\overline{\mathbf{x}}]^{T}$ if $\|\mathbf{X}\|_{0}<0.5(1+1/\mu(D)).$

VI-A Stability

In the previous section we proved condition for exact recovery for $l_{1}$ norm now we move further and prove the condition for stability. When we apply an undetermined Linear transform $H$ to a signal $x_{0}$ there could be some error induced in the undetermined linear system or x0 could have noise in it . In such cases we need to guarantee that the x calculated from measurement y will not blow up in comparison to x0 For this there are two set of theorem defining the upper limit of error between original signal and the recovered signal. we consider that y is contaminated and incomplete observation $y=Hx_{0}+e$, e is the error term .

VII Algorithm

In this section we demonstrate the use of QBP. We draw $\mathbf{B}\in\mathbb{C}^{N\times n}$, and $\mathbf{Q}_{i}\in\mathbb{C}^{n\times n},\;i=1,\ldots N$ from a complex Gaussian random variable. We set $\mathbf{C}=0$ and $\mathbf{a}=0$, $N=18$, $n=20$, and $\mathbf{x}$ is real and shown in Fig. 1. All scenarios are the noiseless case.

Fig. 2 shows the output of using the basis pursuit de-noising algorithm which takes a similar form as (5). The optimization problem that was solved instead of (5) was,

[TABLE]

This optimization problem was solved rather than (5) because it was found that solving (5) often led to the result that the problem is infeasible. The regularization parameter $\lambda$ was chosen to be 50. As expected since basis pursuit algorithms can only account for the linear portion of the model it performs badly and the estimate is far from the truth.

Fig. 3 shows the result of the QBP with $\lambda=50$, the result is a perfect reconstruction and QBP can account for both the quadratic and the linear terms.

Fig. 4 shows the result of the QBP without applying a regularization parameter or setting $\lambda=0$, this is the same as not enforcing the sparsity constraint. It is seen that without the sparsity constraint we recover a very dense estimate that and many solutions can exist.

VIII Results

IX Conclusion

This paper helped understand that geometrically description of a basic Pursuit problem and I can use a similar theory for making a $l_{q}$ minimisation problem have a global minima.