# Reconstruction of jointly sparse vectors via manifold optimization

Armenak Petrosyan<EMAIL_ADDRESS>Hoang Tran<EMAIL_ADDRESS>Clayton
Webster<EMAIL_ADDRESS>Department of Computational and Applied
Mathematics, Oak Ridge National Laboratory, 1 Bethel Valley Road, P.O. Box
2008, Oak Ridge TN 37831-6164, USA  Department of Mathematics, University of
Tennessee, 1403 Circle Drive, Knoxville, TN 37916, USA

###### Abstract

In this paper, we consider the challenge of reconstructing jointly sparse
vectors from linear measurements. Firstly, we show that by utilizing the rank
of the output data matrix we can reduce the problem to a full column rank
case. This result reveals a reduction in the computational complexity of the
original problem and enables a simple implementation of joint sparse recovery
algorithms for full-rank setting. Secondly, we propose a new method for joint
sparse recovery in the form of a non-convex optimization problem on a non-
compact Stiefel manifold. In our numerical experiments our method outperforms
the commonly used $\ell_{2,1}$ minimization in the sense that much fewer
measurements are required for accurate sparse reconstructions. We postulate
this approach possesses the desirable rank aware property, that is, being able
to take advantage of the rank of the unknown matrix to improve the recovery.

###### keywords:

joint sparsity, manifold optimization , non-compact Stiefel manifold , non-
convex optimization , Huber regularization

###### MSC:

[2010] 65F99, 68P30 , 94A12

##  1 Introduction

The sparse sampling and recovery problem has been a subject of intensive
research for the last two decades, especially after major advances in the
field of compressed sensing, through the seminal works [5, 11]. In the
classical problem, the goal is to find the vector $x\in\mathbb{R}^{N}$ from
the system $Ax=y$ where $A\in\mathbb{R}^{M\times N}$, under the assumption
that the vector $x$ is sparse, i.e., has only few non-zero entries. Here, $N$
represents the dimensionality of the unknown vector $x$ and $M$ is the number
of linear measurements performed on $x$. A typical difficulty in solving the
problem stems from the fact that one often has only a limited number of
measurements, i.e., $M<N$, so the system is underdetermined.

In the joint sparse recovery problem, also known as multiple measurements
vectors (MMV), simultaneous sparse approximation, or collaborative sparse
coding, we instead solve a matrix equation $AX=Y$ where the unknown matrix
$X\in\mathbb{R}^{N\times K}$ has only few non-zero rows. $X$ can be viewed as
a collection of $K$ unknown vectors in $\mathbb{R}^{N}$ which share the same
sparsity support. In this case, it has been demonstrated [30, 14, 26] that
recovery of multiple vectors simultaneously can take advantage of their joint
sparsity, which helps reduce the number of measurements compared to individual
reconstructing each at a time. The joint sparse recovery problem has
applications in hyper-spectral imaging [7, 13, 35, 37], source localization
[2, 25] and many other fields. In these problems, the joint sparse recovery is
often conducted under additional positivity constraints, which are out of the
scope of this paper. Our proposed method therefore is not applicable directly,
but is an important block towards these applications. Another example is the
CLSUnSAL algorithm [19] for hyperspectral unmixing, which is based on a joint
sparse regularization.

In this paper, we first utilize the rank of the output data matrix $Y$ to
reduce the problem to the case where unknown matrix has full column rank. More
specifically, we prove the following theorem (Section 3).

######  Theorem 1.

Denote by $\Sigma_{s,r}$ the set of all matrices in $\mathbb{R}^{N\times K}$
of rank at least $r$ and having at most $s$ non-zero rows ($s$-row sparse).
Let $A\in\mathbb{R}^{M\times N}$ be such that the induced map $Z\mapsto AZ$ is
injective on $\Sigma_{s,r}$, $X\in\mathbb{R}^{N\times K}$ be $s$-row sparse,
$Y=AX$ and $\textnormal{rank\,}(Y)=r$. Then

  1. 1.

$Y$ can be written as a product $Y=VU$ where $V\in\mathbb{R}^{M\times r}$ with
$\textnormal{rank\,}(V)=r$ and $U\in\mathbb{R}^{r\times K}$ satisfying
$UU^{T}=I_{r}$ where $I_{r}$ is the identity matrix in $\mathbb{R}^{r\times
r}$.

  2. 2.

The problem

$\operatorname*{arg\,min}_{Z\in\mathbb{R}^{N\times r}}\|Z\|_{0}\ \ \text{ s.t.
}\ \ AZ=V$ (1)

has a unique solution $W\in\mathbb{R}^{N\times r}$. This solution is $s$-row
sparse and has full column rank.

  3. 3.

The original matrix $X$ can be computed by $X=WU.$

Here, $\|Z\|_{0}$ is the number of non-zero rows of matrix $Z$.

This result helps to reduce the computational complexity of the original
problem and allows a simple implementation of joint sparse recovery algorithms
for full-rank setting such as the MUSIC algorithm [30, 14, 15] to the rank
deficient case. We note that a similar reduction approach has been presented
in [26, Section III] and [25]. However, in that work, the authors split the
problem into two parts: first finding the support by solving the reduced
problem, and then finding the unknown $X$ using the knowledge of the support,
whereas we have more complete result utilizing the rank.

As a relaxation of the minimization problem (1) in the theorem above, we
propose a new method for the recovery of jointly sparse vectors in the form of
a non-convex optimization problem on the non-compact Stiefel manifold:

$\operatorname*{arg\,min}_{Z\in\mathbb{R}^{N\times
r},\;\textnormal{rank\,}(Z)=r}\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{2,1}\ \ \text{
s.t. }\ \ AZ=V.$ (2)

For the definition of $\ell_{2,1}$ norm, see (5). Our method generalizes the
$\ell_{1}/\ell_{2}$ minimization method studied in [12, 17, 21, 29, 36]. We
demonstrate in our numerical examples that it outperforms the $\ell_{2,1}$
minimization commonly used to solve the joint sparse recovery problem in the
sense that much fewer measurements are required for accurate sparse
reconstructions.

For numerical examples, we consider the unconstrained version of (2)

$\operatorname*{arg\,min}_{Z\in\mathbb{R}^{N\times
r},\;\textnormal{rank\,}(Z)=r}\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{2,1}+\frac{\lambda}{2}\|AZ-V\|^{2}_{F},$
(3)

where $\|\cdot\|_{F}$ is the Frobenius norm. Note that (3) is the Lagrangian
form of the equivalent problem

$\operatorname*{arg\,min}_{Z\in\mathbb{R}^{N\times
r},\;\textnormal{rank\,}(Z)=r}\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{2,1}\text{ s.t.
}\|AZ-V\|^{2}_{F}\leq\epsilon$

for some other parameter $\epsilon>0$, that measures the noise level in the
measurements. Just like the $\epsilon$, the choices of tuning parameter
$\lambda>0$ are problem dependent, with small values of $\epsilon$
corresponding to large values of $\lambda$, and often determined empirically
(see [9, 23] for more detailed discussion). By choosing large $\lambda$,
ensures equation (3) gives the same solution as (2). Currently, there are
several packages available for manifold optimization (see [18] for a detailed
list). We give preference to the Pymanopt package for its automatic
differentiation capabilities [31]. It is the Python version of the original
Manopt package [4, 34] written for MATLAB. We also smooth out the $\ell_{2,1}$
norm with the Huber function to be able to apply the conjugate gradient
algorithm in Pymanopt. For users of the Manopt package, calculation of
derivatives is needed, and we will explicitly compute the Euclidean gradient
of the objective function in (3) as well as the gradient of its Huber
regularized version.

The paper is organized as follows. For the rest of Section 1, we review the
single sparse and joint sparse recovery problems. Section 2 provides
background results on the rank-awareness of the $\ell_{0}$ minimization
problem for matrices. Section 3 discusses a simple strategy to reduce the rank
deficient problems to the full-rank setting. We relax the $\ell_{0}$
minimization to a manifold optimization problem in Section 4. Section 5
discusses the Huber regularization. Calculation of the derivative of our
objective functions will be presented in Section 6. Finally, numerical example
is given in Section 7.

###  1.1 Sparse vector recovery problem

Many signals of interest have (almost) sparse representation in some basis.
For instance, images typically have sparse representations in wavelet bases,
digital audio signals have sparse representation in Gabor systems, etc. Let
$x\in\mathbb{R}^{N}$ be the coefficients of the signal in such a basis and
$y\in\mathbb{R}^{m}$ be a vector of linear measurements of the signal. Via the
Riesz theorem, the relationship between $x$ and $y$ can be written as $y=Ax$
where the matrix $A$ can be explicitly computed.

In classical sparse sampling and recovery problem, the goal is to find the
vector $x\in\mathbb{R}^{N}$ from the system $Ax=y$ where
$A\in\mathbb{R}^{M\times N}$, under the assumption that the vector $x$ is
sparse. Even though the number $s$ of non-zero entries in $x$ can be much
smaller than the dimension $N$, the indices of these entries are unknown, and
directly solving the system by first finding the non-zero entry indices of $x$
is an NP-hard problem [27]. Other methods have been suggested for solving it,
one of the most popular is the $\ell_{1}$ minimization problem

$\operatorname*{arg\,min}_{z\in\mathbb{R}^{N}}\|z\|_{1}\ \ \text{ s.t. }\ \
Az=y,$

or its unconstrained version

$\operatorname*{arg\,min}_{z\in\mathbb{R}^{N}}\|z\|_{1}+\frac{\lambda}{2}\|Az-y\|^{2}_{2}.$

It has been proved that under some conditions on the matrix $A$, specifically
the null space property (NSP) and the restricted isometry property (RIP),
these methods are able to recover the vector $x$ exactly and stably.

###  1.2 The joint sparse recovery problem

In the joint sparse recovery problem we want to solve a matrix equation of the
form

$AX=Y,\;\;A\in\mathbb{R}^{M\times N},\;\;X\in\mathbb{R}^{N\times K}\text{ and
}Y\in\mathbb{R}^{M\times K},$ (4)

given that $X\in\mathbb{R}^{N\times K}$ is a $s$-row sparse matrix, i.e., if
$X=(X[1],\dots,X[N])^{T}$ then $X[i]$’s are all zero for all indices $i$’s
outside some subset of $\\{1,\dots,N\\}$ of size $s$.

Our goal is to find a method to exactly recover $X$ from (4), where the matrix
$A\in\mathbb{R}^{M\times N}$ is chosen with small $M$. One expects that when
increasing the number $K$ of unknown jointly sparse vectors, a high rank of
$X$ will result in a reduction of $M$ which corresponds to the number of
linear measurements done on the vectors. Define the $\ell_{2,1}$ norm of $X$

$\|X\|_{2,1}=\sum_{n=1}^{N}\|X[n]\|_{2}.$ (5)

Similar to the case of single vector recovery, the following convex
minimization problem has been widely used in place of directly solving for
sparsest solution of (4)

$\operatorname*{arg\,min}_{Z\in\mathbb{R}^{N\times K}}\|Z\|_{2,1}\text{ s.t.
}AZ=Y.$ (6)

In the worst case (uniform recovery) scenario, however, this method does not
allow measurement reduction when the number of vectors is increased. As the
next theorem shows, null space property is a necessary and sufficient
condition for both the single sparse recovery via $\ell_{1}$ minimization and
the joint sparse recovery via $\ell_{2,1}$ minimization. Specifically, if $A$
fails to recover an $s$-sparse vector with $\ell_{1}$ minimization, then it
will fail to recover some $s$-row sparse matrix via (6), no matter how high
the rank of the unknown matrix is. As such, the $\ell_{2,1}$ minimization is
known to be a rank blind algorithm.

######  Theorem 2 ([22, 10]).

Given a matrix $A\in\mathbb{R}^{M\times N}$, every s-row sparse matrix
$X\in\mathbb{R}^{N\times K}$ can be exactly recovered from (6) if and only if
$A$ satisfies the NSP:

$\|x_{I}\|_{1}<\|x_{I^{c}}\|_{1},\ \forall
x\in\ker(A)\setminus\\{\mathbf{0}\\},\ \forall I\subset\\{1,\dots,N\\},\text{
with }|I|\leq s,$

where $x_{I}$ is the vector that we get after nullifying all the entries of
$x$ except the ones with indexes in $I$. Furthermore, if there exists an
$x\in\ker(A)$ such that

$\|x_{I}\|_{1}>\|x_{I^{c}}\|_{1},\;|I|\leq s,$

then, for any $1\leq r\leq\min\\{s,K\\}$, there is an s-row sparse matrix
$X\in\mathbb{R}^{N\times K}$ with $\textnormal{rank\,}(X)=r$ which cannot be
recovered from (6).

There have been several efforts to exploit the rank of the unknown matrix and
design rank aware algorithms to improve the joint sparse recovery results.
Examples of such algorithms include those based on MUSIC [15, 24] and
orthogonal matching pursuit [10]. On the other hand, we are not aware of
available rank aware algorithms based on (convex) functional optimization,
CoSaMP or thresholding approaches.

One may think of resolving the rank blindness issue of $\ell_{2,1}$
minimization by replacing the $\ell_{2,1}$ norm with a different functional,
i.e., consider the problem

$\operatorname*{arg\,min}_{Z\in\mathbb{R}^{N\times K}}\phi(Z)\;\;\text{ s.t.
}AZ=Y$ (7)

where $\phi$ is a convex functional in $\mathbb{R}^{N\times K}$. However, as
we point out below, no such convex functional like $\|\cdot\|_{2,1}$ will be
rank aware, and we may need to change the geometry or look into non-convex
optimization for such desirable property.

First, if (7) has a unique solution for every $Y=AX$ with $X\in\Sigma_{s,r}$
then, for any $r\leq r^{\prime}\leq s$ and $X\in\Sigma_{s,r^{\prime}}$, it has
a unique solution too. Hence (7) is rank aware if it allows unique recovery in
the class $\Sigma_{s,r}$ ($r\leq s$) and fails to recover a matrix $X_{0}$ in
$\Sigma_{s,r^{\prime}}$ for some $r^{\prime}<r$. Observe that, every
$X\in\Sigma_{s,r}$ is the unique solution of (7) if and only if

$\phi(X+W)>\phi(X),\ \ \forall X\in\Sigma_{s,r},\ \forall
W\in(\ker(A))^{K}\setminus\\{\mathbf{0}\\}.$ (8)

And (7) strongly fails to recover an $X_{0}\in\Sigma_{s,r^{\prime}}$ if and
only if there exists a $W_{0}\in(\ker(A))^{K}\setminus\\{\mathbf{0}\\}$ such
that $\phi(X_{0}+W_{0})<\phi(X_{0})$. However, by the continuity of $\phi$ and
the fact that $\Sigma_{s,r}$ is dense in $\Sigma_{s,r^{\prime}}$ when $r\geq
r^{\prime}$, (8) implies that $\phi(X_{0}+W_{0})\geq\phi(X_{0})$, a
contradiction. As a result, (7) cannot recover any matrix in $\Sigma_{s,r}$
while strongly fail to recover in the class $\Sigma_{s,r^{\prime}}$ for
$r^{\prime}\leq r$.

In this paper, we utilize the rank of $Y$ to extract information about $X$ and
reduce the problem to a smaller dimensional problem in full-rank setting. As a
first step towards a rank aware optimization approach for uniform recovery, we
develop a new non-convex optimization method on manifold, which achieves a
significant reduction in the number of measurements in joint sparse recovery,
compared to $\ell_{2,1}$ minimization.

##  2 Preliminary results

This section revisits the $\ell_{0}$ minimization for matrices, based on
results in [6, 10]. We present a uniqueness condition for the sparse recovery
on $\Sigma_{s,r}$, that is, the map $Z\mapsto AZ$ is injective on
$\Sigma_{s,r}$. Before we characterize this condition, let us make a few
observations.

######  Definition 1.

The spark of the matrix $A$, denoted by $\textnormal{spark\,}(A)$, is the
smallest number $j$ such that there exists $j$ linearly dependent columns in
$A$.

Note that any $\textnormal{spark\,}(A)-1$ columns of $A$ are linearly
independent.

######  Proposition 1.

If $A:\Sigma_{s,r}\to\mathbb{R}^{M\times K}$ is injective then for every
$X\in\Sigma_{s}$, $\textnormal{rank\,}(X)=\textnormal{rank\,}(AX).$

###### Proof.

We prove this proposition into two steps:

  1. 1.

if $A:\Sigma_{s,r}\to\mathbb{R}^{M\times K}$ is injective then
$\textnormal{spark\,}(A)>s$,

  2. 2.

if $\,\textnormal{spark\,}(A)>s$ then
$\textnormal{rank\,}(X)=\textnormal{rank\,}(AX)$, $\forall X\in\Sigma_{s}$.

To prove i), assume otherwise ($\textnormal{spark\,}(A)\leq s$), there exists
an $I\subseteq\\{1,\dots,N\\}$, $|I|=s$ and vectors
$v_{1},\dots,v_{r}\in\mathbb{R}^{N}$ which are linearly independent and
supported on $I$ such that $Av_{1}=0$. Take the matrices

$X=(v_{1},\cdots,v_{r},\overbrace{0,\cdots,0}^{K-r}),\;\;\widetilde{X}=(-v_{1},\cdots,v_{r},\overbrace{0,\cdots,0}^{K-r}).$

Then $X,\widetilde{X}\in\Sigma_{s,r}$, $X\neq\widetilde{X}$ and
$AX=A\widetilde{X}$, which contradicts the injectivity of $A$.

For ii), suppose $\textnormal{rank\,}(X)>\textnormal{rank\,}(AX)$ for
$X\in\Sigma_{s}$. Then, there exists a vector $x\neq 0$ in the span of columns
of $X$ such that $Ax=0$. $x$ is $s$-sparse, hence we found a non-trivial
combination of $s$ columns of $A$ that is equal to $0$, contradicting the
assumption that $\textnormal{spark\,}(A)>s$. ∎

The uniqueness condition for the sparse recovery on $\Sigma_{s,r}$ is
characterized in the following theorem.

######  Theorem 3.

The following statements are equivalent:

  1. 1.

$A:\Sigma_{s,r}\to\mathbb{R}^{M\times K}$ is injective,

  2. 2.

$s\leq\dfrac{\textnormal{spark\,}(A)-1+r}{2}$,

  3. 3.

For every $X\in\Sigma_{s,r}$, $X$ is the unique solution of the minimization
problem

$\min_{Z\in\mathbb{R}^{N\times K}}\|Z\|_{0}\text{ s.t. }AZ=Y,$ (9)

where $Y=AX$.

###### Proof.

$(1)\Leftrightarrow(2)$ follows as an intermediate consequence of [10, Theorem
2].

$(1)\Rightarrow(3)$ Suppose $A:\Sigma_{s,r}\to\mathbb{R}^{M\times K}$ is
injective. From Proposition 1,
$\textnormal{rank\,}(Y)=\textnormal{rank\,}(X)\geq r$. Let $\widetilde{X}$ be
a solution of (9), then $\|\widetilde{X}\|_{0}\leq\|X\|_{0}=s$. Also, since
$A\widetilde{X}=Y$,
$\textnormal{rank\,}(\widetilde{X})\geq\textnormal{rank\,}(Y)\geq r$, hence
$\widetilde{X}\in\Sigma_{s,r}$. From the injectivity of $A$, we have
$X=\widetilde{X}$, and $X$ is the unique solution of (9).

$(3)\Rightarrow(1)$ Suppose $X,\widetilde{X}\in\Sigma_{s,r}$ and
$AX=A\widetilde{X}=Y$. Then both of them are the unique solutions of the same
minimization problem, therefore $X=\widetilde{X}$. ∎

Theorem 3 shows that the $\ell_{0}$ minimization (9) is rank aware. Indeed,
matrices $X$ with larger rank allow $A$ with smaller
$\textnormal{spark\,}(A)$, therefore can be reconstructed from fewer
measurements. Alternatively, for a fixed $A$, the upper bound of the sparsity
of unknown matrices that can be exactly recovered grows linearly with its
rank. In general, the $\ell_{0}$ norm as a functional is not suitable for
computations. We present a relaxation of this problem using manifold
optimization approach in Section 4. Immediately below is a useful idea to
reduce matrix equations to full-rank setting.

##  3 Reduction to full-rank setting

We present a simple way to transform the matrix equation to be solved into one
where the unknown is full-rank matrix. At the expense of an additional
decomposition of output data, the benefit of this practice is twofold. First,
we reduce the dimension of the matrix equation and now find solutions in a
potentially much lower dimensional space $\mathbb{R}^{N\times r}$, compared to
the original one in $\mathbb{R}^{N\times K}$. More importantly, this step
allows us to apply and exploit the strength of several joint sparse recovery
methods which are particularly powerful in the full-rank case to the more
general setting with possible rank defect.

Our main result in this section is the following.

######  Theorem 4.

Let $X\in\mathbb{R}^{N\times K}$ be $s$-row sparse, $A\in\mathbb{R}^{M\times
N}$ such that $A:\Sigma_{s,r}\to\mathbb{R}^{M\times K}$ is injective and
$Y=AX$ with $\textnormal{rank\,}(Y)=r$. Then

  1. 1.

$Y$ can be written as a product $Y=VU$ where $V\in\mathbb{R}^{M\times r}$ with
$\textnormal{rank\,}(V)=r$ and $U\in\mathbb{R}^{r\times K}$ satisfying
$UU^{T}=I_{r}$ where $I_{r}$ is the identity matrix in $\mathbb{R}^{r\times
r}$.

  2. 2.

The problem

$\operatorname*{arg\,min}_{Z\in\mathbb{R}^{N\times r}}\|Z\|_{0}\text{ s.t.
}AZ=V$ (10)

has a unique solution $W\in\mathbb{R}^{N\times r}$. This solution is $s$-row
sparse and has full column rank.

  3. 3.

The original matrix $X$ can be computed by $X=WU.$

###### Proof.

Since $\textnormal{rank\,}(X)=r$, we have from Proposition 1 that
$\textnormal{rank\,}(Y)=r$. The decomposition of $Y$ in i) is well-known. We
note that this decomposition is not unique but one such representation can be
computed from the compact singular value decomposition of $Y$ or the RQ
decomposition.

For ii), using the fact that $UU^{T}=I_{r}$, we get $V=YU^{T}=AXU^{T}$. Also,
$XU^{T}$ is $s$-row sparse and
$\textnormal{rank\,}(XU^{T})=\textnormal{rank\,}(X)=r$. From Theorem 3,
$W=XU^{T}$ is the unique solution to (10).

For iii), $AWU=VU=Y$ and $\|WU\|_{0}=\|X\|_{0}$, hence $X=WU$. ∎

######  Remark 3.1.

Other decompositions $Y=VU$ of $Y$ can also be considered. As long as the
columns of $V$ span the same subspace, we will have the same result as above.

Theorem 4 shows that instead of directly solving problem (9), one can solve
the reduced problem (10), where the unknown matrix has full rank. Then, the
solution of original problem can be obtained with one matrix multiplication.
Note that for actual data that contains noise, we can compute the SVD
decomposition of $Y$ and threshold the eigenvalues to get the essential rank
of the matrix. Moreover, if $N$ is too large, and we have an estimate on the
sparsity $s$, we may consider the truncated SVD decomposition instead with
first $s$ largest singular values. This may reduce the computational cost
dramatically when $s\ll N$.

Notice that if $\widetilde{W}$ is an approximate solution of (10), since the
columns of $U$ make up a Parseval frame, we have
$\|WU-\widetilde{W}U\|_{F}=\|W-\widetilde{W}\|_{F}$. Thus the approximation
error in solving the reduced problem does not result in a larger (in
$\ell_{2}$ sense) error for the solution of the original problem.

##  4 Non-convex manifold optimization for joint sparse recovery

In this section, we present a novel relaxation to $\ell_{0}$ minimization in
form of an optimization problem on manifold.

Consider problem (10). As known from the previous section that its solution
must have full rank, we can restrict this problem to the space of all full-
rank matrices in $\mathbb{R}^{N\times r}$. This adds a rank constraint to the
joint sparsity recovery, and our idea is to recast the optimization problem
with full-rank condition in the Euclidean space into one on a manifold that
encodes this constraint (i.e., a non-compact Stiefel manifold), [1].
Furthermore, we relax the $\ell_{0}$ norm of the solution with the concept of
orthogonal factor of matrices.

######  Definition 2.

Let $Z$ be a matrix in $\mathbb{R}^{N\times r}$ with
$\textnormal{rank\,}(Z)=r$, we call $Z(Z^{T}\\!Z)\\!^{-\frac{1}{2}}$ its
orthogonal factor.

The orthogonal factor is the product of orthogonal factors in the compact SVD
decomposition of $Z$. Note that $\|Z\|_{0}=\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{0}$,
it is quite natural to relax $\|Z\|_{0}$ with
$\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{2,1}$. Define the non-compact Stiefel manifold

$\mathbb{S}(N,r)=\\{Z\in\mathbb{R}^{N\times r}:\textnormal{rank\,}(Z)=r\\},$

we seek to approximate the solution to (10) via the following manifold
optimization problem

$\operatorname*{arg\,min}_{Z\in\,\mathbb{S}(N,r)}\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{2,1}\
\ \text{ s.t. }\ \ AZ=V.$ (11)

######  Remark 4.1.

  1. 1.

$\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{2,1}$ can also be interpreted as follows: if we
take any orthonormal basis for the range of $Z$, then the $\ell_{2,1}$ norm of
that basis is invariant of the choice of the basis and exactly equal to
$\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{2,1}$.

  2. 2.

$\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{2,1}$ is a bounded functional and therefore not
geodesically convex. A bounded functional is convex on a Riemannian manifold
if and only if it is constant ([33, Corollary 2.5]).

  3. 3.

Since proposed manifold optimization method is non-convex, the outcomes of
gradient based optimization methods for solving it largely depend on the
quality of the initialization. As such, the convex $\ell_{2,1}$ minimization
can be used to obtain a good initial value for the manifold method, which will
be then used to arrive at a better reconstruction using the conjectured rank-
awareness property of the manifold method.

###  4.1 Particular cases

Some particular cases of (11) are of interest on their own. For $r=s$, there
holds $\|W(W^{T}W)^{-\frac{1}{2}}\|_{2,1}=s$ with $W$ being the true solution
to (11). Indeed, from Remark 4.1, $\|W(W^{T}W)^{-\frac{1}{2}}\|_{2,1}$ is the
$\ell_{2,1}$ norm of $s$ orthonormal basis vectors that have joint sparsity
$s$. Since each of the non-zero rows of this matrix forms a unit vector in
$\ell_{2}$ norm, its $\ell_{2,1}$ norm is equal to $s$.

For $r=1$, (11) is equivalent to minimizing the $\ell_{1}/\ell_{2}$
functional:

$\operatorname*{arg\,min}_{z\in\mathbb{R}^{N}}\|z\|_{1}/\|z\|_{2}\ \ \text{
s.t. }\ \ Az=v.$

This functional has been shown to enhance sparsity of the solutions and
outperform $\ell_{1}$ on many test problems [17, 21, 12, 36]. However,
rigorous theory on whether $\ell_{1}/\ell_{2}$ is superior to $\ell_{1}$ in
uniform recovery, specifically, whether the NSP is sufficient for the exact
recovery using $\ell_{1}/\ell_{2}$ is largely open. One particular challenge
is that $\ell_{1}/\ell_{2}$ penalty is neither convex nor concave nor
separable, therefore could not be treated under current analysis, see [32].
Extending to the general case ($r>1$), we conjecture that the NSP for the
manifold optimization problem (11) is less demanding than that for
$\ell_{2,1}$ minimization (which is identical to standard NSP, see Theorem 2).
The theoretical verification of this postulation is a subject of ongoing
research.

We remark that other functionals are possible. In this paper, we mainly
demonstrate the effective performance of the manifold optimization method,
leaving open the question of optimality of the proposed functionals.

##  5 Huber regularization

We consider the following unconstrained problem, which produces the same
solution as (11) for large $\lambda$:

$\operatorname*{arg\,min}_{Z\in\,\mathbb{S}(N,r)}\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{2,1}+\frac{\lambda}{2}\|AZ-V\|^{2}_{F}.$
(12)

Satisfactory choices of $\lambda$ are problem dependent and often determined
empirically, as have been done herein. With notice that $\ell_{2,1}$ norm is
non-smooth, we further relax $\|\cdot\|_{2,1}$ with Huber regularization to be
able to apply manifold optimization methods, which have been mostly developed
for smooth optimization. In particular, we replace the $\|X[n]\|_{2}$ with
$H_{\delta}(\|X[n]\|_{2})$ where $H_{\delta}\in C^{1}$ is the Huber function

$H_{\delta}(x)=\begin{cases}x-\delta/2,&x\geq\delta,\\\
\frac{x^{2}}{2\delta},&0\leq x<\delta,\end{cases}$

(see Figure 1) and consider the problem

$\operatorname*{arg\,min}_{Z\in\,\mathbb{S}(N,r)}F_{\delta}(Z)+\frac{\lambda}{2}\|AZ-V\|^{2}_{F},$
(13)

where

$F_{\delta}(Z)=\sum_{n=1}^{N}\left[H_{\delta}(\|Z(Z^{T}Z)^{-\frac{1}{2}}[n]\|_{2})+\frac{\delta}{2}\right].$

The following proposition reveals that (13) is an effective surrogate for
problem (12).

######  Proposition 2.

If (12) has a unique solution and there exists $\delta_{0}>0$ such that (13)
has unique solution for all $\delta<\delta_{0}$, then the solution to (13)
converges to the solution to (12) as $\delta$ converges to $0$.

###### Proof.

$\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{2,1}$ is a continuous functional on the compact
Stiefel manifold and $F_{\delta}(Z)$ is decreasing with respect to $\delta$
and pointwise converges to $\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{2,1}$. Thus, from
[8, Proposition 5.7], $F_{\delta}(Z)$ $\Gamma$-converges to
$\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{2,1}$. ∎

Unlike the $\lambda$ parameter, there is no trade-off for the parameter
$\delta$ in the Huber function per se. The smaller it is, the closer the
recovered solution will be to the true solution. However the trade-off will be
introduced by the fact that we are using a smooth method (conjugate gradient)
to solve a non-smooth problem. If the parameter is too small we may end up
with a non-smooth problem which we are trying to avoid.

We remark that it is possible to obtain a differentiable objective from (12)
with other strategies, such as using dummy variables, see [12] for a
development of this approach for $\ell_{1}/\ell_{2}$ minimization. The
conjugate gradient (CG) algorithm is applied to solve the Huber regularization
problem (13). Since this is a non-convex problem, the iteration may be trapped
at some local minimum. For example, when $X$ is of full column rank ($r=s$),
we know from Section 4.1 that the global minimum of the objective function is
equal to $s$; however, the CG algorithm fails to find it in some cases. To
deal with this difficulty, we run the algorithm with different starting points
that are randomly selected in $\mathbb{S}(N,r)$ using random normally
distributed matrices. In our numerical tests, the global minima can be reached
this way after only a few choices of initial values. Certainly, for large-
scale problems, advanced methods like simulated annealing [20, 16] can be
employed to reduce the number of starting values required.

Figure 1: Huber function for $\delta=0.3$ and the $y=x$ function.  Figure 2:
Derivative of Huber function.

##  6 Explicit formulae for the Euclidean gradients of cost functions

The Euclidean gradients are the fundamental building block for the
optimization of cost functions on manifolds. Such derivatives are often taxing
to calculate, posing a major hurdle for implementing manifold optimization
methods for many problems, especially where large matrices or vectors are
involved. Interestingly, for our objective functions (12) and (13), it is
possible to derive the exact, explicit formulae for their gradients.
Expectedly, these formulae save the users a considerable amount of time on
computing either the numerical or symbolic differentiation.

Throughout this section, $\nabla\equiv\nabla_{Z}$, the derivative with respect
to $Z$. Recall $I_{N}$ is the identity matrix in $\mathbb{R}^{N\times N}$. It
can be easily checked that

$\nabla\|AZ-V\|_{F}^{2}=2A^{T}(AZ-V),$

therefore we will focus on computing the gradients
$\nabla\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{2,1}$ and $\nabla F_{\delta}(Z)$.

Let us introduce some additional notation. For a matrix $Z\in\mathbb{S}(N,r)$,
denote by $D(Z)$ and $D_{\delta}(Z)$ the diagonal matrices in
$\mathbb{R}^{N\times N}$ with

$D(Z)_{i,i}=\\!\begin{cases}\frac{1}{\|Z[i,:]\|_{2}},&\|Z[i,:]\|_{2}\neq 0,\\\
0,&\text{otherwise,}\end{cases}$ $\ \ \text{ and }\ \
D_{\delta}(Z)_{i,i}=\\!\begin{cases}\frac{H_{\delta}^{\prime}(\|Z[i,:]\|_{2})}{\|Z[i,:]\|_{2}},&\|Z[i,:]\|_{2}\neq
0,\\\ 0,&\text{otherwise,}\end{cases}$

where $H^{\prime}_{\delta}(x)$ is the derivative of the Huber function

$H^{\prime}_{\delta}(x)=\begin{cases}1,&x\geq\delta,\\\
\frac{x}{\delta},&0\leq x<\delta.\end{cases}$

Also denote by $\odot$ the Schur-Hadamard (a.k.a., elementwise) product
between two matrices of the same size. The following proposition provides an
explicit formula for the cost function (12).

######  Proposition 3.

Let $Z=U\Lambda V$ be the compact singular value decomposition of
$Z\in\mathbb{S}(N,r)$. Then

$\nabla\|Z(Z^{T}Z)^{-\frac{1}{2}}\|_{2,1}=(I_{N}-UU^{T})D(UV)U\Lambda^{-1}V.$

###### Proof.

For a scalar function $f:\mathbb{R}\to\mathbb{R}$, by an abuse of notation,
the matrix function induced by $f$ is defined on $\mathbb{S}(N,r)$ as
$f(Z)=U\Lambda_{f}V$, where $\Lambda_{f}$ is the $r\times r$ diagonal matrix
such that

$(\Lambda_{f})_{ii}=f(\Lambda_{ii}).$

Note that $f(Z)=Z(Z^{T}Z)^{-\frac{1}{2}}$ for $f\equiv 1$. Let
$g(Z)=\|Z\|_{2,1}$, our goal is to compute $\nabla(g\circ f)$.

Let $L_{g}(Z,E)$ denote the Fréchet derivative of $g$ at $Z$ in the direction
$E\in\mathbb{R}^{N\times r}$ (and similarly for $f$, $g\circ f$, etc.). We
will compute $\nabla(g\circ f)$ via

$\textnormal{tr\,}(\nabla(g\circ f)\cdot E^{T})=L_{g\circ f}(Z,E).$ (14)

First, by the chain rule

$L_{g\circ f}(Z,E)=L_{g}(f(Z),L_{f}(Z,E)).$

Notice that $\nabla g=D(Z)\cdot Z$, hence

$L_{g}(f(Z),L_{f}(Z,E))=\textnormal{tr\,}(D(f(Z))\cdot f(Z)\cdot
L_{f}(Z,E)^{T})=\textnormal{tr\,}(D(UV)\cdot f(Z)\cdot L_{f}(Z,E)^{T}).$

Let $F(Z)\in\mathbb{R}^{r\times r}$ and $G(Z)\in\mathbb{R}^{(N-r)\times r}$ be
defined by

$F(Z)_{i,j}=\begin{cases}\frac{1}{\lambda_{i}+\lambda_{j}},&i\neq j,\\\
0,&i=j,\end{cases}\ \ \ \ \text{ and }\ \ \ \
G(Z)_{i,j}=\frac{1}{\lambda_{j}},$

where $\lambda_{j}$ is the $j$-th singular value of $Z$. Let
$Z=\left[\begin{matrix}U\;W\end{matrix}\right]\left[\begin{matrix}\Lambda\\\
0\end{matrix}\right]V$ be the full SVD of $Z$. From [28, Corollary 3.12],

$\begin{split}L_{f}(Z,E)&=\left[\begin{matrix}U\;W\end{matrix}\right]\left(\left[\begin{matrix}F\\\
G\end{matrix}\right]\odot\left[\begin{matrix}U^{T}EV^{T}\\\
W^{T}EV^{T}\end{matrix}\right]\ +\left[\begin{matrix}-F\\\
0\end{matrix}\right]\odot\left[\begin{matrix}VE^{T}U\\\
W^{T}EV^{T}\end{matrix}\right]\ \right)V\\\
&=U(F\odot(U^{T}EV^{T}-VE^{T}U))V+W(G\odot(W^{T}EV^{T}))V\\\
&=U(F\odot(U^{T}EV^{T}-VE^{T}U))V+WW^{T}EV^{T}\Lambda^{-1}V\\\
&=U(F\odot(U^{T}EV^{T}-VE^{T}U))V+(I_{N}-UU^{T})EV^{T}\Lambda^{-1}V,\end{split}$

using the fact that $G\odot A=A\Lambda^{-1}$, $\forall
A\in\mathbb{R}^{(N-r)\times r}$ and $WW^{T}=I_{N}-UU^{T}$. Thus we have that

$\begin{split}L_{g\circ f}(Z,E)&=\textnormal{tr\,}(D(UV)UV\cdot
V^{T}[F\odot(U^{T}EV^{T}-VE^{T}U)]^{T}U^{T})\\\
&\qquad\qquad\qquad+\textnormal{tr\,}(D(UV)UVV^{T}\Lambda^{-1}VE^{T}(I_{N}-UU^{T}))\\\
&=\textnormal{tr\,}(U^{T}D(UV)U[F\odot(U^{T}EV^{T}-VE^{T}U)]^{T})\\\
&\qquad\qquad\qquad+\textnormal{tr\,}(D(UV)U\Lambda^{-1}VE^{T}(I_{N}-UU^{T}))\\\
&=\textnormal{tr\,}((U^{T}EV^{T}-VE^{T}U)\cdot[F\odot(U^{T}D(UV)U)]^{T})\\\
&\qquad\qquad\qquad+\textnormal{tr\,}((I_{N}-UU^{T})D(UV)U\Lambda^{-1}VE^{T}).\end{split}$
(15)

Here we used the fact that $\textnormal{tr\,}(A\cdot(B\odot
C)^{T})=\textnormal{tr\,}(C\cdot(B\odot A)^{T}).$ Continuing our computations

$\begin{split}&\textnormal{tr\,}((U^{T}EV^{T}-VE^{T}U)\cdot[F\odot(U^{T}D(UV)U)]^{T})\\\
=\
&\textnormal{tr\,}(U^{T}EV^{T}\cdot[F\odot(U^{T}D(UV)U)]^{T})-\textnormal{tr\,}(VE^{T}U\cdot[F\odot(U^{T}D(UV)U)]^{T})\\\
=\ &\textnormal{tr\,}([F\odot(U^{T}D(UV)U)]\cdot
VE^{T}U)-\textnormal{tr\,}(U[F\odot(U^{T}D(UV)U)]^{T}VE^{T})\\\ =\
&\textnormal{tr\,}(U([F\odot(U^{T}D(UV)U)]-[F\odot(U^{T}D(UV)U)]^{T})VE^{T})=0,\end{split}$

which, in combining with (14)–(15), proves our assertion. ∎

Explicit formula for the cost function (13) can be easily derived from the
above proof.

######  Proposition 4.

Let $Z=U\Lambda V$ be the compact singular value decomposition of
$Z\in\mathbb{S}(N,r)$. Then

$\nabla F_{\delta}(Z)=(I_{N}-UU^{T})D_{\delta}(UV)U\Lambda^{-1}V.$

###### Proof.

Notice that the functional

$g_{\delta}(Z)=\sum_{n=1}^{N}\left[H_{\delta}(\|Z[n]\|_{2})+\frac{\delta}{2}\right]$

has the gradient $\nabla g_{\delta}=D_{\delta}(Z)\cdot Z$ and the rest of the
computations follows the same arguments as in Proposition 3. ∎

##  7 Numerical experiments with Gaussian data and measurements

We use the Pymanopt package [31], which is the Python version of the original
Manopt package [4, 34] for MATLAB, in these numerical tests. The conjugate
gradient method with backtracking Armijo line-search provided in Pymanopt is
used to solve (13). For automatic differentiation, we employ the Autograd
package. In each iteration, an $r\times r$ SVD is computed. For sparse
recovery, since $r\leq s$, this step is not numerically expensive. The initial
value is chosen randomly.

In this test we compare the number of measurements required by $\ell_{2,1}$
minimization and our method in sparse matrix recovery. We generate a random
Gaussian matrix $A$ of size $80\times 300$, and form measurement matrices
$A_{k}$ by taking the first $k$ rows of $A$. We also generate a random
Gaussian matrix $X$ of sparsity $s=30$ of size $300\times 70$ for the solution
and test the recovery of $X$ given the output $Y_{k}=A_{k}X$, for each
$k\in\\{38,40,42,\dots,80\\}$. (Hence, $k$ stands for the number of
measurements). The test is rerun 22 times for different realizations of $A$
and $X$. We perform the $\ell_{2,1}$ minimization using the spgl1 package [3]
(which is specialized to solve $\ell_{1}$ and related regularized least square
problems), and our method using Pymanopt, with $\delta=0.001$ and $\lambda=9$,
until relative reduction of the gradient norm is $10^{-8}$ or 1000 iterations
have been run. For each value of $k$, we plot the scattered errors and median
of recovery errors in log scale over $22$ trials. The results of the
experiment are presented in Figure 3. As can be seen, our manifold
optimization algorithm starts converging with much fewer measurements compared
to the Euclidean minimization with the $\ell_{2,1}$ norm.

Figure 3: Scattered errors and median of errors in log scale vs. number of
measurements for $\ell_{2,1}$ (circular labels) and manifold optimization (13)
(triangular labels).

##  8 Conclusion and future work

We study the joint sparse recovery problem. Our result here allows to reduce
the problem to a potentially much smaller dimensional problem where the
unknown matrix has full column rank. We also offer a new functional
minimization method on the manifold of full-rank matrices for the recovery of
the unknown matrix. As pointed out in Section 1.2, we cannot expect to have a
rank sensitive convex optimization method on the Euclidean space, so changing
the geometry and looking into manifold optimization is probably a good
direction to go. The functional we suggest does encourage sparsity and in our
numerical experiment outperforms the Euclidean $\ell_{2,1}$ minimization. We
pose the following questions as potential directions for future investigation.

  1. 1.

If $A\in\mathbb{R}^{M\times N}$ satisfies the NSP condition, then is it true
that any $s$-sparse $x\in\mathbb{R}^{N}\setminus\\{0\\}$ can be uniquely
recovered by solving

$\operatorname*{arg\,min}_{z\in\mathbb{R}^{N}}\|z\|_{1}/\|z\|_{2}\text{ s.t.
}Az=Ax?$

More generally, does NSP imply unique recovery of every $s$-row sparse
$W\in\mathbb{R}^{N\times r}$ from (11), for any $1\leq r\leq s$?

  2. 2.

Is (11) rank aware? More specifically, are there measurement matrices $A$ that
recover every $W\in\Sigma_{s,r}$ by solving (11) for large $r$’s but not for
small $r$’s?

  3. 3.

Does there exist a sparsity-promoting geodesically convex manifold
optimization method for joint sparse recovery?

  4. 4.

Investigation of robustness and stability of our method in presence of noise.

##  9 Acknowledgments

We want to thank Richard Barnard (Western Washington University) for
suggesting the use of Huber regularization and Carola-Bibiane Schönlieb
(University of Cambridge) for pointing out the $\Gamma$-convergence argument
for Huber regularization to us. We also thank the anonymous referee for
valuable suggestions.

This work is supported by the National Science Foundation, Division of
Mathematical Sciences, Computational Mathematics program under contract number
DMS1620280; the U.S. Department of Energy, Office of Science, Office of
Advanced Scientific Computing Research, Applied Mathematics program under
contract number ERKJ314; Scientific Discovery through Advanced Computing
(SciDAC) program through the FASTMath Institute under Contract No. DE-
AC02-05CH11231; and, the Laboratory Directed Research and Development program
at the Oak Ridge National Laboratory, which is operated by UT-Battelle, LLC.,
for the U.S. Department of Energy under Contract DE-AC05-00OR22725.

## References

  * Absil et al. [2009] P.A. Absil, R. Mahony, R. Sepulchre, Optimization algorithms on matrix manifolds, Princeton University Press, 2009\.
  * Van den Berg and Friedlander [2011] E. Van den Berg, M. Friedlander, Sparse optimization with least-squares constraints, SIAM Journal on Optimization 21 (2011) 1201–1229.
  * van den Berg and Friedlander [2008] E. van den Berg, M.P. Friedlander, Probing the pareto frontier for basis pursuit solutions, SIAM Journal on Scientific Computing 31 (2008) 890–912.
  * Boumal et al. [2014] N. Boumal, B. Mishra, P.A. Absil, R. Sepulchre, Manopt, a Matlab toolbox for optimization on manifolds, Journal of Machine Learning Research 15 (2014) 1455–1459.
  * Candès et al. [2006] E. Candès, J. Romberg, T. Tao, Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information, IEEE Transactions on information theory 52 (2006) 489–509.
  * Chen and Huo [2006] J. Chen, X. Huo, Theoretical results on sparse representations of multiple-measurement vectors, IEEE Trans. Signal Process. 54 (2006) 4634–4643.
  * Chen et al. [2011] Y. Chen, N.M. Nasrabadi, T. Tran, Hyperspectral image classification using dictionary-based sparse representation, IEEE Transactions on Geoscience and Remote Sensing 49 (2011) 3973–3985.
  * Dal Maso [2012] G. Dal Maso, An introduction to $\Gamma$-convergence, volume 8, Springer Science & Business Media, 2012\.
  * Dalalyan et al. [2017] A.S. Dalalyan, M. Hebiri, J. Lederer, et al., On the prediction performance of the lasso, Bernoulli 23 (2017) 552–581.
  * Davies and Eldar [2012] M.E. Davies, Y.C. Eldar, Rank awareness in joint sparse recovery, IEEE Trans. Inform. Theory 58 (2012) 1135–1146.
  * Donoho [2006] D.L. Donoho, Compressed sensing, IEEE Trans. Inform. Theory 52 (2006) 1289–1306.
  * Esser et al. [2013] E. Esser, Y. Lou, J. Xin, A method for finding structured sparse solutions to non-negative least squares problems with applications, SIAM J. Imaging Sciences 6 (2013) 2010–2046.
  * Fang et al. [2014] L. Fang, S. Li, X. Kang, J. Benediktsson, Spectral–spatial hyperspectral image classification via multiscale adaptive sparse representation, IEEE Transactions on Geoscience and Remote Sensing 52 (2014) 7738–7749.
  * Feng [1998] P. Feng, Universal Minimum-Rate Sampling and Spectrum-Blind Reconstruction for Multiband Signals, Ph.D. thesis, University of Illinois, Urbana-Champaign, 1998.
  * Feng and Bresler [1996] P. Feng, Y. Bresler, Spectrum-blind minimum-rate sampling and reconstruction of multiband signals, Proc. ICASSP 3 (1996) 1688–1691.
  * Harjek [1988] B. Harjek, Cooling schedules for optimal annealing, Mathematics of Operation Research 13 (1988) 311–329.
  * Hoyer [2004] P. Hoyer, Non-negative matrix factorization with sparseness constraints, Journal of Machine Learning Research 12 (2004) 1457–1469.
  * Huang et al. [2018] W. Huang, P.A. Absil, K.A. Gallivan, P. Hand, Roptlib: an object-oriented c++ library for optimization on riemannian manifolds, ACM Transactions on Mathematical Software (TOMS) 44 (2018) 43.
  * Iordache et al. [2014] M.D. Iordache, J.M. Bioucas-Dias, A. Plaza, Collaborative sparse regression for hyperspectral unmixing, IEEE Transactions on Geoscience and Remote Sensing 52 (2014) 341–354.
  * Kirkpatrick et al. [1983] S. Kirkpatrick, C.D. Gelatt, M.P. Vecchi, Optimization by simulated annealing, Science 220 (1983) 671–680.
  * Krishnan et al. [2011] D. Krishnan, T. Tay, R. Fergus, Blind deconvolution using a normalized sparsity measure, Proc. of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (2011).
  * Lai and Liu [2011] M.J. Lai, Y. Liu, The null space property for sparse recovery from multiple measurement vectors, Applied and Computational Harmonic Analysis 30 (2011) 402 – 406.
  * Lederer and Müller [2015] J. Lederer, C. Müller, Don’t fall for tuning parameters: tuning-free variable selection in high dimensions with the trex, Twenty-Ninth AAAI Conference on Artificial Intelligence (2015).
  * Lee et al. [2012] K. Lee, Y. Bresler, M. Junge, Subspace methods for joint sparse recovery, IEEE Transactions on Information Theory 58 (2012) 3613–3641.
  * Malioutov et al. [2005] D. Malioutov, M. Cetin, A. Willsky, A sparse signal reconstruction perspective for source localization with sensor arrays, IEEE transactions on signal processing 53 (2005) 3010–3022.
  * Mishali and Eldar [2008] M. Mishali, Y.C. Eldar, Reduce and boost: Recovering arbitrary sets of jointly sparse vectors, IEEE Transactions on Signal Processing 56 (2008) 4692–4702.
  * Natarajan [1995] B.K. Natarajan, Sparse approximate solutions to linear systems, SIAM J. Comput. 24 (1995) 227–234.
  * Noferini [2017] V. Noferini, A formula for the fréchet derivative of a generalized matrix function, SIAM Journal on Matrix Analysis and Applications 38 (2017) 434–457.
  * Rahimi et al. [2018] Y. Rahimi, C. Wang, H. Dong, Y. Lou, A scale invariant approach for sparse signal recovery, arXiv preprint arXiv:1812.08852 (2018).
  * Schmidt [1979] R.O. Schmidt, Multiple emitter location and signal parameter estimation, Proc. RADC Spectral Estimation Workshop (1979) 243–258.
  * Townsend et al. [2016] J. Townsend, N. Koep, S. Weichwald, Pymanopt: A python toolbox for optimization on manifolds using automatic differentiation, Journal of Machine Learning Research 17 (2016) 1–5.
  * Tran and Webster [2017] H. Tran, C. Webster, Unified sufficient conditions for uniform recovery of sparse signals via nonconvex minimizations, arXiv:1710.07348 (2017).
  * Udriste [1994] C. Udriste, Convex Functions and Optimization Methods on Riemannian Manifolds, volume 297 of Mathematics and Its Applications, Springer Netherlands, 1 edition, 1994\.
  * Vandereycken [2013] B. Vandereycken, Low-rank matrix completion by Riemannian optimization, SIAM Journal on Optimization 23 (2013) 1214–1236.
  * Wang et al. [2017] Z. Wang, J. Liu, J. Xue, Joint sparse model-based discriminative k-svd for hyperspectral image classification, Signal Processing 133 (2017) 144–155.
  * Yin et al. [2014] P. Yin, E. Esser, J. Xin, Ratio and difference of $\ell_{1}$ and $\ell_{2}$ norms and sparse representation with coherent dictionaries, Commun. Inform. Systems 14 (2014) 87–109.
  * Zhu et al. [2014] F. Zhu, Y. Wang, S. Xiang, B. Fan, C. Pan, Structured sparse method for hyperspectral unmixing, ISPRS Journal of Photogrammetry and Remote Sensing 88 (2014) 101 – 118.