Density Corrected Sparse Recovery when R.I.P. Condition Is Broken

1 Introduction

In high dimensional statistics, an important problem is to recover a sparse vector from a small number of observations whose number is usually much smaller than the nominal dimension of the observation matrix. [Candes et al., 2006] proves that sparse recovery is possible as long as the feature matrix satisfies the Restricted Isometric Property (R.I.P.) condition. Informally, the R.I.P. condition requires any set of columns of the feature matrix with size less than the sparsity to be linearly independent. More likely than not, the R.I.P. condition is too strong to be useful in real world applications. Many R.I.P.-like relaxations have been proposed in the past decades [Van De Geer and Bhlmann, 2009; Foucart, 2012]. These alternatives usually have better constants, or more general forms of the restricted isometric inequalities. See [Van De Geer and Bhlmann, 2009] for a nice survey.

Although there is a vast literature studying relaxed R.I.P.- like conditions, their core requirement is that the feature matrix must be well-conditioned on the sparse subspace. In other words, their largest and smallest singular values must be close to unit. For many real applications, such assumptions are still too strong to be useful. Considering a simple case where we duplicate some columns of the feature matrix, these R.I.P.- like conditions are no longer true because we can find a subset of columns whose smallest singular value is exactly zero.

In this paper we consider a prevalent non-R.I.P. setting in which the features form cluster structures, as can be seen in many machine learning [Lehiste, 1976] and computer vision problems [Lan et al., 2013; Lowe, 2004]. Due to the fact that many features extractors are similar to each others and they reflect the characteristics of the same image, vision features are often correlated and have cluster structures. This correlation is even stronger in those systems that have thousands to millions of features [Lan et al., 2013; Gan et al., 2015a; 2015b].

Under this setting, instead of requiring the observed feature matrix to be R.I.P., we only have the same requirement on the cluster center matrix, which is much easier to satisfy. A trivial solution to perform sparse recovery under this setting is by first clustering the features and then applying the conventional sparse recovery methods on the clustered feature matrix. However, due to the random perturbation, the estimated cluster center is the biased version of the real cluster centers. Unlike instance clustering where we can improve the estimation accuracy by having more samples, the bias in feature clustering process cannot be asymptotically removed even if infinite training instances were given. This inconsistency happens in all conventional methods, including convex/non-convex sparse regularizers , greedy methods and their stochastic or adaptive variants [Agarwal et al., 2012; Shalev-Shwartz and Tewari, 2011; Ghadimi and Lan, 2013; Ji Liu, 2013; Jin et al., 2013; Zhaoran Wang, 2013; Lin et al., 2014b; Lin and Xiao, 2014; Yang et al., 2014; Lin et al., 2014a].

In this paper, under cluster assumption, we develop a consistent sparse recovery method, called density correction. The key idea of density correction is that although we cannot eliminate the bias in feature clustering, we are able to estimate the variance of this bias via the cluster density. We can correct the bias in the gradient of the loss function with cluster density in the sparse recovery process. When combining the density correction with hard iterative thresholding, we obtain a consistent sparse estimator.

The remaining paper is organized as following. In section 2 we briefly review several closely related works. We introduce our notations and backgrounds of our method in section 3.

Section 4 gives the details of our algorithm. In section 5 we prove the geometrical convergence rate and the consistency of our method. We conclude our paper with discussions in section 6.

2 RelatedWork

In this section, we briefly review several closely related works. Our review is by no means to be comprehensive but to capture the backgrounds of our approaches.

Sparse recovery methods can be roughly grouped into three categories: convex regularizer, non-convex regularizer and greedy methods. The convex regularizer mainly bases on ℓ1-norm LASSO problems [Shah, 2012; Xu et al., 2010; Tibshirani, 1996]. The non-convex regularizer is proved to be more accurate than convex methods in several cases [Zhaoran Wang, 2013; Xiang et al., 2013; Loh and Wainwright, 2013; Gong et al., 2013; Zhang, 2012]. Recently, greedy methods are rediscovered in the sparse community due to its simplicity and efficiency [Xiaotong Yuan, 2014; Ji Liu, 2013; Blumensath, 2012; Foucart, 2011]. However, all these studies are based on R.I.P.-like assumptions therefore fail to cover the non-R.I.P. problems discussed in this paper.

The gradient correction technique dates back to the corrected least square regression [Markovsky and Van Huffel, 2007]. The same technique is considered in sparse regression under feature corruption settings [Chen et al., 2013; Loh and Wainwright, 2012]. These methods require that the feature matrix to be R.I.P. when noise is small. In our setting, the observed feature matrix becomes non-R.I.P. when the perturbation is small. Also, traditional methods require the noise level to be known while our method not.

3 Notations and Backgrounds

The sparse recovery problem aims to recover a sparse vector ¯w∗ ∈ Rd from feature matrix ¯X ∈ Rn×d and label y ∈ Rn . The sparsity of ¯w∗ is denoted by s and its support set is denoted as S∗ . We assume that y is generated by a linear model

where ξ is additive subgaussian noise with noise level

We denote ¯Xs the submatrix of ¯X whose columns are indexed by set S. we denote [w]F the truncation of w whose elements are zero outside F. For simplicity, let ∇Fℓ(w) ≜ [∇ℓ(w)]F be the truncated gradient. ∇sℓ(w) is the truncated gradient with largest s elements in absolute value. In this paper we assume that both S∗ and ¯X are unknown. We can only observe feature matrix X whose columns are randomly sampled from a clustering model with cluster center ¯X. That is,

where P is called duplication matrix. The random noise ϵ is subgaussian satisfying

D(·) is the diagonal function. In the i-th cluster, there are ki features. We denote σ =maxi σi .

The feature cluster model defined in Eq. (2) coincides with many real world applications. For example, in computer vision, we usually need to encode low level features by clustering, which is called "feature encoding" in literature. Some features are built on other features. For example, semantic features are built from low-level features hence are correlated to low-level features. In those scenarios, for set S with size s, the submatrix XS is usually ill-conditioned.

The conventional sparse recovery methods assume the feature matrix X to satisfy the R.I.P. [Candes and Plan, 2011] or its variants [Trzasko and Manduca, 2009]. Roughly speaking it requires XS to be well-conditioned for |S| ≤ s. Then we can estimate ¯w∗ via the following optimization problem:

To recover ¯w∗ ,

However in Eq. (2), when ϵ → 0, the R.I.P.-like condition clearly does not hold hence conventional methods P1 are no longer applicable. To address this issue, a naive approach is to cluster the feature matrix and then apply P1 on the clustered features. Denote ˆX the clustered feature matrix. The naive approach can be formulated as

After clustering, ˆX is expected to be R.I.P. if ¯X is R.I.P. . A main drawback of P2 is that w(2) ∗ is not a consistent estimator. The inconsistency of w(2) ∗ comes from the uncertainty of ˆX. Recall that ˆX is an estimation of ¯X by matrix column clustering, from concentration inequality, it is easy to see that with high probability,

Therefore, ˆX is a biased estimation of ¯X. This bias cannot be eliminated by more training instances since the right side of Eq. (3) does not depend on n. As a consequence, P2 which is based on ˆX is biased. Even if we already know the support set S∗ , w(2) ∗ is still biased due to the uncertainty of ˆX. This is known as the least square inconsistency when we are given the support set S∗ . In short, we have the following proposition. We omit the proof since it is trivial once we realize the uncertainty in Eq. (3).

Proposition 1. No matter how many training instances were given, P1 and P2 cannot estimate ¯w∗ consistently even if we know the support set S∗ and the duplication matrix P .

Proposition 1 shows that the conventional sparse recovery methods are far from optimal because they are all inconsistent even in the ideal case. Note that their suboptimality does not inherit from the way they solve the problem P1 or P2. We can use any algorithm to solve the above two problems. For example, directly optimize the ℓ0 norm constraint by greedy methods [Tropp and Gilbert, 2007; Blumensath and Davies, 2008; Xiaotong Yuan, 2014], or use convex [Tibshirani, 1996] or non-convex relaxation methods [Zhang, 2012; Zhaoran Wang, 2013; Lin et al., 2014a]. All these methods cannot eliminate the uncertainty in ˆX. By proposition 1, no consistent estimation is possible by these methods.

A direct corollary of proposition 1 is that it is impossible to consistently estimate ¯X by any other method. Based on Eq. (3), even when we know the duplication matrix P, the best possible accuracy in estimating ¯Xi is O(σi/√ki). This upper bound is based on the concentration of mean value of ki random variables. Because ki is finite, it is impossible to improve this estimation unless we find a better way to estimate the mean of random variables.

In next section, we will show that by using the proposed density correction, we can design a consistent sparse estimator under model Eq. (2) even if we cannot eliminate the uncertainty in ˆX.

4 Density Corrected Sparse Recovery

From Eq. (3), it is impossible to eliminate the uncertainty in ˆX. To motivate an alternative approach, we examine the gradient of loss function in P2,

Denote ˆX = ¯X + θ, we have

Therefore, the bias in the gradient comes from E1 and E2. Since the perturbation ϵ is zero mean subgaussian random variable, θ is also zero mean subgaussian. When n is large enough, E2/n will converge to ¯XTy/n thus is consistent. However, E1 will not converge to zero. This inspires us to eliminate E1 without eliminating θ itself. In other words, although it is impossible to eliminate the bias term θ, we can still get a consistent estimator if we are able to estimate the variance of this bias term. Based on this intuition, we propose our Density Corrected Sparse Recovery (DCSR) in Algorithm 1.

Algorithm 1 is based on Hard Iterative Thresholding (HIT) [Blumensath and Davies, 2009]. It first clusters features into r groups. r could be potentially much larger than the sparsity s. Although we can use any clustering algorithm to cluster features, we prove that the pairwise tree clustering is sufficient to recover the cluster center ¯X with the optimal bias O(σi/√ki) up to constants. For the i-th cluster, its cluster density is stored in νi. Actually, νi is an unbiased estimation of the variance of perturbation term ξ in model Eq. (2). Then the algorithm does a gradient descent with step size η. The

Algorithm 1 Density Corrected Sparse Recovery (DCSR)

• 1: Input: X is the feature matrix . r is the number of columns of ¯X. s is the sparsity of ¯w∗ . T is the total iteration number. γ is the step size in gradient descent.
• 2: Use pairwise tree clustering to cluster columns of X into r groups. The distance between i-th and j-th column is defined by

Denote the index set of i-th cluster as Ci , i = 1, 2, · · · r.

• 3: Denote ˆXi the cluster center and νi the density of Ci. For |Ci| = ki, ˆXi and νi is computed by

• 6: for t = 1 to T do
• 7: Gradient descent with density correction:

• 8: Hard iterative sparse thresholding: let St be the set of indexes of elements in ˆwt of s largest amplitude ,

• 9: end for
• 10: Output: Cluster index Ci and sparse vector wT .

range of η is given in the next section. The key difference is that the gradient is corrected by the cluster density νi in Eq. (7). Eq. (7) can be reformulated as a gradient descent step for the following objective function

This is similar to Elastic-Net [Zou and Hastie, 2005]. In Elastic-Net, an ℓ2-norm regularizer is added in the sparse regression objective function to ensure a faster convergence rate. While in density correction, we subtract an ℓ2-norm regularizer adaptively to remove the bias in the gradient.

After gradient descent, the algorithm truncates ˆwt greedily. Note that the density correction is a principled method. We can apply density correction in various sparse regression methods such as LASSO and adaptive non-convex regularizer. We choose HIT because it is recently proven to be efficient in sparse recovery [Xiaotong Yuan, 2014]. In each iteration, the intermediate solution wt is always s sparse.

5 Theoretical Analysis

In this section, we give theoretical guarantees for Algorithm 1. Our analysis relies on the following matrix concentration facts [Tropp, 2012].

Lemma 1. Let ϵ ∈ Rn×d be a subgaussian random matrix with covariance

The eigenvalue of 1n AAT is bounded by [µ, β]. Then with probability at least 1- δ, denote ˆ A = A + ϵ,

where c is at least 1/4,

We assume the feature cluster center matrix ¯X to be τdistinguishable:

Definition 1. Feature matrix ¯X is τ-distinguishable, if for a positive constant τ ,

Clearly, feature clustering is possible if and only if ¯X is τdistinguishable for τ large enough. In the following theorem, we prove that the pairwise clustering step correctly groups features in X if ¯X is sufficiently τ-distinguishable.

Theorem 1. When ¯X is at least τ-distinguishable with

then with probability at least 1-2rδ, ∀i ∈ {1, 2, · · · , r}, the cluster center ˆXi is an inconsistent estimator of ¯Xi ,

Proof. First we prove that the columns in X can be separated when noise ϵ is not too large. For ¯Xt and its cluster index set Ct, we require

With probability at least 1 - 2rδ,

And

Therefore, Eq. (12) is satisfied when

From concentration,

Theorem 1 claims that when different features in ¯X is dissimilar to each other, we can correctly cluster the feature matrix X. However, the estimated cluster center is not consistent. The squared bias of the cluster center estimation is on order of O(σ2 i /ki +σ2/√n). Even n is infinity, this bias will never converge to zero. Theorem 1 shows that it is impossible to consistently estimate ¯X, no matter how many training instances are given.

As conventional sparse analysis, we assume that ¯X satisfies R.I.P. condition. Note that our definition is slightly different from conventional one in order to simplify the notation.

Definition 2 (δs-RIP). A matrix ¯X is δs-R.I.P. if there is a constant 0 < δs < 1, for any s sparse vector w whose support set is F,

Remark 1 If ¯X is δs-R.I.P. , it must be τ-distinguishable. To see this, take any two 1-sparse vector ei, ej where ei is the unit vector along the i-th coordinate. When s ≥ 2, we have

Therefore ¯X is at least 2√1- δs distinguishable. The converse is not true. Clearly ¯XF can be well separated but is low rank. Here we want to emphasize that we should treat the constants δs and τ independently. This is because although ¯X is at least 2√1- δs distinguishable, its optimal τ may be significantly larger than 2√1 - δs.

Following Definition 2, it is easy to check the next lemma. We omit the proof to save space.

Lemma 2. Assume that u, v are s sparse vectors and ¯X is δ2s-R.I.P. . Denote set F = supp(u)∪supp(v), then for any positive η ≤ 1/(1 - δ2s),

where ρ = 1 - η(1 - δ2s).

Lemma 2 claims that the gradient descent operator is restricted isometric when ¯X is R.I.P. . Lemma 2 can be extended to a more generalized version of R.I.P. condition, as shown in [Xiaotong Yuan, 2014]. Lemma 2 is critical in the convergence analysis because when ρ < 1, Eq. (13) is a contraction map which indicates a geometrical convergence rate.

With the help of the above discussion, we give the convergence rate of Algorithm 1.

Theorem 2. Under the same assumptions in Theorem 1 and ¯X is δ3s-R.I.P. . Algorithm 1 converges to ¯w∗ geometrically,

where

provide that 0 ≤ ρ < 1 .

Theorem 2 claims that the convergence rate of Algorithm 1 is linear. The ρ controls the exponential convergence rate which depends on the step size and R.I.P. condition. ∆3s is the bias term due to clustering which converges to zero at speed O(1/√n). b is the estimation error due to the noise in label y and the uncertainty in the clustering step. We see that when the noise level ξ is zero and n is sufficiently large, the estimation error b converges to zero, which indicates the consistency of Algorithm 1.

The novelty of Theorem 2 is its consistency. It is easy to verify that without density correction, sparse recovery methods based on P1 and P2 can only estimate ¯w∗ up to accuracy O(σi/√ki), no matter how many training instances are given. By density correction, we improve the accuracy to O(σi/√nki). This is significantly better than traditional methods when n is large enough. However, this power does not come free. In R.I.P. sparse recovery, the sampling complexity is O(s log d). From Theorem 2, in non-R.I.P. sparse recovery, our sampling complexity is O(s2 log d), which is slightly larger than the R.I.P. one. This extra s in the sampling complexity comes from the uncertainty of clustering step. In our proof, we prove that with a fix support set we can bound the gradient after density correction and this bound holds for any set no larger than s. This step brings us an extra s in sampling complexity. At the time when we write this paper, we are not aware of any published results showing a sharper bound.

We give a sketch of our proof of Theorem 2. First we show that in Eq. (7) the gradient descent after correction is consistent.

Lemma 3. For any w supported on set F, |F| ≤ s, with probability at least 1-δ,

where

Proof. For a fix set F, similar to standard online learning analysis,

From matrix concentration Lemma 1, with probability at least 1 - δ,

Now we consider the all possible set |F| ≤ s in d dimensional space. Clearly there are no more than (ds) ≤ (ed/s)s such sets. We have

Combining all the above together, the proof is done.

In Lemma 3, ∆s is the bias of the corrected gradient ∇F ˆ ℓ(w) + [D(ν)]Fw. ∆s converges to zero at speed O(1/√n). The complete proof of Theorem 1 follows immediately. Our proof is based on the proof given in [Yuan et al., 2013, Theorem 2 part b]. It is worth to mention that their model cannot deal with non-R.I.P. feature matrix.

proof of Theorem 2. Denote F = S(t) ∪ S(t-1) ∪ S∗ . So |F| ≤ 3s. We try to bound estimation accuracy of the 3s sparse vector u = [ ˆwt]F . We have

Fro R.I.P. condition,

(1) comes from Lemma 3. Because St is the largest s elements in ˆwt, we have

Then we get the recursive inequality about wt. Therefore we have

where

To ensure the convergence, we require ρ < 1.

Since y is generated by model Eq. (1),

6 Conclusion

We propose a density correction sparse recovery algorithm for non-R.I.P. problems. We first show that the conventional methods cannot consistently recover the sparse vector due to feature correlation and cluster center uncertainty. Then we propose the density correction that combines with hard iterative thresholding to recover the sparse vector consistently. The proposed algorithm has geometrical convergence rate. It adaptively removes the variance of the bias term in the clustered feature matrix.

The density correction has potentially more applications other than sparse recovery. It shows that directly learning on the clustered features may not be optimal. The cluster center is usually biased therefore the estimation is usually inconsistent. The density correction improves the estimation to be consistent without estimating the bias term itself, which could be rather difficult in many applications. Instead it exploits the cluster structure to estimate the second order statistics: the variance of the bias term. Considering the fact that clustering is widely used in machine learning, the density correction might be applied in various circumstances.

We mainly focus on least square loss function in this paper. For general convex loss function, it is much more difficult to derive the density correction. A key difficulty is that the expectation of the bias term in the gradient is no longer the cluster density. We need more sophisticated statistics to estimate the variance of this bias term. For example, recent development in convex total least square [Malioutov and Slavov, 2014] extends the total least square estimation to general convex loss functions. We open this topic for future research.

Acknowledgement

This paper was partially supported by the Intelligence Advanced Research Projects Activity (IARPA) via Department of Interior National Business Center contract number D11PC20068 and by the National Science Foundation under Grant No. IIS-1251187. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. Disclaimer: The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of IARPA, DOI/NBC, the National Science Foundation or the U.S. Government.