$$
\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:
- Direct Complex Domain Solution
Empoying Non Smooth Convex Optimization methods directly on the Complex Domain problem.
- 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.
Sub 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 $.
Proximal Gradient Method
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
Proximal Gradient Method
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.
Sub Gradient Method
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
- Add Accelerated (FISTA Style) versions of the Sub Gradient Methods and Proximal Gradient Method.
