Noisy matrix decomposition via convex relaxation: Optimal rates in high dimensions

Introduction

Problems of matrix decomposition are motivated by a variety of applications. Many classical methods for dimensionality reduction, among them factor analysis and principal components analysis (PCA), are based on estimating a low-rank matrix from data. Different forms of robust PCA can be formulated in terms of matrix decomposition using the matrix ${\Gamma^{\star}}$ to model the gross errors . Similarly, certain problems of robust covariance estimation can be described using matrix decompositions with a column/row-sparse structure, as we describe in this paper. The problem of low rank plus sparse matrix decomposition also arises in Gaussian covariance selection with hidden variables , in which case the inverse covariance of the observed vector can be decomposed as the sum of a sparse matrix with a low rank matrix. Matrix decompositions also arise in multi-task regression , which involve solving a collection of regression problems, referred to as tasks, over a common set of features. For some features, one expects their weighting to be preserved across features, which can be modeled by a low-rank constraint, whereas other features are expected to vary across tasks, which can be modeled by a sparse component . See Section 2.1 for further discussion of these motivating applications.

Most past work on the model (1) has focused on the noiseless setting ($W=0$), and for the identity observation operator (so that $\mathfrak{X}({\Theta^{\star}}+{\Gamma^{\star}})={\Theta^{\star}}+{\Gamma^{\star}}$). Chandrasekaran et al. studied the case when ${\Gamma^{\star}}$ is assumed to sparse, with a relatively small number $s\ll{d_{1}}{d_{2}}$ of non-zero entries. In the noiseless setting, they gave sufficient conditions for exact recovery for an adversarial sparsity model, meaning the non-zero positions of ${\Gamma^{\star}}$ can be arbitrary. Subsequent work by Candes et al. analyzed the same model but under an assumption of random sparsity, meaning that the non-zero positions are chosen uniformly at random. In very recent work, Xu et al. have analyzed a different model, in which the matrix ${\Gamma^{\star}}$ is assumed to be columnwise sparse, with a relatively small number $s\ll{d_{2}}$ of non-zero columns. Their analysis guaranteed approximate recovery for the low-rank matrix, in particular for the uncorrupted columns. After initial posting of this work, we became aware of recent work by Hsu et al. , who derived Frobenius norm error bounds for the case of exact elementwise sparsity. As we discuss in more detail in Section 3.4, in this special case, our bounds are based on milder conditions, and yield sharper rates for problems where the rank and sparsity scale with the dimension.

In addition, the error bounds obtained by our analysis are sharp, and cannot be improved in general. More precisely, for the case of stochastic noise matrices and the identity observation operator, we prove that the squared Frobenius errors achieved by our estimators are minimax-optimal (see Theorem 2). An interesting feature of our analysis is that, in contrast to previous work , we do not impose incoherence conditions on the singular vectors of ${\Theta^{\star}}$; rather, we control the interaction with a milder condition involving the dual norm of the regularizer. In the special case of elementwise sparsity, this dual norm enforces an upper bound on the “spikiness” of the low-rank component, and has proven useful in the related setting of noisy matrix completion . This constraint is not strong enough to guarantee identifiability of the models (and hence exact recovery in the noiseless setting), but it does provide a bound on the degree of non-identifiability. We show that this same term arises in both the upper and lower bounds on the problem of approximate recovery that is of interest in the noisy setting.

The remainder of the paper is organized as follows. In Section 2, we set up the problem in a precise way, and describe the estimators. Section 3 is devoted to the statement of our main result on achievability, as well as its various corollaries for special cases of the matrix decomposition problem. We also state a matching lower bound on the minimax error for matrix decomposition with stochastic noise. In Section 4, we provide numerical simulations that illustrate the sharpness of our theoretical predictions. Section 5 is devoted to the proofs of our results, with certain more technical aspects of the argument deferred to the appendices, and we conclude with a discussion in Section 6.

Convex relaxations and matrix decomposition

We begin with some motivating applications for the general linear observation model with noise (1).

2 Convex relaxation for noisy matrix decomposition

Given the observation model $Y=\mathfrak{X}({\Theta^{\star}}+{\Gamma^{\star}})+W$, it is natural to consider an estimator based on solving the regularized least-squares program

Here $(\lambda_{d},\mu_{d})$ are non-negative regularizer parameters, to be chosen by the user. Our theory also provides choices of these parameters that guarantee good properties of the associated estimator. Although this estimator is reasonable, it turns out that an additional constraint yields an equally simple estimator that has attractive properties, both in theory and in practice.

In order to understand the need for an additional constraint, it should be noted that without further constraints, the model (1) is unidentifiable, even in the noiseless setting ($W=0$). Indeed, as has been discussed in past work , no method can recover the components $({\Theta^{\star}},{\Gamma^{\star}})$ unless the low-rank component is “incoherent” with the matrix ${\Gamma^{\star}}$. For instance, supposing for the moment that ${\Gamma^{\star}}$ is a sparse matrix, consider a rank one matrix with $\Theta^{\star}_{11}\neq 0$, and zeros in all other positions. In this case, it is clearly impossible to disentangle ${\Theta^{\star}}$ from a sparse matrix. Past work on both matrix completion and decomposition has ruled out these types of troublesome cases via conditions on the singular vectors of the low-rank component ${\Theta^{\star}}$, and used them to derive sufficient conditions for exact recovery in the noiseless setting (see the discussion following Example 4 for more details).

In this paper, we impose a related but milder condition, previously introduced in our past work on matrix completion , with the goal of performing approximate recovery. To be clear, this condition does not guarantee identifiability, but rather provides a bound on the radius of non-identifiability. It should be noted that non-identifiability is a feature common to many high-dimensional statistical models.For instance, see the paper for discussion of non-identifiability in high-dimensional sparse regression. Moreover, in the more realistic setting of noisy observations and/or matrices that are not exactly low-rank, such approximate recovery is the best that can be expected. Indeed, one of our main contributions is to establish minimax-optimality of our rates, meaning that no algorithm can be substantially better over the matrix classes that we consider.

For a given regularizer $\mathcal{R}$, we define the quantity $\kappa_{d}(\mathcal{R}):\,=\sup_{V\neq 0}|\!|\!|V|\!|\!|_{{\operatorname{F}}}/\mathcal{R}(V)$, which measures the relation between the regularizer and the Frobenius norm. Moreover, we define the associated dual norm

More specifically, we analyze the family of estimators

subject to $\varphi_{\mathcal{R}}(\Theta)\leq\alpha$ for some fixed parameter $\alpha$.

3 Some examples

Let us consider some examples to provide intuition for specific forms of the estimator (7), and the role of the additional constraint.

With this choice, it is straightforward to verify that

and moreover, that $\kappa_{d}(\mathcal{R}^{*})=\sqrt{{d_{1}}{d_{2}}}$. Consequently, in this specific case, the general convex program (7) takes the form

The constraint involving $\|\Theta\|_{\infty}$ serves to control the “spikiness” of the low rank component, with larger settings of $\alpha$ allowing for more spiky matrices. Indeed, this type of spikiness control has proven useful in analysis of nuclear norm relaxations for noisy matrix completion . To gain intuition for the parameter $\alpha$, if we consider matrices with $|\!|\!|\Theta|\!|\!|_{{\operatorname{F}}}\approx 1$, as is appropriate to keep a constant signal-to-noise ratio in the noisy model (1), then setting $\alpha\approx 1$ allows only for matrices for which $|\Theta_{jk}|\approx 1/\sqrt{{d_{1}}{d_{2}}}$ in all entries. If we want to permit the maximally spiky matrix with all its mass in a single position, then the parameter $\alpha$ must be of the order $\sqrt{{d_{1}}{d_{2}}}$. In practice, we are interested in settings of $\alpha$ lying between these two extremes.

Other applications involve models in which ${\Gamma^{\star}}$ has a relatively small number $s\ll{d_{2}}$ of non-zero columns (or a relatively small number $s\ll{d_{1}}$ of non-zero rows). Such applications include the multi-task regression problem from Example 2, the robust covariance problem from Example 3, as well as a form of robust PCA considered by Xu et al. . In this case, it is natural to constrain $\Gamma$ via the $(2,1)$-norm regularizer

where $\Gamma_{k}$ is the $k^{th}$ column of $\Gamma$ (or the $(1,2)$-norm regularizer that enforces the analogous constraint on the rows of $\Gamma$). For this choice, it can be verified that

where $U_{k}$ denotes the $k^{th}$ column of $U$, and that $\kappa_{d}(\mathcal{R}^{*})=\sqrt{{d_{2}}}$. Consequently, in this specific case, the general convex program (7) takes the form

Main results and their consequences

The notion of decomposability is defined in terms of a pair of subspaces, which (in general) need not be orthogonal complements. Here we consider a special case of decomposability that is sufficient to cover the examples of interest in this paper:

Similarly, the columnwise $(2,1)$-norm is also decomposable with respect to appropriately defined subspaces, indexed by subsets $C\subseteq\{1,2,\ldots,{d_{2}}\}$ of column indices. Indeed, using $V_{k}$ to denote the $k^{th}$ column of the matrix $V$, define

2 Restricted strong convexity

Given a loss function, the general notion of strong convexity involves establishing a quadratic lower bound on the error in the first-order Taylor approximation . In our setting, the loss is the quadratic function $\mathcal{L}(\Omega)=\frac{1}{2}|\!|\!|Y-\mathfrak{X}(\Omega)|\!|\!|_{{\operatorname{F}}}^{2}$ (where we use $\Omega=\Theta+\Gamma$), so that the first-order Taylor series error at $\Omega$ in the direction of the matrix $\Delta$ is given by

Consequently, strong convexity is equivalent to a lower bound of the form $\frac{1}{2}\|\mathfrak{X}(\Delta)\|_{2}^{2}\geq\frac{\gamma}{2}|\!|\!|\Delta|\!|\!|_{{\operatorname{F}}}^{2}$, where $\gamma>0$ is the strong convexity constant.

Restricted strong convexity is a weaker condition that also involves a norm defined by the regularizers. In our case, for any pair $(\mu_{d},\lambda_{d})$ of positive numbers, we first define the weighted combination of the two regularizers—namely

For a given matrix $\Delta$, we can use this weighted combination to define an associated norm

corresponding to the minimum value of $\mathcal{Q}(\Theta,\Gamma)$ over all decompositions of $\Delta$Defined this way, $\Phi(\Delta)$ is the infimal-convolution of the two norms $|\!|\!|\cdot|\!|\!|_{{\operatorname{N}}}$ and $\mathcal{R}$, which is a very well studied object in convex analysis (see e.g. ).

Note that if condition (22) holds with ${\tau_{n}}=0$ and any $\gamma>0$, then we recover the usual definition of strong convexity (with respect to the Frobenius norm). In the special case of the identity operator (i.e., $\mathfrak{X}(\Theta)=\Theta$), such strong convexity does hold with $\gamma=1$. More general observation operators require different choices of the parameter $\gamma$, and also non-zero choices of the tolerance parameter ${\tau_{n}}$.

While RSC establishes a form of (approximate) identifiability in general, here the error $\Delta$ is a combination of the error in estimating ${\Theta^{\star}}$ ($\Delta^{\Theta}$) and ${\Gamma^{\star}}$ ($\Delta^{\Gamma}$). Consequently, we will need a further lower bound on $|\!|\!|\Delta|\!|\!|_{{\operatorname{F}}}$ in terms of $|\!|\!|\Delta^{\Theta}|\!|\!|_{{\operatorname{F}}}$ and $|\!|\!|\Delta^{\Gamma}|\!|\!|_{{\operatorname{F}}}$ in the proof of our main results to demonstrate the (approximate) identifiability of our model under the RSC condition 22.

3 Results for general regularizers and noise

We begin by stating a result for a general observation operator $\mathfrak{X}$, a general decomposable regularizer $\mathcal{R}$ and a general noise matrix $W$. In later subsections, we specialize this result to particular choices of observation operator, regularizers, and stochastic noise matrices. In all our results, we measure error using the squared Frobenius norm summed across both matrices

With this notation, the following result applies to the observation model $Y=\mathfrak{X}({\Gamma^{\star}}+{\Theta^{\star}})+W$, where the low-rank matrix satisfies the constraint $\varphi_{\mathcal{R}}({\Theta^{\star}})\leq\alpha$. Our upper bound on the squared Frobenius error consists of three terms

As will be clarified shortly, these three terms correspond to the errors associated with the low-rank term ($\mathcal{K}_{\Theta^{\star}}$), the sparse term ($\mathcal{K}_{\Gamma^{\star}}$), and additional error ($\mathcal{K}_{\tau_{n}}$) associated with a non-zero tolerance ${\tau_{n}}\neq 0$ in the RSC condition (22).

Suppose that the observation operator $\mathfrak{X}$ satisfies the RSC condition (22) with curvature $\gamma>0$, and a tolerance ${\tau_{n}}$ such that there exist integers $r=1,2,\ldots,\min\{{d_{1}},{d_{2}}\}$, for which

Then if we solve the convex program (7) with regularization parameters $(\lambda_{d},\mu_{d})$ satisfying

Let us make a few remarks in order to interpret the meaning of this claim.

Let us focus first on the term $\mathcal{K}_{\Theta^{\star}}$, which corresponds to the complexity of estimating the low-rank component. It is further sub-divided into two terms, with the term $\lambda_{d}^{2}\>r$ corresponding to the estimation error associated with a rank $r$ matrix, whereas the term $\lambda_{d}\>\sum_{j=r+1}^{d}\sigma_{j}({\Theta^{\star}})$ corresponds to the approximation error associated with representing ${\Theta^{\star}}$ (which might be full rank) by a matrix of rank $r$. A similar interpretation applies to the two components associated with ${\Gamma^{\star}}$, the first of which corresponds to a form of estimation error, whereas the second corresponds to a form of approximation error.

where the $\precsim$ notation indicates that we ignore constant factors.

Consider an observation operator $\mathfrak{X}$ that satisfies the RSC condition (22) with $\gamma>0$ and ${\tau_{n}}=0$. Suppose that we solve the convex program (10) with regularization parameters $(\lambda_{d},\mu_{d})$ such that

Then there are universal constants $c_{j}$ such that for any matrix pair $({\Theta^{\star}},{\Gamma^{\star}})$ with $\|{\Theta^{\star}}\|_{\infty}\leq\frac{\alpha}{\sqrt{{d_{1}}{d_{2}}}}$ and for all integers $r=1,2,\ldots,\min\{{d_{1}},{d_{2}}\}$, and $s=1,2,\ldots,({d_{1}}{d_{2}})$, we have

where $S$ is an arbitrary subset of matrix indices of cardinality at most $s$.

It is worth noting the inequality (27) corresponds to a family of upper bounds indexed by $r$ and the subset $S$. For any fixed integer $s\in\{1,2,\ldots,({d_{1}}{d_{2}})\}$, it is natural to let $S$ index the largest $s$ values (in absolute value) of ${\Gamma^{\star}}$. Moreover, the choice of the pair $(r,s)$ can be further adapted to the structure of the matrix. For instance, when ${\Theta^{\star}}$ is exactly low rank, and ${\Gamma^{\star}}$ is exactly sparse, then one natural choice is $r=\operatorname{rank}({\Theta^{\star}})$, and $s=|\operatorname{supp}({\Gamma^{\star}})|$. With this choice, both the approximation terms vanish, and Corollary 1 guarantees that any solution $(\widehat{\Theta},\widehat{\Gamma})$ of the convex program (10) satisfies

Further specializing to the case of noiseless observations ($W=0$), yields a form of approximate recovery—namely

This guarantee is weaker than the exact recovery results obtained in past work on the noiseless observation model with identity operator ; however, these papers imposed incoherence requirements on the singular vectors of the low-rank component ${\Theta^{\star}}$ that are more restrictive than the conditions of Theorem 1.

[Unimprovability for elementwise sparse model] Consider a given sparsity index $s\in\{1,2,\ldots,({d_{1}}{d_{2}})\}$, where we may assume without loss of generality that $s\leq{d_{2}}$. We then form the matrix

4.1 Results for stochastic noise matrices

Our discussion thus far has applied to general observation operators $\mathfrak{X}$, and general noise matrices $W$. More concrete results can be obtained by assuming particular forms of $\mathfrak{X}$, and that the noise matrix $W$ is stochastic. Our first stochastic result applies to the identity operator $\mathfrak{X}=I$ and a noise matrix $W$ generated with i.i.d. $N(0,\nu^{2}/({d_{1}}{d_{2}}))$ entries.To be clear, we state our results in terms of the noise scaling $\nu^{2}/({d_{1}}{d_{2}})$ since it corresponds to a model with constant signal-to-noise ratio when the Frobenius norms of ${\Theta^{\star}}$ and ${\Gamma^{\star}}$ remain bounded, independently of the dimension. The same results would hold if the noise were not rescaled, modulo the appropriate rescalings of the various terms.

Suppose $\mathfrak{X}=I$, the matrix ${\Theta^{\star}}$ has rank at most $r$ and satisfies $\|{\Theta^{\star}}\|_{\infty}\leq\frac{\alpha}{\sqrt{{d_{1}}{d_{2}}}}$, and ${\Gamma^{\star}}$ has at most $s$ non-zero entries. If the noise matrix $W$ has i.i.d. $N(0,\nu^{2}/({d_{1}}{d_{2}}))$ entries, and we solve the convex program (10) with regularization parameters

then with probability greater than 1-\exp\big{(}-2\log({d_{1}}{d_{2}})\big{)}, any optimal solution $(\widehat{\Theta},\widehat{\Gamma})$ satisfies

In the statement of this corollary, the settings of $\lambda_{d}$ and $\mu_{d}$ are based on upper bounding $\|W\|_{\infty}$ and $|\!|\!|W|\!|\!|_{{\operatorname{op}}}$, using large deviation bounds and some non-asymptotic random matrix theory. With a slightly modified argument, the bound (35) can be sharpened slightly by reducing the logarithmic term to $\log(\frac{{d_{1}}{d_{2}}}{s})$. As shown in Theorem 2 to follow in Section 3.7, this sharpened bound is minimax-optimal, meaning that no estimator (regardless of its computational complexity) can achieve much better estimates for the matrix classes and noise model given here.

It is also worth observing that both terms in the bound (35) have intuitive interpretations. Considering first the term $\mathcal{K}_{\Theta^{\star}}$, we note that the numerator term $r({d_{1}}+{d_{2}})$ is of the order of the number of free parameters in a rank $r$ matrix of dimensions ${d_{1}}\times{d_{2}}$. The multiplicative factor $\frac{\nu^{2}}{{d_{1}}{d_{2}}}$ corresponds to the noise variance in the problem. On the other hand, the term $\mathcal{K}_{\Gamma^{\star}}$ measures the complexity of estimating $s$ non-zero entries in a ${d_{1}}\times{d_{2}}$ matrix. Note that there are ${{d_{1}}{d_{2}}\choose s}$ possible subsets of size $s$, and consequently, the numerator includes a term that scales as $\log{{d_{1}}{d_{2}}\choose s}\approx s\log({d_{1}}{d_{2}})$. As before, the multiplicative pre-factor $\frac{\nu^{2}}{{d_{1}}{d_{2}}}$ corresponds to the noise variance. Finally, the second term within $\mathcal{K}_{\Gamma^{\star}}$—namely the quantity $\frac{\alpha^{2}\,s}{{d_{1}}{d_{2}}}$—arises from the non-identifiability of the model, and as discussed in Example 6, it cannot be avoided without imposing further restrictions on the pair $({\Gamma^{\star}},{\Theta^{\star}})$.

Consider the factor analysis model with $n\geq d$ samples, and regularization parameters

Then with probability greater than 1-c_{2}\exp\big{(}-c_{3}\log(d)\big{)}, any optimal solution $(\widehat{\Theta},\widehat{\Gamma})$ satisfies

We note that the condition $n\geq d$ is necessary to obtain consistent estimates in factor analysis models, even in the case with ${\Gamma^{\star}}=I_{d\times d}$ where PCA is possible (e.g., see Johnstone ). Again, the terms in the bound have a natural interpretation: since a matrix of rank $r$ in $d$ dimensions has roughly $rd$ degrees of freedom, we expect to see a term of the order $\frac{rd}{n}$. Similarly, since there are $\log{d^{2}\choose s}\approx s\log d$ subsets of size $s$ in a $d\times d$ matrix, we also expect to see a term of the order $\frac{s\log d}{n}$. Moreover, although we have stated our choices of regularization parameter in terms of $|\!|\!|\Sigma|\!|\!|_{{2}}$ and $\rho(\Sigma)$, these can be replaced by the analogous versions using the sample covariance matrix $\widehat{\Sigma}$. (By the concentration results that we establish, the population and empirical versions do not differ significantly when $n\geq d$.)

4.2 Comparison to Hsu et al. [14]

Apart from these minor differences, there are two major differences between our results, and those of Hsu et al. First of all, their analysis involves three quantities ($\alpha$, $\beta$, $\gamma$) that measure singular vector incoherence, and must satisfy a number of inequalities. In contrast, our analysis is based only on a single condition: the “spikiness” condition on the low-rank component ${\Theta^{\star}}$. As we have seen, this constraint is weaker than singular vector incoherence, and consequently, unlike the result of Hsu et al., we do not provide exact recovery guarantees for the noiseless setting. However, it is interesting to see (as shown by our analysis) that a very simple spikiness condition suffices for the approximate recovery guarantees that are of interest for noisy observation models. Given these differing assumptions, the underlying proof techniques are quite distinct, with our methods leveraging the notion of restricted strong convexity introduced by Negahban et al. .

The second (and perhaps most significant) difference is in the sharpness of the results for the noisy setting, and the permissible scalings of the rank-sparsity pair $(r,s)$. As will be clarified in Section 3.7, the rates that we establish for low-rank plus elementwise sparsity for the noisy Gaussian model (Corollary 2) are minimax-optimal up to constant factors. In contrast, the upper bounds in Theorem 3 of Hsu et al. involve the product $rs$, and hence are sub-optimal as the rank and sparsity scale. These terms appear only additively both our upper and minimax lower bounds, showing that an upper bound involving the product $rs$ is sub-optimal. Moreover, the bounds of Hsu et al. (see Section IV.D) are limited to matrix decompositions for which the rank-sparsity pair $(r,s)$ are bounded as

This bound precludes many scalings that are of interest. For instance, if the sparse component ${\Gamma^{\star}}$ has a nearly constant fraction of non-zeros (say $s\asymp\frac{{d_{1}}{d_{2}}}{\log({d_{1}})\,\log({d_{2}})}$ for concreteness), then the bound (37) restricts to ${\Theta^{\star}}$ to have constant rank. In contrast, our analysis allows for high-dimensional scaling of both the rank $r$ and sparsity $s$ simultaneously; as can be seen by inspection of Corollary 2, our Frobenius norm error goes to zero under the scalings $s\asymp\frac{{d_{1}}{d_{2}}}{\log({d_{1}})\,\log({d_{2}})}$ and $r\asymp\frac{{d_{2}}}{\log({d_{2}})}$.

4.3 Results for multi-task regression

Suppose that the matrix ${\Theta^{\star}}$ has rank at most $r$ and satisfies $\|{\Theta^{\star}}\|_{\infty}\leq\frac{\alpha}{\sqrt{{d_{1}}{d_{2}}}}$, and the matrix ${\Gamma^{\star}}$ has at most $s$ non-zero entries. If the entries of $W$ are i.i.d. $N(0,\nu^{2})$, and we solve the convex program (10) with regularization parameters

then with probability greater than 1-\exp\big{(}-2\log({d_{1}}{d_{2}})\big{)}, any optimal solution $(\widehat{\Theta},\widehat{\Gamma})$ satisfies

We see that the results presented above are analogous to those presented in Corollary 2. However, in this setting, we leverage large deviations results in order to find bounds on $\|\mathfrak{X}^{*}(W)\|_{\infty}$ and $|\!|\!|\mathfrak{X}^{*}(W)|\!|\!|_{{\operatorname{op}}}$ that hold with high probability given our observation model.

5 An alternative two-step method

In detail, let us consider the following two-step estimator:

Estimate the sparse component ${\Gamma^{\star}}$ by solving

As is well-known, this convex program has an explicit solution based on soft-thresholding the entries of $Y$.

Given the estimate $\widehat{\Gamma}$, estimate the low-rank component ${\Theta^{\star}}$ by solving the convex program

Interestingly, note that this method can be understood as the first two steps of a blockwise co-ordinate descent method for solving the convex program (10). In step (a), we fix the low-rank component, and minimize as a function of the sparse component. In step (b), we fix the sparse component, and then minimize as a function of the low-rank component. The following result that these two steps of co-ordinate descent achieve the same rates (up to constant factors) as solving the full convex program (10):

Given observations $Y$ from the model $Y={{\Theta^{\star}}+{\Gamma^{\star}}}+W$ with $\|{\Theta^{\star}}\|_{\infty}\leq\frac{\alpha}{\sqrt{{d_{1}}{d_{2}}}}$, consider the two-step procedure (40) and (41) with regularization parameters $(\lambda_{d},\mu_{d})$ such that

Then the error bound (30) from Corollary 1 holds with $\gamma=1$.

Consequently, in the special case that $\mathfrak{X}=I$, then there is no need to solve the convex program (10) to optimality; rather, two steps of co-ordinate descent are sufficient.

On the other hand, the simple two-stage method will not work for general observation operators $\mathfrak{X}$. As shown in the proof of Proposition 1, the two-step method relies critically on having the quantity $\|\mathfrak{X}({\Theta^{\star}}+W)\|_{\infty}$ be upper bounded (up to constant factors) by $\max\{\|{\Theta^{\star}}\|_{\infty},\|W\|_{\infty}\}$. By triangle inequality, this condition holds trivially when $\mathfrak{X}=I$, but can be violated by other choices of the observation operator, as illustrated by the following example.

Recall the multi-task observation model first introduced in Example 2. In Corollary 4, we showed that the general estimator (10) will recover good estimates under certain assumptions on the observation matrix. In this example, we provide an instance for which the assumptions of Corollary 4 are satisfied, but on the other hand, the two-step method will not return a good estimate.

We now verify that the conditions of Corollary 4 are satisfied. Letting $\sigma_{\min}$ and $\sigma_{\max}$ denote (respectively) the smallest and largest singular values of $X$, we have $\sigma_{\min}=1$ and $\sigma_{\max}\leq 2$. Moreover, letting $X_{j}$ denote the $j^{th}$ column of $X$, we have $\max_{j=1,\ldots,d}\|X_{j}\|_{2}\leq 2$. Consequently, if we consider rescaled observations with noise variance $\nu^{2}/d$, the conditions of Corollary 4 are all satisfied with constants (independent of dimension), so that the $M$-estimator (10) will have good performance.

where step (i) exploits Jensen’s inequality, and step (ii) uses the fact that

For any noise matrix $W$ with reasonable tail behavior, the variable $\|X({\Theta^{\star}}+W)\|_{\infty}$ will concentrate around its expectation, showing that $\|X({\Theta^{\star}}+W)\|_{\infty}$ will be larger than $\|{\Theta^{\star}}\|_{\infty}$ by an order of magnitude (factor of $\sqrt{d}$). Consequently, the two-step method will have much larger estimation error in this case. $\clubsuit$

6 Results for ∥⋅∥2,1\|\cdot\|_{2,1} regularization

Let us return again to the general Theorem 1, and illustrate some more of its consequences in application to the columnwise $(2,1)$-norm previously defined in Example 5, and methods based on solving the convex program (14). As before, specializing Theorem 1 to this decomposable regularizer yields a number of guarantees. In order to keep our presentation relatively brief, we focus here on the case of the identity observation operator $\mathfrak{X}=I$.

Suppose that we solve the convex program (14) with regularization parameters $(\lambda_{d},\mu_{d})$ such that

Then there is a universal constant $c_{1}$ such that for any matrix pair $({\Theta^{\star}},{\Gamma^{\star}})$ with $\|{\Theta^{\star}}\|_{2,\infty}\leq\frac{\alpha}{\sqrt{{d_{2}}}}$ and for all integers $r=1,2,\ldots,d$ and $s=1,2,\ldots,{d_{2}}$, we have

where $C\subseteq\{1,2,\ldots,{d_{2}}\}$ is an arbitrary subset of column indices of cardinality at most $s$.

As before, if we assume that ${\Theta^{\star}}$ has exactly rank $r$ and ${\Gamma^{\star}}$ has at most $s$ non-zero columns, then both approximation error terms in the bound (44) vanish, and we recover an upper bound of the form $|\!|\!|\widehat{\Theta}-{\Theta^{\star}}|\!|\!|_{{\operatorname{F}}}^{2}+|\!|\!|\widehat{\Gamma}-{\Gamma^{\star}}|\!|\!|_{{\operatorname{F}}}^{2}\lesssim\lambda_{d}^{2}r+\mu_{d}^{2}s$. If we further specialize to the case of exact observations ($W=0$), then Corollary 5 guarantees that

The following example shows, that given our conditions, even in the noiseless setting, no method can recover the matrices to precision more accurate than $\alpha^{2}s/{d_{2}}$.

In order to demonstrate that the term $\alpha^{2}s/{d_{2}}$ is unavoidable, it suffices to consider a slight modification of Example 6. In particular, let us define the matrix

Suppose ${\Theta^{\star}}$ has rank at most $r$ and satisfies $\|{\Theta^{\star}}\|_{2,\infty}\leq\frac{\alpha}{\sqrt{{d_{2}}}}$, and ${\Gamma^{\star}}$ has at most $s$ non-zero columns. If the noise matrix $W$ has i.i.d. $N(0,\nu^{2}/({d_{1}}{d_{2}}))$ entries, and we solve the convex program (14) with regularization parameters $\lambda_{d}=\frac{8\nu}{\sqrt{{d_{1}}}}+\frac{8\nu}{\sqrt{{d_{2}}}}$ and

then with probability greater than 1-\exp\big{(}-2\log({d_{2}})\big{)}, any optimal solution $(\widehat{\Theta},\widehat{\Gamma})$ satisfies

Note that the setting of $\lambda_{d}$ is the same as in Corollary 2, whereas the parameter $\mu_{d}$ is chosen based on upper bounding $\|W\|_{2,\infty}$, corresponding to the dual norm of the columnwise $(2,1)$-norm. With a slightly modified argument, the bound (46) can be sharpened slightly by reducing the logarithmic term to $\log(\frac{{d_{2}}}{s})$. As shown in Theorem 2 to follow in Section 3.7, this sharpened bound is minimax-optimal.

As with Corollary 2, both terms in the bound (46) are readily interpreted. The term $\mathcal{K}_{\Theta^{\star}}$ has the same interpretation, as a combination of the number of degrees of freedom in a rank $r$ matrix (that is, of the order $r({d_{1}}+{d_{2}})$) scaled by the noise variance $\frac{\nu^{2}}{{d_{1}}{d_{2}}}$. The second term $\mathcal{K}_{\Gamma^{\star}}$ has a somewhat more subtle interpretation. The problem of estimating $s$ non-zero columns embedded within a ${d_{1}}\times{d_{2}}$ matrix can be split into two sub-problems: first, the problem of estimating the $s{d_{1}}$ non-zero parameters (in Frobenius norm), and second, the problem of column subset selection—i.e., determining the location of the $s$ non-zero parameters. The estimation sub-problem yields the term $\frac{\nu^{2}s{d_{1}}}{{d_{1}}{d_{2}}}$, whereas the column subset selection sub-problem incurs a penalty involving $\log{{d_{2}}\choose s}\approx s\log{d_{2}}$, multiplied by the usual noise variance. The final term $\alpha^{2}s/{d_{2}}$ arises from the non-identifiability of the model. As discussed in Example 8, it is unavoidable without further restrictions.

We now turn to some consequences for the problem of robust covariance estimation formulated in Example 3. As seen from equation (4), the disturbance matrix in this setting can be written as a sum $({\Gamma^{\star}})^{T}+{\Gamma^{\star}}$, where ${\Gamma^{\star}}$ is a column-wise sparse matrix. Consequently, we can use a variant of the estimator (14), in which the loss function is given by $|\!|\!|Y-\{{\Theta^{\star}}+({\Gamma^{\star}})^{T}+{\Gamma^{\star}}\}|\!|\!|_{{F}}^{2}$. The following result summarizes the consequences of Theorem 1 in this setting:

Consider the problem of robust covariance estimation with $n\geq d$ samples, based on a matrix ${\Theta^{\star}}$ with rank at most $r$ that satisfies $\|{\Theta^{\star}}\|_{2,\infty}\leq\frac{\alpha}{\sqrt{d}}$, and a corrupting matrix ${\Gamma^{\star}}$ with at most $s$ rows and columns corrupted. If we solve SDP (14) with regularization parameters

then with probability greater than 1-c_{2}\exp\big{(}-c_{3}\log(d)\big{)}, any optimal solution $(\widehat{\Theta},\widehat{\Gamma})$ satisfies

Some comments about this result: with the motivation of being concrete, we have given an explicit choice (47) of the regularization parameters, involving the operator norm $|\!|\!|{\Theta^{\star}}|\!|\!|_{{\operatorname{op}}}$, but any upper bound would suffice. As with the noise variance in Corollary 6, a typical strategy would choose this pre-factor by cross-validation.

7 Lower bounds

For the case of i.i.d Gaussian noise matrices, Corollaries 2 and 6 provide results of an achievable nature, namely in guaranteeing that our estimators achieve certain Frobenius errors. In this section, we turn to the complementary question: what are the fundamental (algorithmic-independent) limits of accuracy in noisy matrix decomposition? One way in which to address such a question is by analyzing statistical minimax rates.

More formally, given some family $\mathcal{F}$ of matrices, the associated minimax error is given by

where the infimum ranges over all estimators $(\widetilde{\Theta},\widetilde{\Gamma})$ that are (measurable) functions of the data $Y$, and the supremum ranges over all pairs $({\Theta^{\star}},{\Gamma^{\star}})\in\mathcal{F}$. Here the expectation is taken over the Gaussian noise matrix $W$, under the linear observation model (1).

Given a matrix ${\Gamma^{\star}}$, we define its support set $\operatorname{supp}({\Gamma^{\star}}):\,=\{(j,k)\,\mid\,{\Gamma^{\star}}_{jk}\neq 0\}$, as well as its column support \operatorname{colsupp}({\Gamma^{\star}}):\,=\{k\,\mid\,{\Gamma^{\star}}_{k}\neq 0\big{\}}, where ${\Gamma^{\star}}_{k}$ denotes the $k^{th}$ column. Using this notation, our interest centers on the following two matrix families:

By construction, Corollaries 2 and 6 apply to the families $\mathcal{F}_{\operatorname{sp}}$ and $\mathcal{F}_{\operatorname{col}}$ respectively.

The following theorem establishes lower bounds on the minimax risks (in squared Frobenius norm) over these two families for the identity observation operator:

Consider the linear observation model (1) with identity observation operator: $\mathfrak{X}(\Theta+\Gamma)=\Theta+\Gamma$. There is a universal constant $c_{0}>0$ such that for all $\alpha\geq 32\sqrt{\log({d_{1}}{d_{2}})}$, we have

Note the agreement with the achievable rates guaranteed in Corollaries 2 and 6 respectively. (As discussed in the remarks following these corollaries, the sharpened forms of the logarithmic factors follow by a more careful analysis.) Theorem 2 shows that in terms of squared Frobenius error, the convex relaxations (10) and (14) are minimax optimal up to constant factors.

In addition, it is worth observing that although Theorem 2 is stated in the context of additive Gaussian noise, it also shows that the radius of non-identifiability (involving the parameter $\alpha$) is a fundamental limit. In particular, by setting the noise variance to zero, we see that under our milder conditions, even in the noiseless setting, no algorithm can estimate to greater accuracy than $c_{0}\frac{\alpha^{2}\,s}{{d_{1}}{d_{2}}}$, or the analogous quantity for column-sparse matrices.

Simulation results

We have implemented the $M$-estimators based on the convex programs (10) and (14), in particular by adapting first-order optimization methods due to Nesterov . In this section, we report simulation results that demonstrate the excellent agreement between our theoretical predictions and the behavior in practice. In all cases, we used square matrices ($d={d_{1}}={d_{2}}$), and a stochastic noise matrix $W$ with i.i.d. $N(0,\frac{\nu^{2}}{d^{2}})$ entries, with $\nu^{2}=1$. For any given rank $r$, we generated ${\Theta^{\star}}$ by randomly choosing the spaces of left and right singular vectors. We formed random sparse (elementwise or columnwise) matrices by choosing the positions of the non-zeros (entries or columns) uniformly at random.

Note that under this scaling, Corollary 2 predicts that the squared Frobenius error should be upper bounded as $c_{1}\gamma+c_{2}\beta\log(1/\beta)$, for some universal constants $c_{1}$, $c_{2}$. Figure 1(a) provides experimental confirmation of the accuracy of these theoretical predictions: varying $\gamma$ (with $\beta$ fixed) produces linear growth of the squared error as a function of $\gamma$. In Figure 1(b), we study the complementary scaling, with the rank ratio $\gamma$ fixed and the sparsity ratio $\beta$ varying in the interval $[.01,.1]$. Since $\beta\log(1/\beta)\approx\Theta(\beta)$ over this interval, we should expect to see roughly linear scaling. Again, the plot shows good agreement with the theoretical predictions.

Now recall the estimator (14) from Example 5, designed for estimating a low-rank matrix plus a columnwise sparse matrix. We have observed similar linear dependence on the analogs of the parameters $\gamma$ and $\beta$, as predicted by our theory. In the interests of exhibiting a different phenomenon, here we report its performance for matrices of varying dimension, in all cases with ${\Gamma^{\star}}$ having $s=3r$ non-zero columns. Figure 2(a) shows plots of squared Frobenius error versus the dimension for two choices of the rank ($r=10$ and $r=15$), and the matrix dimension varying in the range $d\in\{100:25:300\}$. As predicted by our theory, these plots decrease at the rate $1/d$. Indeed, this scaling is revealed by replotting the inverse squared error versus $d$, which produces the roughly linear plots shown in panel (b). Moreover, by comparing the relative slopes of these two curves, we see that the problem with rank $r=15$ requires roughly a dimension that is roughly $\frac{3}{2}$ larger than the problem with $r=10$ to achieve the same error. Again, this linear scaling in rank is consistent with Corollary 6.

Proofs

In this section, we provide the proofs of our main results, with the proofs of some more technical lemmas deferred to the appendices.

For the reader’s convenience, let us recall here the two assumptions on the regularization parameters:

We now turn to a lemma that deals with the behavior of the error matrices $(\widehat{\Delta}^{\Theta},\widehat{\Delta}^{\Gamma})$ when measured together using a weighted sum of the nuclear norm and regularizer $\mathcal{R}$. In order to state the following lemma, let us recall that for any positive $(\mu_{d},\lambda_{d})$, the weighted norm is defined as $\mathcal{Q}(\Theta,\Gamma):\,=|\!|\!|\Theta|\!|\!|_{{\operatorname{N}}}+\frac{\mu_{d}}{\lambda_{d}}\mathcal{R}(\Gamma)$.

With this notation, we have the following:

For any $r=1,2,\ldots,d$, there is a decomposition $\widehat{\Delta}^{\Theta}=\widehat{\Delta}^{\Theta}_{A}+\widehat{\Delta}^{\Theta}_{B}$ such that:

The difference $\mathcal{Q}({\Theta^{\star}},{\Gamma^{\star}})-\mathcal{Q}({\Theta^{\star}}+\widehat{\Delta}^{\Theta},{\Gamma^{\star}}+\widehat{\Delta}^{\Gamma})$ is upper bounded by

Under conditions (52) on $\mu_{d}$ and $\lambda_{d}$, the error matrices $\widehat{\Delta}^{\Theta}$ and $\widehat{\Delta}^{\Gamma}$ satisfy

See Appendix A for the proof of this result.

Our second lemma guarantees that the cost function $\mathcal{L}(\Theta,\Gamma)=\frac{1}{2}|\!|\!|Y-\mathfrak{X}(\Theta+\Gamma)|\!|\!|_{{\operatorname{F}}}^{2}$ is strongly convex in a restricted set of directions. In particular, if we let $\delta\mathcal{L}(\widehat{\Delta}^{\Theta},\widehat{\Delta}^{\Gamma})$ denote the error in the first-order Taylor series expansion around $({\Theta^{\star}},{\Gamma^{\star}})$, then some algebra shows that

The following lemma shows that (up to a slack term) this Taylor error is lower bounded by the squared Frobenius norm.

Under the conditions of Theorem 1, the first-order Taylor series error (55) is lower bounded by

Using Lemmas 1 and 2, we can now complete the proof of Theorem 1. By the optimality of $(\widehat{\Theta},\widehat{\Gamma})$ and the feasibility of $({\Theta^{\star}},{\Gamma^{\star}})$, we have

Recalling that $Y=\mathfrak{X}({\Theta^{\star}}+{\Gamma^{\star}})+W$, and re-arranging in terms of the errors $\widehat{\Delta}^{\Theta}=\widehat{\Theta}-{\Theta^{\star}}$ and $\widehat{\Delta}^{\Gamma}=\widehat{\Gamma}-{\Gamma^{\star}}$, we obtain

where the weighted norm $\mathcal{Q}$ was previously defined (20).

We now substitute inequality (53) from Lemma 1 into the right-hand-side of the above equation to obtain

Some algebra and an application of Hölder’s inequality and the triangle inequality allows us to obtain the upper bound

Recalling conditions (52) for $\mu_{d}$ and $\lambda_{d}$, we obtain the inequality

Using inequality (56) from Lemma 2 to lower bound the right-hand side, and then rearranging terms yields

Substituting this upper bound into equation (57) yields

Substituting the above inequality into equation (58) and rearranging the terms involving $e^{2}(\widehat{\Delta}^{\Theta},\,\widehat{\Delta}^{\Gamma})$ yields the claim.

2 Proof of Corollaries 2 and 4

from which we conclude that the stated choices of $(\lambda_{d},\mu_{d})$ are valid with high probability. Turning now to the RSC condition, we note that in the case of multivariate regression, we have

showing that the RSC condition holds with $\gamma=\sigma_{\min}^{2}$.

In order to obtain the sharper result for $X=I_{{d_{1}}\times{d_{1}}}$ in Corollary 2—in which $\log({d_{1}}{d_{2}})$ is replaced by the smaller quantity $\log(\frac{{d_{1}}{d_{2}}}{s})$— we need to be more careful in upper bounding the noise term $\langle\!\langle{W},\;{\widehat{\Delta}^{\Gamma}}\rangle\!\rangle$. We refer the reader to Appendix C.1 for details of this argument.

3 Proof of Corollary 3

For this model, the noise matrix is recentered Wishart noise—namely, $W=\frac{1}{n}\sum_{i=1}^{n}Z_{i}Z_{i}^{T}-\Sigma$, where each $Z_{i}\sim N(0,\Sigma)$. Letting $U_{i}\sim N(0,I_{d\times d})$ be i.i.d. Gaussian random vectors, we have

where the final bound holds with probability greater than $1-2\exp(-c_{1}d)$, using standard tail bounds on Gaussian random matrices . Thus, we see that the specified choice (36) of $\lambda_{d}$ is valid for Theorem 1 with high probability.

which shows that the specified choice of $\mu_{d}$ is also valid with high probability.

4 Proof of Proposition 1

Our proof of Proposition 1 is based on two lemmas, of which the first provides control on the error $\widehat{\Delta}^{\Gamma}$ in estimating the sparse component.

Under the assumptions of Proposition 1, for any subset $S$ of matrix indices of cardinality at most $s$, the sparse error $\widehat{\Delta}^{\Gamma}$ in any solution of the convex program (40) satisfies the bound

Since $\widehat{\Gamma}$ and ${\Gamma^{\star}}$ are optimal and feasible (respectively) for the convex program (40), we have

Re-writing this inequality in terms of the error $\widehat{\Delta}^{\Gamma}$ and re-arranging terms yields

where the second step is based on two applications of the triangle inequality. Now by applying Hölder’s inequality and the triangle inequality to the first term on the right-hand side, we obtain

where the final inequality follows from our stated choice (42) of the regularization parameter $\mu_{d}$. Since $\|\widehat{\Delta}^{\Gamma}_{S}\|_{1}\leq\sqrt{s}|\!|\!|\widehat{\Delta}_{S}|\!|\!|_{{\operatorname{F}}}\leq\sqrt{s}|\!|\!|\widehat{\Delta}^{\Gamma}|\!|\!|_{{\operatorname{F}}}$, the claim (59) follows with some algebra. ∎

Our second lemma provides a bound on the low-rank error $\widehat{\Delta}^{\Theta}$ in terms of the sparse matrix error $\widehat{\Delta}^{\Gamma}$.

If in addition to the conditions of Proposition 1, the sparse erorr matrix is bounded as $|\!|\!|\widehat{\Delta}^{\Gamma}|\!|\!|_{{\operatorname{F}}}\leq\delta$, then the low-rank error matrix is bounded as

As the proof of this lemma is somewhat more involved, we defer it to Appendix D. Finally, combining the low-rank bound (60) with the sparse bound (59) from Lemma 3 yields the claim of Proposition 1.

5 Proof of Corollary 6

For this corollary, we have $\mathcal{R}(\cdot)=\|\cdot\|_{2,1}$ and $\mathcal{R}^{*}(\cdot)=\|\cdot\|_{2,\infty}$. In order to establish the claim, we need to show that the conditions of Corollary 5 on the regularization pair $(\lambda_{d},\mu_{d})$ hold with high probability. The setting of $\lambda_{d}$ is the same as Corollary 2, and is valid by our earlier argument. Hence, in order to complete the proof, it remains to establish an upper bound on $\|W\|_{2,\infty}$.

Let $W_{k}$ be the $k^{th}$ column of the matrix. Noting that the function $W_{k}\mapsto\|W_{k}\|_{2}$ is Lipschitz, by concentration of measure for Gaussian Lipschitz functions , we have

As before, a sharper bound (with $\log{d_{2}}$ replaced by $\log({d_{2}}/s)$) can be obtained by a refined argument; we refer the reader to Appendix C.2 for the details.

6 Proof of Corollary 7

For this model, the noise matrix takes the form $W:\,=\frac{1}{n}\sum_{i=1}^{n}U_{i}U_{i}^{T}-{\Theta^{\star}}$, where $U_{i}\sim N(0,{\Theta^{\star}})$. Since ${\Theta^{\star}}$ is positive semidefinite with rank at most $r$, we can write

7 Proof of Theorem 2

Our lower bound proofs are based on a standard reduction from estimation to a multiway hypothesis testing problem over a packing set of matrix pairs. In particular, given a collection $\{(\Theta^{j},\Gamma^{j}),j=1,2,\ldots,M\}$ of matrix pairs contained in some family $\mathcal{F}$, we say that it forms a $\delta$-packing in Frobenius norm if, for all distinct pairs $i,j\in\{1,2,\ldots,M\}$, we have

Given such a packing set, it is a straightforward consequence of Fano’s inequality that the minimax error over $\mathcal{F}$ satisfies the lower bound

We begin by proving the lower bound (50) for matrix decompositions over the family $\mathcal{F}_{\operatorname{sp}}(r,s,\alpha)$.

Let us first establish the lower bound involving the radius of non-identifiability, namely the term scaling as $\frac{\alpha^{2}s}{{d_{1}}{d_{2}}}$ in the case of $s$-sparsity for ${\Theta^{\star}}$. Recall from Example 6 the “bad” matrix (33), which we denote here by $B^{*}$. By construction, we have $|\!|\!|B^{*}|\!|\!|_{{F}}^{2}=\frac{\alpha^{2}s}{{d_{1}}{d_{2}}}$. Using this matrix, we construct a very simple packing set with $M=4$ matrix pairs $(\Theta,\Gamma)$:

Each one of these matrix pairs $(\Theta,\Gamma)$ belongs to the set $\mathcal{F}_{\operatorname{sp}}(1,s,\alpha)$, so it can be used to establish a lower bound over this set. (Moreover, it also yields a lower bound over the sets $\mathcal{F}_{\operatorname{sp}}(r,s,\alpha)$ for $r>1$, since they are supersets.) It can also be verified that for any two distinct pairs of matrices in the set (62), they differ in squared Frobenius norm by at least $\delta^{2}=\frac{1}{2}|\!|\!|B^{*}|\!|\!|_{{F}}^{2}=\frac{1}{2}\frac{\alpha^{2}s}{{d_{1}}{d_{2}}}$. Let $J$ be a random index uniformly distributed over the four possible models in our packing set (62). By construction, for any matrix pair $(\Theta,\Gamma)$ in the packing set, we have $\Theta+\Gamma=0$. Consequently, for any one of these models, the observation matrix $Y$ is simply equal to pure noise $W$, and hence $I(Y;J)=0$. Putting together the pieces, the Fano bound (61) implies that

We now describe the construction of a packing set for lower bounding the estimation error. In this case, our construction is more subtle, based on the the Cartesian product of two components, one for the low rank matrices, and the other for the sparse matrices. For the low rank component, we re-state a slightly modified form (adapted to the setting of non-square matrices) of Lemma 2 from the paper :

For ${d_{1}},{d_{2}}\geq 10$, a tolerance $\delta>0$, and for each $r=1,2,\ldots,d$, there exists a set of ${d_{1}}\times{d_{2}}$-dimensional matrices $\{\Theta^{1},\ldots,\Theta^{M}\}$ with cardinality M\geq\frac{1}{4}\exp\big{(}\frac{r{d_{1}}}{256}+\frac{r{d_{2}}}{256}\big{)} such that each matrix has rank $r$, and moreover

As for the sparse matrices, the following result is a modification, so as to apply to the matrix setting of interest here, of Lemma 5 from the paper :

For any $\delta>0$, and for each integer $s<{d_{1}}{d_{2}}$, there exists a set of matrices $\{\Gamma^{1},\ldots,\Gamma^{N}\}$ with cardinality N\geq\exp\big{(}\frac{s}{2}\log\frac{{d_{1}}{d_{2}}-s}{s/2}\big{)} such that

and such that each $\Gamma^{j}$ has at most $s$ non-zero entries.

We now have the necessary ingredients to prove the lower bound (50). By combining Lemmas 5 and 6, we conclude that there exists a set of matrices with cardinality

where the bound (i) follows from the condition (66b). Combined with lower bound (65), we see that it suffices to choose $\delta$ such that

For ${d_{1}},{d_{2}}$ larger than a finite constant (to exclude degenerate cases), we see that the choice

for a suitably small constant $c_{0}>0$ is sufficient, thereby establishing the lower bound (50).

7.2 Lower bounds for columnwise sparsity

The lower bound (51) for columnwise follows from a similar argument. The only modifications are in the packing sets.

In order to establish a lower bound of order $\frac{\alpha^{2}s}{d_{2}}$, recall the “bad” matrix (45) from Example 8, which we denote by $B^{*}$. By construction, it has squared Frobenius norm $|\!|\!|B^{*}|\!|\!|_{{F}}^{2}=\frac{\alpha^{2}s}{{d_{2}}}$. We use it to form the packing set

Each one of these matrix pairs $(\Theta,\Gamma)$ belongs to the set $\mathcal{F}_{\operatorname{col}}(1,s,\alpha)$, so it can be used to establish a lower bound over this set. (Moreover, it also yields a lower bound over the sets $\mathcal{F}_{\operatorname{col}}(r,s,\alpha)$ for $r>1$, since they are supersets.) It can also be verified that for any two distinct pairs of matrices in the set (69), they differ in squared Frobenius norm by at least $\delta^{2}=\frac{1}{2}|\!|\!|B^{*}|\!|\!|_{{F}}^{2}=\frac{1}{2}\frac{\alpha^{2}s}{{d_{2}}}$. Consequently, the same argument as before shows that

We now describe packings for the estimation error terms. For the low-rank packing set, we need to ensure that the $(2,\infty)$-norm is controlled. From the bound (63c), we have the guarantee

The following lemma characterizes a suitable packing set for the columnwise sparse component:

For all ${d_{2}}\geq 10$ and integers $s$ in the set $\{1,2,\ldots,{d_{2}}-1\}$, there exists a family ${d_{1}}\times{d_{2}}$ matrices $\{\Gamma^{k},k=1,2,\ldots N\}$ with cardinality

and such that each $\Gamma^{j}$ has at most $s$ non-zero columns.

Using this lemma and the packing set for the low-rank component and following through the Fano construction yields the claimed lower bound (50) on the minimax error for the class $\mathcal{F}_{\operatorname{col}}(r,s,\alpha)$, which completes the proof of Theorem 2.

Discussion

In this paper, we analyzed a class of convex relaxations for solving a general class of matrix decomposition problems, in which the goal is recover a pair of matrices, based on observing a noisy contaminated version of their sum. Since the problem is ill-posed in general, it is essential to impose structure, and this paper focuses on the setting in which one matrix is approximately low-rank, and the second has a complementary form of low-dimensional structure enforced by a decomposable regularizer. Particular cases include matrices that are elementwise sparse, or columnwise sparse, and the associated matrix decomposition problems have various applications, including robust PCA, robustness in collaborative filtering, and model selection in Gauss-Markov random fields. We provided a general non-asymptotic bound on the Frobenius error of a convex relaxation based on a regularizing the least-squares loss with a combination of the nuclear norm with a decomposable regularizer. When specialized to the case of elementwise and columnwise sparsity, these estimators yield rates that are minimax-optimal up to constant factors.

Various extensions of this work are possible. We have not discussed here how our estimator would behave under a partial observation model, in which only a fraction of the entries are observed. This problem is very closely related to matrix completion, a problem for which recent work by Negahban and Wainwright shows that a form of restricted strong convexity holds with high probability. This property could be adapted to the current setting, and would allow for proving Frobenius norm error bounds on the low rank component. Finally, although this paper has focused on the case in which the first matrix component is approximately low rank, much of our theory could be applied to a more general class of matrix decomposition problems, in which the first component is penalized by a decomposable regularizer that is “complementary” to the second matrix component. It remains to explore the properties and applications of these different forms of matrix decomposition.

Acknowledgements

All three authors were partially supported by grant AFOSR-09NL184. In addition, SN and MJW were partially supported by grant NSF-CDI-0941742, and AA was partially supported a Microsoft Research Graduate Fellowship. All three authors would like to acknowledge the Banff International Research Station (BIRS) in Banff, Canada for hospitality and work facilities that stimulated and supported this collaboration.

Appendix A Proof of Lemma 1

The decomposition described in part (a) was established by Recht et al. , so that it remains to prove part (b). With the appropriate definitions, part (b) can be recovered by exploiting Lemma 1 from Negahban et al. . Their lemma applies to optimization problems of the general form

We now discuss how this lemma can be applied in our special case. Here the relevant parameters are of the form $\theta=(\Theta,\Gamma)$, and the loss function is given by

The sample size $n=d^{2}$, since we make one observation for each entry of the matrix. On the other hand, the regularizer is given by the function

coupled with the regularization parameter $\gamma_{n}=\lambda_{d}$. By assumption, the regularizer $\mathcal{R}$ is decomposable, and as shown in the paper , the nuclear norm is also decomposable. Since $\mathcal{Q}$ is simply a sum of these decomposable regularizers over separate matrices, it is also decomposable.

It remains to compute the gradient $\nabla\mathcal{L}({\Theta^{\star}},{\Gamma^{\star}})$, and evaluate the dual norm. A straightforward calculation yields that $\nabla\mathcal{L}({\Theta^{\star}},{\Gamma^{\star}})=\begin{bmatrix}W&W\end{bmatrix}^{T}$. In addition, it can be verified by standard properties of dual norms

Thus, it suffices to choose the regularization parameter such that

meaning that it suffices to have $\lambda_{d}\geq 4|\!|\!|W|\!|\!|_{{\operatorname{op}}}$, as stated in the second part of condition (52).

Appendix B Proof of Lemma 2

where the second inequality follows by the definitions (20) and (21) of $\mathcal{Q}$ and $\Phi$ respectively. We now derive a lower bound on $|\!|\!|\widehat{\Delta}^{\Theta}+\widehat{\Delta}^{\Gamma}|\!|\!|_{{F}}$, and an upper bound on $\mathcal{Q}^{2}(\widehat{\Delta}^{\Theta},\widehat{\Delta}^{\Gamma})$. Beginning with the former term, observe that

so that it suffices to upper bound $\gamma|\langle\!\langle{\widehat{\Delta}^{\Theta}},\;{\widehat{\Delta}^{\Gamma}}\rangle\!\rangle|$. By the duality of the pair $(\mathcal{R},\mathcal{R}^{*})$, we have

Now since $\widehat{\Theta}$ and ${\Theta^{\star}}$ are both feasible for the program (7) and recalling that $\widehat{\Delta}^{\Theta}=\widehat{\Theta}-{\Theta^{\star}}$, an application of triangle inequality yields

where inequality (i) follows from our choice of $\mu_{d}$. Putting together the pieces, we have shown that

Since the quantity $\lambda_{d}|\!|\!|\widehat{\Delta}^{\Theta}|\!|\!|_{{\operatorname{N}}}\geq 0$, we can write

where the latter equality follows by the definition (20) of $\mathcal{Q}$.

Next we turn to the upper bound on $\mathcal{Q}(\widehat{\Delta}^{\Theta},\widehat{\Delta}^{\Gamma})$. By the triangle inequality, we have

Furthermore, substituting in equation (53) into the above equation yields

The claim then follows by substituting the above equation into equation (73), and then substituting the result into the earlier inequality (72).

Appendix C Refinement of achievability results

In this appendix, we provide refined arguments that yield sharpened forms of Corollaries 2 and 6. These refinements yield achievable bounds that match the minimax lower bounds in Theorem 2 up to constant factors. We note that these refinements are significantly different only when the sparsity index $s$ scales as $\Theta({d_{1}}{d_{2}})$ for Corollary 2, or as $\Theta({d_{2}})$ for Corollary 6.

First, suppose that $\|\widehat{\Delta}^{\Gamma}\|_{1}\leq\sqrt{s}|\!|\!|\widehat{\Delta}^{\Gamma}|\!|\!|_{{\operatorname{F}}}$. In this case, we have the upper bound

It remains to upper bound the random variable $Z(s)$. Viewed as a function of $W$, it is a Lipschitz function with parameter $\frac{\nu}{\sqrt{{d_{1}}{d_{2}}}}$, so that

Setting $\delta^{2}=\frac{4s\nu^{2}}{{d_{1}}{d_{2}}}\log(\frac{{d_{1}}{d_{2}}}{s})$, we have

with probability greater than 1-\exp\big{(}-2s\log(\frac{{d_{1}}{d_{2}}}{s})\big{)}.

It remains to upper bound the expected value. In order to do so, we apply Theorem 5.1(ii) from Gordon et al. with $(q_{0},q_{1})=(1,2)$, $n={d_{1}}{d_{2}}$ and $t=\sqrt{s}$, thereby obtaining

With this bound, proceeding through the remainder of the proof yields the claimed rate.

Alternatively, we must have $\|\widehat{\Delta}^{\Gamma}\|_{1}>\sqrt{s}|\!|\!|\widehat{\Delta}^{\Gamma}|\!|\!|_{{\operatorname{F}}}$. In this case, we need to show that the stated choice (74) of $\mu_{d}$ satisfies $\mu_{d}\|\widehat{\Delta}^{\Gamma}\|_{1}\geq 2|\langle\!\langle{W},\;{\widehat{\Delta}^{\Gamma}}\rangle\!\rangle|$ with high probability. As can be seen from examining the proofs, this condition is sufficient to ensure that Lemma 1 and Lemma 2 all hold, as required for our analysis.

where for any radius $t>0$, we define the random variable

For each fixed $t$, the same argument as before shows that $Z(t)$ is concentrated around its expectation, and Theorem 5.1(ii) from Gordon et al. with $(q_{0},q_{1})=(1,2)$, $n={d_{1}}{d_{2}}$ yields

Setting $\delta^{2}=\frac{4t^{2}\nu^{2}}{{d_{1}}{d_{2}}}\log(\frac{{d_{1}}{d_{2}}}{s})$ in the concentration bound, we conclude that

with high probability. A standard peeling argument (e.g., ) can be used to extend this bound to a uniform one over the choice of radii $t$, so that it applies to the random one $t=\frac{\|\widehat{\Delta}^{\Gamma}\|_{1}}{|\!|\!|\widehat{\Delta}^{\Gamma}|\!|\!|_{{\operatorname{F}}}}$ of interest. (The only changes in doing such a peeling are in constant terms.) We thus conclude that

with high probability. Since $\|\widehat{\Delta}^{\Gamma}\|_{1}>\sqrt{s}|\!|\!|\widehat{\Delta}^{\Gamma}|\!|\!|_{{\operatorname{F}}}$, we have $\frac{1}{\|\widehat{\Delta}^{\Gamma}\|_{1}^{2}/|\!|\!|\widehat{\Delta}^{\Gamma}|\!|\!|_{{\operatorname{F}}}^{2}}\leq\frac{1}{s}$, and hence

with high probability. With this bound, the remainder of the proof proceeds as before. In particular, the refined choice (74) of $\mu_{d}$ is adequate.

C.2 Refinement of Corollary 6

As in the refinement of Corollary 2 from Appendix C.1, we need to be more careful in controlling the noise term $\langle\!\langle{W},\;{\widehat{\Delta}^{\Gamma}}\rangle\!\rangle$. For this corollary, we make the refined choice of regularizer

As in Appendix C.1, we split our analysis into two cases.

First, suppose that $\|\widehat{\Delta}^{\Gamma}\|_{2,1}\leq\sqrt{s}|\!|\!|\widehat{\Delta}^{\Gamma}|\!|\!|_{{\operatorname{F}}}$. In this case, we have

with probability greater than 1-\exp\big{(}-2s\log(\frac{{d_{2}}}{s})\big{)}.

It remains to upper bound the expectation. Applying the Cauchy-Schwarz inequality to each column, we have

Now the variable $V_{k}$ is zero-mean, and sub-Gaussian with parameter $\frac{\nu}{\sqrt{{d_{1}}{d_{2}}}}$, again using concentration of measure for Lipschitz functions of Gaussians . Consequently, by setting $\delta_{k}=\|\Delta_{k}\|_{2}$, we can write

Applying Theorem 5.1(ii) from Gordon et al. with $(q_{0},q_{1})=(1,2)$, $n={d_{2}}$ and $t=4\sqrt{s}$ then yields

which combined with the concentration bound (76) yields the refined claim.

Alternatively, we may assume that $\|\widehat{\Delta}^{\Gamma}\|_{2,1}>\sqrt{s}|\!|\!|\widehat{\Delta}^{\Gamma}|\!|\!|_{{\operatorname{F}}}$. In this case, we need to verify that the choice (75) $\mu_{d}$ satisfies $\mu_{d}\|\widehat{\Delta}^{\Gamma}\|_{2,1}\geq 2|\langle\!\langle{W},\;{\widehat{\Delta}^{\Gamma}}\rangle\!\rangle|$ with high probability. We have the upper bound

where for any radius $t>0$, we define the random variable

Following through the same argument as in Case 2 of Appendix C.1 yields that for any fixed $t>0$, we have

with high probability. As before, this can be extended to a uniform bound over $t$ by a peeling argument, and we conclude that

with high probability. Since $\frac{1}{\|\widehat{\Delta}^{\Gamma}\|_{2,1}^{2}/|\!|\!|\widehat{\Delta}^{\Gamma}|\!|\!|_{{\operatorname{F}}}^{2}}\leq\frac{1}{s}$ by assumption, the claim follows.

Appendix D Proof of Lemma 4

Since $\widehat{\Theta}$ and ${\Theta^{\star}}$ are optimal and feasible (respectively) for the convex program (41), we have

Recalling that $Y={\Theta^{\star}}+{\Gamma^{\star}}+W$ and re-writing in terms of the error matrices $\widehat{\Delta}^{\Gamma}=\widehat{\Gamma}-{\Gamma^{\star}}$ and $\widehat{\Delta}^{\Theta}=\widehat{\Theta}-{\Theta^{\star}}$, we find that

Expanding the Frobenius norm and reorganizing terms yields

From Lemma 1 in the paper , there exists a decomposition $\widehat{\Delta}^{\Theta}=\widehat{\Delta}^{\Theta}_{A}+\widehat{\Delta}^{\Theta}_{B}$ such that the rank of $\widehat{\Delta}^{\Theta}_{A}$ upper-bounded by $2\,r$ and

where step (i) follows by triangle inequality; step (ii) by the Cauchy-Schwarz and Hölder inequality, and our assumed bound $|\!|\!|\widehat{\Delta}^{\Gamma}|\!|\!|_{{\operatorname{F}}}\leq\delta$; and step (iii) follows by substituting $\widehat{\Delta}^{\Theta}=\widehat{\Delta}^{\Theta}_{A}+\widehat{\Delta}^{\Theta}_{B}$ and applying triangle inequality.

Since we have chosen $\lambda_{d}\geq|\!|\!|W|\!|\!|_{{\operatorname{op}}}$, we conclude that

where the second inequality follows since $|\!|\!|\widehat{\Delta}^{\Theta}_{A}|\!|\!|_{{\operatorname{N}}}\leq\sqrt{2r}|\!|\!|\widehat{\Delta}^{\Theta}_{A}|\!|\!|_{{\operatorname{F}}}\leq\sqrt{2r}|\!|\!|\widehat{\Delta}^{\Theta}|\!|\!|_{{\operatorname{F}}}$. We have thus obtained a quadratic inequality in $|\!|\!|\widehat{\Delta}^{\Theta}|\!|\!|_{{\operatorname{F}}}$, and applying the quadratic formula yields the claim.

References