# Compressive Sensing over the Complex Domain

I'm fairly new to compressive sensing, and I have been looking for a MATLAB implementation of the problem:

$$A x = b$$

where $A \in \mathbb{C}^{m \times n}$, $b \in \mathbb{C}^{m}$ and $x \in \mathbb{C}^{n}$.

Until now I was solving this problem with Least Squares framework:

$$\arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2}$$

Recently I have been having problems with ill conditioned matrices.
Hence I wanted to give ${L}_{1}$ minimization a try:

$$\arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + \lambda {\left\| x \right\|}_{1}$$

Anything you can recommend how to solve such problem?

• $\ell_1$-Magic and SPGL1 are good places to start. Commented Jun 30, 2015 at 21:19
• My package CVX can handle this as well, but you do have to make sure to declare the variable to be convex. It won't be as fast as a dedicated compressed sensing package though. I'm not so sure $\ell_1$-Magic is a good choice; its algorithms are a bit dated. But SPGL1 is great. Commented Jul 2, 2015 at 14:19
• @nabla, Have you looked on my answer?
– Royi
Commented Mar 18, 2020 at 9:17

$$\DeclareMathOperator{\sign}{sign} \DeclareMathOperator{\prox}{prox} \newcommand{\MyParen}[1]{\left( #1 \right)} \newcommand{\MyBrack}[1]{\left\lbrack #1 \right\rbrack} \newcommand{\MyBrace}[1]{\left\lbrace #1 \right\rbrace} \newcommand{\MyNorm}[2]{{\left\| #1 \right\|}_{#2}} \newcommand{\MyAbs}[1]{\left| #1 \right|} \newcommand{\MyNormTwo}[1]{\MyNorm{#1}{2}} \newcommand{\MyCeil}[1]{\lceil #1 \rceil} \newcommand{\MyProd}[2]{\langle #1, #2 \rangle} \newcommand{\MyUndBrace}[2]{\underset{#2}{\underbrace{#1}}} \newcommand{\RR}[1]{\mathds{R}^{#1}}$$

Remark
This is a summary of a work I did for my self out of curiosity for the same question.

# Problem Formulation

The Basis Pursuit Problem:

\begin{align*} \arg \min_{x} \quad & {\left\| x \right\|}_{1} \\ \text{subject to} \quad & A x = b \end{align*}

The ${L}_{1}$ Regularized Least Squares:

\begin{align*} \arg \min_{x} \quad & \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + \lambda {\left\| x \right\|}_{1} \\ \end{align*}

Where $A \in \mathbb{C}^{m \times n}$, $b \in \mathbb{C}^{m}$ and $x \in \mathbb{C}^{n}$.

The ${\ell}_{1}$ norm of a complex vector is defined as:

$$\MyNorm{z}{1} = \sum_{k = 1}^{n} \sqrt{ \MyParen{ \Re \MyParen{ {z}_{k} } }^{2} + \MyParen{ \Im \MyParen{ {z}_{k} } }^{2} }$$

Since the 2 are connected (See Compressive Sensing Results for Unconstrained Form) it enough to deal with one of them.
In this case we will use the second formulation.

# Approaches

There are 2 main approaches for solving this:

1. Direct Complex Domain Solution
Empoying Non Smooth Convex Optimization methods directly on the Complex Domain problem.
2. Transformation into Real Domain Problem
Transforming the proolem into equivalent Real Domain problem and employing optimization methods on it.

I will try showing example of solution of each approach (Though there are many ways to solve this).

## Transformation into Real Domain

By defining $A = B + i C$, $x = y + i z$ and $b = c + i d$ one could rewrite the problem as:

\begin{align*} f \left( x \right) & = \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + \lambda {\left\| x \right\|}_{1} = \frac{1}{2} \MyNormTwo{ \MyParen{ B + i C } \MyParen{y + i z} - c - i d }^{2} + \lambda \MyNorm{y + i z}{1} \\ & = \frac{1}{2} \MyNormTwo{ \MyParen{B y - C z - c} + i \MyParen{C y + B z - d} }^{2} + \lambda \sum_{j = 1}^{n} \sqrt{ {y}_{j}^{2} + {z}_{j}^{2} } \\ & = \frac{1}{2} \MyNormTwo{ B y - C z - c }^{2} + \frac{1}{2} \MyNormTwo{ C y + B z - d }^{2} + \lambda \sum_{j = 1}^{n} \sqrt{ {y}_{j}^{2} + {z}_{j}^{2} } \\ & = \frac{1}{2} \MyNormTwo{ \begin{bmatrix} B & -C \\ C & B \end{bmatrix} \begin{bmatrix} y \\ z \end{bmatrix} - \begin{bmatrix} c \\ d \end{bmatrix} }^{2} + \lambda \sum_{j = 1}^{n} \sqrt{ {y}_{j}^{2} + {z}_{j}^{2} } \\ & = \frac{1}{2} \MyNormTwo{ \begin{bmatrix} B & -C \\ C & B \end{bmatrix} \begin{bmatrix} y \\ z \end{bmatrix} - \begin{bmatrix} c \\ d \end{bmatrix} }^{2} + \lambda \sum_{j = 1}^{n} \MyNormTwo{\begin{bmatrix} \boldsymbol{e}_{j}^{T} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{e}_{j}^{T} \end{bmatrix} \begin{bmatrix} y \\ z \end{bmatrix}} \end{align*}

By setting $\hat{A} = \begin{bmatrix} B & -C \\ C & B \end{bmatrix}$ and $\hat{b} = \begin{bmatrix} c \\ d \end{bmatrix}$ the above can be written as:

$$f \left( y, z \right) = \frac{1}{2} \MyNormTwo{ \hat{A} \begin{bmatrix} y \\ z \end{bmatrix} - \hat{b} }^{2} + \lambda \sum_{j = 1}^{n} \sqrt{ {y}_{j}^{2} + {z}_{j}^{2} }$$

Given the above one could employ many methods to solve the problem such as: * Sub Gradient Method. * Proximal Gradient Method.

The Sub Gradient is given by:

$${\nabla}_{y, z} f \left( y, z \right) = {\hat{A}}^{T} \hat{A} \begin{bmatrix} y \\ z \end{bmatrix} - {\hat{A}}^{T} \hat{b} + \lambda w$$

Where $w \in {\mathbb{R}}^{2 n}$, the Sub Gradient of the last term, is given by:

$${w}_{i} = \begin{cases} \frac{ {y}_{i} }{ \sqrt{{y}_{i}^{2} + {z}_{i}^{2}} } & \text{ if } {y}_{i}^{2} + {z}_{i}^{2} \neq 0, \; i \leq n \\ \frac{ {z}_{i - n} }{ \sqrt{{y}_{i - n}^{2} + {z}_{i - n}^{2}} } & \text{ if } {y}_{i - n}^{2} + {z}_{i - n}^{2} \neq 0, \; > n \\ 0 & \text{ if } {y}_{i} + {z}_{i} = 0 \text{ or } {y}_{i - n} + {z}_{i - n} = 0 \end{cases}$$

Using the above one could apply Sub Gradient Method to find the optimal values of $y$ and $z$ and then compose $x = y + i z$.

Defining $g \MyParen{x} = \sum_{j = 1}^{n} \sqrt{ {x}_{j}^{2} + {x}_{j + n}^{2} }$ where $x \in {\mathbb{R}}^{2 n}$. Then the $\prox$ Operator is defined as:

\begin{align*} \prox_{\lambda g \MyParen{\cdot}} \MyParen{ x } = \arg \min_{u} \MyParen{ \frac{1}{2} \MyNormTwo{u - x}^{2} + \lambda \sum_{j = 1}^{n} \sqrt{ {u}_{j}^{2} + {u}_{j + n}^{2} } } \end{align*}

Using the Sub Gradient defined above the First Optimality Condition becomes:

\begin{align*} 0 \in u - x + \lambda \partial g \MyParen{u} & \Leftrightarrow \begin{cases} {u}_{i} - {x}_{i} + \lambda \frac{ {u}_{i} }{ \sqrt{{u}_{i}^{2} + {u}_{i + n}^{2}} } = 0 & \text{ if } i \leq n \\ {u}_{i} - {x}_{i} + \lambda \frac{ {u}_{i} }{ \sqrt{{u}_{i}^{2} + {u}_{i - n}^{2}} } = 0 & \text{ if } i > n \end{cases} \\ & \Leftrightarrow {u}_{i} = \begin{cases} {x}_{i} \MyParen{1 - \frac{\lambda}{ \sqrt{{x}_{i}^{2} + {x}_{i + n}^{2}} }}_{+} & \text{ if } i \leq n \\ {x}_{i} \MyParen{1 - \frac{\lambda}{ \sqrt{{x}_{i}^{2} + {x}_{i - n}^{2}} }}_{+} & \text{ if } i > n \\ \end{cases} \end{align*}

Where $\MyParen{t}_{+} = \max \MyParen{0, t}$.

Using the above definition one could use Proximal Gradient Method to find the optimal values of $y$ and $z$ and then compose $x = y + i z$.

Remark
The above $\prox$ can actually be decomposed as a $\prox$ of the ${\ell}_{2}$ norm of two elements of the vector $x$ jointly.

## Direct Complex Domain Solution

In order to apply the Proximal Sub Gradient method one should know how to solve the following minimization problem:

$$\tau_{ \lambda } \MyParen{x} = \prox_{ \lambda \MyNorm{\cdot}{1} } \MyParen{x} = \arg \min_{u} \MyBrace{ \MyNorm{u}{1} + \frac{1}{2 \lambda} \MyNormTwo{u - x}^{2} }, \; u, x \in {\mathbb{C}}^{n}$$

Since $\MyNorm{u}{1} = \max_{\MyNorm{p}{\infty} \leq 1} {p}^{H} u, \; p \in {\mathbb{C}}^{n}$ one could write:

\begin{align*} \min_{u} \MyBrace{ \MyNorm{u}{1} + \frac{1}{2 \lambda} {\MyNormTwo{u - x}}^{2} } & = \min_{u} \MyBrace{ \max_{\MyNorm{p}{\infty} \leq 1} {p}^{H} u + \frac{1}{2 \lambda} \MyNormTwo{u - x}^{2} } \\ & = \min_{u} \max_{\MyNorm{p}{\infty} \leq 1} \MyBrace{ {p}^{H} u + \frac{1}{2 \lambda} \MyNormTwo{u - x}^{2} } \\ & = \max_{\MyNorm{p}{\infty} \leq 1} \min_{u} \MyBrace{ {p}^{H} u + \frac{1}{2 \lambda} \MyNormTwo{u - x}^{2} } \end{align*}

Since The function is Convex with respect to $u$ and concave with respect to $p$ and Strong Duality holds one could change the order of minimization and maximization and extract optimal $u$ from the solution of $p$.

The problem $\min_{u} \MyBrace{ {p}^{H} u + \frac{1}{2 \lambda} \MyNormTwo{u - x}^{2} }$ is easily solved by $u = x - \lambda p$. Plugging the result yields:

$$\max_{\MyNorm{p}{\infty} \leq 1} \MyBrace{ {p}^{H} x - \frac{\lambda}{2} \MyNormTwo{p}^{2} } \Rightarrow {p}^{\star}_{k} = \sign \MyParen{{x}_{k}} \min \MyParen{1, \frac{\MyAbs{{x}_{k}}}{\lambda}}$$

The problem above is separable hence the solution is given element wise. Using the Strong Duality:

$${u}^{\star}_{k} = {x}_{k} - \lambda {p}^{\star}_{k} = {x}_{k} - \sign \MyParen{{x}_{k}} \min \MyParen{1, \frac{\MyAbs{{x}_{k}}}{\lambda}} = \sign \MyParen{{x}_{k}} \max \MyParen{0, \MyAbs{{x}_{k}} - \lambda}$$

Where the $\sign \MyParen{\cdot}$ above is the Complex Sign Function $\sign \MyParen{z} = \frac{z}{\MyAbs{z}}$.

Remark
The result is actually the same as the result for the problem in Real Domain with the adaption of the $\sign$ function to the complex domain.

Using the identity:

$$x - u \in \partial f \left( u \right), \; u = \prox_{f \MyParen{\cdot}} \MyParen{x}$$

One could see that:

$$w \in \partial \MyNorm{z}{1}, \; {w}_{k} = \begin{cases} \frac{ {w}_{i} }{ \left | {w}_{i} \right | } & \text{ if } {w}_{i} \neq 0 \\ 0 & \text{ if } {w}_{i} = 0 \end{cases}$$

Hence one could utilize the Sub Gradient Method for this case.

Remark
This is basically deriving the $\sign$ operator for the Complex Domain as the Sub Gradient of the $\MyNorm{\cdot}{1}$ over the Complex Domain.

### Alternating Direction Method of Multipliers (ADMM)

The ADMM Objective Function is given by:

$$\arg \min_{x} \left\{ \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + \lambda {\left\| y \right\|}_{1} \; : \; x - y = 0 \right\}$$

The Augmented Lagrangian is given by:

$$\mathcal{L}_{\rho} \MyParen{x, y, z} = \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + \lambda {\left\| y \right\|}_{1} + {z}^{H} \MyParen{x - y} + \frac{\rho}{2} {\left\| x - y \right\|}_{2}^{2}$$

The algorithm is given by:

\begin{align*} {x}^{k + 1} & = \arg \min_{ x } {\mathcal{L}}_{\rho} \MyParen{ x, {y}^{k}, {z}^{k} } = {\MyParen{{A}^{H} A + \rho I}}^{-1} \MyParen{{A}^{H} b + \rho {y}^{k} - {z}^{k}} \\ {y}^{k + 1} & = \arg \min_{ y } {\mathcal{L}}_{\rho} \MyParen{ {x}^{k + 1}, y, {z}^{k} } = {\tau}_{\frac{\lambda}{\rho}} \MyParen{ \frac{ {z}^{k} }{\rho} + {x}^{k + 1} } \\ {z}^{k + 1} & = {z}^{k} + \rho \MyParen{{x}^{k} - {y}^{k}} \end{align*}

### Fixed Point Iteration

Let ${x}^{\star} \in {\mathbb{C}^{n}}$ be optimal solution of Equation. Then by First Order Optimality Condition:

$$0 \in {A}^{H} \MyParen{A {x}^{\star} - b } + \lambda w, \; w \in \partial \MyNorm{{x}^{\star}}{1}$$

As seen above, $w$ is given by:

$${X}^{-1} x, \; X = \begin{bmatrix} \left | {x}_{1} \right | & & & \\ & \left | {x}_{2} \right | & & \\ & & \ddots & \\ & & & \left | {x}_{n} \right | \end{bmatrix} = \operatorname{diag} \MyParen{ \MyAbs{x} }$$

Hence the First Order Optimality Condition can be written as:

$$0 \in {A}^{H} \MyParen{A {x}^{\star} - b } + \lambda w, \; w \in \partial \MyNorm{{x}^{\star}}{1} \Rightarrow {x}^{\star} = \MyParen{{A}^{H} A + \lambda {X}^{-1}}^{-1} {A}^{H} b$$

Hence the Fixed Point Iteration can be written as:

$${x}^{k + 1} = \MyParen{{A}^{H} A + \lambda {\MyParen{{X}^{k}}}^{-1}}^{-1} {A}^{H} b, \; {X}^{k} = \operatorname{diag} \MyParen{ \MyAbs{{x}^{k}} }$$

Yet, one could notice that:

$$\MyParen{ {A}^{H} A + \lambda {X}^{-1} } X {A}^{H} = {A}^{H} \MyParen{ \lambda I + A X {A}^{H} }$$

If both $\MyParen{ {A}^{H} A + \lambda {X}^{-1} }$ and $\MyParen{ \lambda I + A X {A}^{H} }$ invertible then:

$$X {A}^{H} \MyParen{ \lambda I + A X {A}^{H} }^{-1} = \MyParen{ {A}^{H} A + \lambda {X}^{-1} }^{-1} {A}^{H}$$

Which allows inversion of (Sometimes) much smaller matrix:

$${x}^{k + 1} = {X}^{k} {A}^{H} \MyParen{A {X}^{k} {A}^{H} + \lambda I}^{-1} b$$

Moreover, in cases where ${A}^{H} A + \lambda {X}^{-1}$ or $\lambda I + A X {A}^{H}$ is a Toeplitz Matrix (For instance when $A$ is the Fourier Basis) it can be solved efficiently by Levinson Recursion.

Remark
The first form of the Fixed Point iteration can be seen as Iteratively Reweighted Least Squares (IRLS).

# Simulation Results

All methods above were implemented in MATLAB.
The data was built as following:

%% Simulation Parameters

numRows = 64;
numCols = 256;

paramLambda = 4.0;

numIterations = 1000;

methodIdx = 0;

%% Generate Data

mA = randn([numRows, numCols]) + (1i * randn([numRows, numCols]));
vB = randn([numRows, 1]) + (1i * randn([numRows, 1]));

hCalcErrorNorm = @(mX, vXRef) sum((mX - vXRef) .* conj(mX - vXRef));
hCalcObjFunVal = @(mX) 0.5 * sum(((mA * mX) - vB) .* conj((mA * mX) - vB)) + paramLambda * sum(abs(mX), 1);


The reference was created using CVX.

Those are the results:

As can be seen above, the IRLS method is the fastest to converge and since it can do that in few iterations it is the most efficient (As the Run Time if for 1000 iterations).

One could have a look at Basis pursuit denoising (BPDN) and LASSO with a given measurement error? in order to utilize the methods above to solve the Basis Pursuit formulation of the problem.

The full code, including validation using CVX, can be found in my StackExchange Mathematics Q1344369 GitHub Repository.

To Do