Orthogonally weighted $\ell_{2,1}$ regularization for rank-aware joint sparse recovery: algorithm and analysis

The joint sparse recovery or multiple measurement vector (MMV) problem aims to find a solution to the matrix equation $AX=Y$, where $A\in\mathbb{R}^{M\times N}$ and $Y\in\mathbb{R}^{M\times K}$ are given, and $X\in\mathbb{R}^{N\times K}$ is an unknown $s$-row sparse matrix, meaning that at most $s$ rows of $X$ are non-zero.

The task of finding the solution of $AX=Y$ with the fewest non-zero rows in the noiseless case can be formulated as an optimization problem using mixed $\ell_{2,0}$ norms:

where $\ell_{2,0}(Z)$ is the number of non-zero rows in $Z$. However, the main drawback of this method is that the optimization problem is computationally infeasible, as it results in an exhaustive search that is an NP-hard problem. To overcome this limitation, the $\ell_{2,1}$ norm regularization method is commonly used to solve the joint sparse recovery problem [51]

where $\ell_{2,1}(Z)=\lVert Z\rVert_{2,1}=\sum_{n=1}^{N}\sqrt{z_{n}z_{n}^{T}}$, and $z_{n}$ is the $n$-th row of $Z$. This method generalizes the $\ell_{1}$ sparsifying penalty to the multiple vector setting and is computationally feasible. However, this method has the drawback of being rank-blind, meaning it does not take into account the rank of the $s$-row sparse matrix [16]. In real-world scenarios, data typically has a maximum essential rank, equal to its essential sparsity.

In [42], it was argued that a rank-aware regularizer cannot be continuous and, thus, convex in all of $\mathbb{R}^{N\times K}$. A new method for recovery of jointly sparse vectors in the form of a (non-convex) optimization problem was proposed

where the regularizer, here we call the orthogonally weighted $\ell_{2,1}$ or $\mathit{ow}\ell_{2,1}$, is defined as

with $(\cdot)^{\dagger}$ denoting the Moore-Penrose pseudo-inverse and $(\cdot)^{\dagger/2}=((\cdot)^{\dagger})^{1/2}$ its square root for a given symmetric matrix. If the minimization is restricted to matrices of full rank, we obtain the functional $\mathit{ow}\ell_{2,1}(Z)=\ell_{2,1}(Z(Z^{T}Z)^{-1/2})$. The papers [36, 42] also discussed a strategy to reduce the joint sparse problem (1) to a simpler setting where the measurement matrix has full column-rank. This transformation can be achieved through techniques such as Singular Value Decomposition (SVD) or other factorization methods. In particular, if $Y=Y^{\prime}Q$ where $Y^{\prime}\in\mathbb{R}^{M\times r}$ with ${\rm rank\,}(Y^{\prime})=r$ and $Q\in\mathbb{R}^{r\times K}$ with $QQ^{T}=I\in\mathbb{R}^{r\times r}$, then for any solution $Z^{\prime}$ of the optimization problem

the matrix $\bar{Z}=Z^{\prime}Q$ is a solution of original problem (1). This strategy is independent of the chosen regularization functional and can be used for the noisy problem as well to reduce it to an approximate full column rank problem, and the resulting reconstruction error can be controlled; see [42] for more details.

However, [42] neither provided an affirmative conclusion on the rank-awareness, nor provided a general algorithm for implementation of $\mathit{ow}\ell_{2,1}$ regularization. The algorithm proposed in [42] relies on being able to determine the rank of the solution directly from the data and does not directly address the more practical case of noisy data discussed below. In this paper, we aim to address these two important open questions by providing: i) a proof of the rank-awareness of the $\mathit{ow}\ell_{2,1}$ regularization, making it the first regularization-based method possessing this property; and ii) a computationally efficient algorithm for solving the optimization problem, supported by rigorous convergence analysis.

The naming of the functional is derived from the following. Note that, if $Z=U\Sigma V^{T}$ is the compact singular value decomposition of $Z$, i.e., with $U\in\mathbb{R}^{N\times r}$, $V\in\mathbb{R}^{K\times r}$ orthogonal, $r={\rm rank\,}(Z)$, and $\Sigma\in\mathbb{R}^{r\times r}$ diagonal and s.p.d., then $Z(Z^{T}Z)^{\dagger/2}=UV^{T}$, thus it renormalizes $Z$ to its orthogonal component before applying the $\ell_{2,1}$ norm. Also

in other words, $\mathit{ow}\ell_{2,1}(Z)$ is the $\ell_{2,1}$ norm of an orthonormal basis of the column space of $Z$ and its value is independent of the choice of the basis. This and many other properties of the $\mathit{ow}\ell_{2,1}(Z)$ functional, like its locally Lipschitz continuity, are discussed in detail in Appendix A. Consequently, when $Z$ is of $s$-row sparsity and of rank $s$, it can be seen that $\mathit{ow}\ell_{2,1}(Z)=\ell_{2,0}(Z)$, in other words it interpolates the values of $\ell_{2,0}$ on such matrices. However, unlike the latter, it is feasible to develop an efficient algorithm for solving the optimization problem with the $\mathit{ow}\ell_{2,1}$ regularizer, as we shall see in this effort.

In the presence of measurement noise, we consider the regularized formulation of (3)

The term $1/(2\alpha)\|AZ-Y\|_{\textnormal{Fro}}^{2}$ acts as a penalty for deviating from the measurements, where $\|\cdot\|_{\textnormal{Fro}}$ denotes the Frobenius norm of a matrix. The quantity $1/(2\alpha)$ controls the strength of the regularization or can be interpreted as a Lagrange multiplier for the constrained formulation

To simplify notation, we will primarily consider the regularized formulation of the problem as given in (4).

$\ell_{1}/\ell_{2}$ regularization. Our method represents a generalization of $\ell_{1}/\ell_{2}$, which uses the ratio of $\ell_{1}$ and $\ell_{2}$ to model the sparsity [28, 21]. Indeed, in the special case of single vector recovery (i.e. for $K=1$), the problems (3) and (4) are equivalent to the $\ell_{1}/\ell_{2}$ regularized problem due to

Theoretical analysis for $\ell_{1}/\ell_{2}$ has been developed in [56] for non-negative signals, [43] with a local optimality condition, and [46], where an analytic solution for the proximal operator is characterized, among others. Our proof of the existence of the global solution to (3) and (4) here is built on a generalization of the spherical section property, which has been studied in [57, 55] for the global optimality condition of $\ell_{1}/\ell_{2}$. The $\ell_{1}/\ell_{2}$ functional has been shown to enhance sparsity of the solutions and outperform $\ell_{1}$ on many test problems such as nonnegative matrix factorization and blind deconvolution [31, 29]. However, the existing uniform recovery guarantees for $\ell_{1}/\ell_{2}$ remains too strict in the general, noisy setting (i.e., more restrictive than well-known $\ell_{1}$ recovery conditions such as the null space property), despite many recent theoretical advances.

Joint sparse recovery. The joint sparse recovery problem, also known as multiple measurement vectors, distributed compressed sensing or simultaneous sparse approximation, has garnered a lot of attention in the literature. Mathematically, this problem arises when we want to reconstruct multiple vectors which have the same sparsity support using the same set of linear combinations for each. While one can certainly recover each vector independently, approximating them all at once should give superior results because we can take advantage of additional information about their shared support. To realize this potential, however, a joint sparse recovery algorithm is desired to be rank aware, that is, its performance improves with increasing rank of the solution matrix.

Existing approaches for the joint sparse recovery can be roughly categorized into greedy methods [14, 48, 25], and algorithms based on mixed norm optimization [14, 47, 23, 19, 50]. Pioneering joint sparse methods are often straightforward extensions of algorithms for single vector recovery, whose recovery conditions are based on standard compressed sensing properties and typically turn out to be the same as those obtained from a single measurement problem. In such cases, the theoretical results do not show any advantage of the joint sparse approach. One of the first methods exploiting rank information is a discrete version of the MUSIC (MUltiple SIgnal Classification) algorithm [22], where recovery guarantees better than those from single vector problem are proved when the solution matrix has full column-rank. Extensions of this work to non-full-rank cases have been developed in [33] (subspace-augmented MUSIC), [30] (compressive MUSIC) and [54] (semi-supervised MUSIC). Also inspired by the MUSIC technique, [16] proposes the rank-aware orthogonal recursive matching pursuit (RA-ORMP) based on a greedy selection strategy, and further demonstrates the rank awareness and rank blindness for a large class of joint sparse recovery algorithms. That work proves that some of the most popular algorithms, including orthogonal matching pursuit and mixed $\ell_{q,1}$ norm for any $q\geq 1$, are effectively rank blind. More recent rank aware methods include Simultaneous Normalsize Iterative Hard Thresholding (SNIHT) [6] and Rank Aware Row Thresholding (RART) [7]. Compared to the greedy approach, the integration of rank information into regularization-based algorithms seems an open problem. To date, all existing rank-aware algorithms belong to the greedy class.

Applications. The joint sparse recovery problem finds applications in various fields, including data analytics and imaging [34, 24], uncertainty quantification [17], and parameter identification of dynamical systems [44]. Here, we will illustrate one particular area of our interest in detail, which is the Column Subset Selection Problem (CSSP), also known as the feature selection, or the dictionary selection problem. Suppose we have a set of data columns $Y\in\mathbb{R}^{M\times K}$ and a dictionary of features $A\in\mathbb{R}^{M\times N}$ (here $N$ represents the number of features in the dictionary and $M$ is the size of the features), and we want to extract a subset of the columns of $A$ such that $Y$ can be represented in terms of only $s$ columns of $A$. This problem exactly translates to the joint sparse recovery problem $Y=AX$ where $X$ is the matrix of coefficients and the non-zero rows correspond to the columns that have been selected. In the feature selection problem, the data is also the matrix $A$ (i.e., $Y=A$) and the goal is to select the most representative columns (features) of the matrix $A$. Thus, we aim to solve the problem

This problem has garnered extensive attention in machine learning where it was often relaxed by mixed $\ell_{2,1}$ norm [34, 38, 35, 37]. This has also led to various matrix decomposition techniques such as the CUR decomposition [26, 8]. Demonstrations of our $\mathit{ow}\ell_{2,1}$ regularization on this problem on some bioinformatics datasets will be shown in Section 7.

Proximal algorithms. To our knowledge, this is the first paper that provides an iterative solution method for the problem (4) aside from the initial work [42]. The methods we outline in Section 6 will be based on standard proximal gradient based methods for least squares regression with convex regularizers [3, 9], and for more general non-smooth, convex optimization that encompasses joint sparse regression [11, 23, 12, 45, 18]. In order to incorporate the nonconvex non-smooth aspects of the $\mathit{ow}\ell_{2,1}$-regularizer, we employ a first order convex approximation through a model function in every step of the method. This allows us to write the new iteration with the help of a proximal mapping that is simple to calculate. Convex model functions are also used in related methods for nonsmooth nonconvex optimization [40, 39] or Trust-Region methods [13]. Our method employs a specialized variable metric to achieve simple and efficient iterations, which requires a dedicated analysis and leads to a global convergence result.

As already pointed out before, for the scalar case $K=1$ we obtain the $\ell_{1}/\ell_{2}$ regularizer, and optimization methods have been proposed in [56, 43, 53, 55, 46, 57]. For $K=1$, proximal method in this manuscript is similar to the one proposed in [53, Algorithm 2] for the equality constrained formulation (3) and subsequently analyzed and generalized to include the noisy constrained formulation (5) in [57]. In this paper we focus on the noisy regularized formulation (4), which can also be used to approach the solutions of the equality constrained problem for $\alpha\to 0$ and the solutions of the noisy constrained problem for $\alpha\to\alpha_{\delta}$; see Section 6.3. While it is easy to specialize the method proposed in this manuscript to $K=1$ (see Remark 6.5), the generalization to larger $K$ is not obvious, and requires additional effort.

Finally, since the $\mathit{ow}\ell_{2,1}$ functional is discontinuous across different ranks, we can not rely on only gradient descent unless we know or can guess the optimal rank ahead of time [42]. As a remedy, we propose an optional Lipschitz-continuous relaxation of the problem for a parameter $\gamma>0$, which we show converges to the original problem in the sense of $\Gamma$-convergence [15].

Denote for any $r\leq K$

When the columns of $Z$ are linearly independent ($r=K$), we abuse the notation and write $\mathbb{S}(N,r)=\mathbb{S}(N,K,r)$. It is worth noting that the set $\mathbb{S}(N,K,r)$ is often referred to as the Stiefel manifold and it is a key object in differential geometry and optimization, with numerous applications in machine learning and data analysis.

The following theorem [32, 16] provides a necessary and sufficient condition for the exact recovery of sparse matrices from the $\ell_{2,1}$ minimization problem (2):

The following theorem (see [42] for details) highlights the difference in behavior between the $\ell_{2,1}$ norm and the $\ell_{2,0}$ norm regularizers. The $\ell_{2,1}$ norm regularizer is not rank-aware, while the $\ell_{2,0}$ norm regularizer is rank-aware. The smallest number of linearly dependent columns in a matrix $A$ is referred to as the spark of $A$ and is denoted by $\textnormal{spark\,}(A)$.

In the above theorem, the least restrictive condition occurs when the rank of the row-sparse matrix $X$ is equal to its sparsity, i.e., $r=s$. In fact, in practical applications, the row-sparse matrix $X$ in the joint sparse recovery problem is often of maximum rank, i.e., it is approximately $s$-row and also $s$-row sparse. In such a scenario, the condition specified in Theorem 3.2 simplifies to $s<\textnormal{spark\,}(A)$. Hence, the assumption $s<\textnormal{spark\,}(A)$ is consistently made throughout the manuscript.

For a given matrix $A$, if ${\textnormal{spark\,}(A)}/{2}<s<\textnormal{spark\,}(A)$ and $r<s$, it is not necessarily guaranteed that the solution to the joint sparse recovery problem exists or is unique. In this context, Section 4 introduces a condition named spherical section property to ensure the existence of the global solution.

Earlier, we discussed that when the rank is equal to one, the $\mathit{ow}\ell_{2,1}$ reduces to $\ell_{1}/\ell_{2}$. Below we consider examples in which the $\ell_{1}/\ell_{2}$ fails to recover a single vector but if we add another vector along with the corresponding measurements, then we successfully recover both.

The next example, in a similar spirit, demonstrates that the functional can even fail to have a minimizer for large $\textnormal{spark\,}(A)$, if the rank is smaller than the sparsity.

Here we turn to the question of showing that the global minimizers of the noise free (3) and the noisy unconstrained (4) problems indeed exist under certain conditions on $A$. They do not have to exist in general.

To illustrate this fact, let us consider for the moment the case $K=1$, where the functional $\mathit{ow}\ell_{2,1}$ is identical to the $\ell_{1}/\ell_{2}$ functional (6). For the existence of the global minimum of the cost function with the $\ell_{1}/\ell_{2}$ regularizer, it is often required that $A$ satisfies the $s$-spherical section property. An example can be found in [57, Theorem 3.4]. In this section we aim to extend these results to joint sprase recovery problem with the $\mathit{ow}\ell_{2,1}$ regularizer. The failure of minimizers to exist without appropriate assumptions has been illustrated in the previous section.

This set contains all matrices of rank $k$ where some linear combination of the columns is in the kernel of $A$. Notice that the rank-$1$ singularity set for $K=1$ is the nontrivial part of the kernel: For $K=1$ we have $\Omega_{1}(A)=\ker(A)\setminus\{\,0\,\}$. This justifies the following definition as an extension of the $s$-spherical section property (cf. [52, 57]).