OptShrink: An algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage

Introduction

Techniques for low-rank signal matrix extraction from a signal-plus-noise matrix appear prominently in many statistical signal processing , machine learning , estimation and classification applications . In many applications, the low-rank approximation is the first step in an inferential process (see, for e.g. ). These techniques are necessary whenever the $n\times m$ signal-plus-noise data or measurement matrix formed by, for example lining up the $m$ samples or measurements of $n\times 1$ observation vectors alongside each other, can be modeled as

where H denotes the conjugate transpose and $u_{i}$ and $v_{i}$ are left and right “signal” singular vectors associated with singular values $\theta_{i}$ of the signal matrix

and $X$ is the noise-only matrix of random (not necessarily i.i.d.) noises. These models also arise in other graph signal processing type settings; see for example [62, Text before (9)], [63, Section V], [47, Section III.A] or the various models described in .

Relative to this model the objective is to form an estimate of the low-rank signal matrix assuming, for now, that its rank $r$ is known. The truncated singular value decomposition (SVD) plays a prominent role in a widely-used ‘optimal’ solution to a problem that is addressed by the famous Eckart-Young-Mirsky (henceforth, EYM) theorems . Specifically, if $||\cdot||_{F}$ denotes the matrix Frobenius norm then the solution to the constrained optimization problem

where $\widetilde{X}=\sum_{i}\widehat{\sigma}_{i}\widehat{u}_{i}\widehat{v}_{i}^{H}$ is the SVD of $\widetilde{X}$. This is also the maximum likelihood (ML), rank $r$ estimate when $X$ is assumed to be a matrix with i.i.d. Gaussian entries since the negative log-likelihood function is precisely the right hand side of (3). Its use is also justified in the small $n$, large $m$ (or vice versa) regime, whenever local asymptotic normality has ‘kicked in’.

A natural extension is to consider settings where the signal matrix is low rank and has some additional exploitable structure. Examples include low-rank and sparse (see the body of work on sparse principal component analysis. e.g. ), low rank and Toeplitz structured (see e.g. ), low rank and Hankel structured and low rank and nonnegative ; see for an excellent overview of these methods and additional references. As expected, by exploiting structure in the signal matrix we can improve estimation performance relative to the EYM estimator which assumes no structure besides the low-rank condition.

Here we place ourselves in the setting where no structure is assumed in the low-rank signal matrix and ask how the EYM estimator can be improved. The starting point for our investigation is the observation that as formulated in (3), the EYM estimator solves the representation problem of finding the best rank $r$ approximation of the signal-plus-noise measurement matrix. It says nothing about the denoising problem of how to best estimate the low-rank signal matrix, even though practitioners sometimes invoke it as though it does. Thus we should not expect the EYM estimator to be the optimal solution to the denoising problem.

Note that in (4), we are trying to approximate the unknown signal matrix using the singular vectors estimated from the noisy measurement matrix. Our setup is different from other weighted low-rank approximation problems considered in the literature as in , which involve weighted modifications of the problem in (3). In our formulation, setting $w_{i}=\widehat{\sigma}_{i}$ recovers the EYM estimator so that by inspecting the solution we can directly assess when and the extent to which the EYM estimator will be suboptimal.

We prove, using recent results from random matrix theory , that for a large class of noise models, which includes but goes well beyond the i.i.d. Gaussian model, we can compute $w^{{\rm opt}}$ in closed-form in the large matrix limit. The computation shows that $w^{{\rm opt}}$ depends only on (an integral transform of) the limiting singular value distribution of the noise-only matrix $X$. We then exploit this fact to develop a concrete algorithm for computing a consistent (in a sense we make precise) estimate of the limiting oracle solution directly from measurement matrix.

2. Form of the optimal shrinkage-and-thresholding operator

The analysis shows that $w^{{\rm opt}}_{i}$ takes the form of a shrinkage-and-thresholding operator (on the singular values of $\widetilde{X}$) that is completely characterized by the limiting singular value distribution of the noise-only matrix. The resulting shrinkage function is non-convex with $w^{{\rm opt}}_{i}\approx\widehat{\sigma}_{i}(1-O(1/\widehat{\sigma}_{i}^{2}))$ for large $\widehat{\sigma}_{i}$ and $w^{{\rm opt}}_{i}\to 0$ for $\widehat{\sigma}_{i}\leq b+o(1)$ where $b$ is a critical threshold that depends on the limiting noise-only singular value distribution.

The shrinkage portion of the solution arises because $\widehat{\sigma}_{i}$ is positively biased relative to $\theta_{i}$ and because the corresponding singular vectors of $\widetilde{X}$ are biased, noisy estimates of the (true) singular vectors of the latent signal matrix . The thresholding portion of the solution arises because of a phase transition in the ‘informativeness’ of the estimated singular vectors, relative to the latent singular vectors whereby for $\theta_{i}>\theta_{c}$ inner-products of the form $(\widehat{u}_{i}^{H}u_{i})\,$ and $(v_{i}^{H}\widehat{v}_{i})$ are $O(1)$ and tend to a constant, while for $\theta_{i}<\theta_{c}$, inner-products of the form $(\widehat{u}_{i}^{H}u_{i})\,$ and $(v_{i}^{H}\widehat{v}_{i})$ are $o(1)$ and tend to zero.

Our analysis of the structure of the optimal solution 1) brings into sharp focus the form of the optimal shrinkage-and-thresholding operator, 2) provides insight on why the EYM estimator is near optimal in the low noise regime but sub-optimal in the moderate to high noise regime and 3) explains why we can expect that soft thresholding (of singular value) operators with convex penalty functions (such as the nuclear norm ) that are tuned to be near-optimal in the small $\theta_{i}$ regime will be suboptimal in the large $\theta_{i}$ regime (and vice versa).

3. Mitigating the effect of rank over-estimation

It is a delightful fact that even though the optimization problem in (4) is unobservable, because it depends on the unknown matrix we are trying to estimate, the optimal solution itself is computable. We assume no structure, other than low rank, on the signal matrix; the exploitable structure is present in the ‘noise portion’ of the eigen-spectrum, i.e., the $\min(m,n)-r$ singular values of $\widetilde{X}$.

This makes contact with the important question of how to estimate $r$ in (1) so that one may distinguish the ‘signal portion’ of the eigen-spectrum from the ‘noise portion’. The problem has been completely solved for the setting where $X$ has i.i.d. Gaussian entries. In this setting, the recentering and rescaling constants that must be applied to the largest eigenvalue of $XX^{H}$ to produce the Tracy-Widom distribution can be precisely characterized and used to set the appropriate threshold; see . Recent work on the universality of this limiting distribution provides a rigorous justification for using essentially the same method in the non-Gaussian setting.

Similarly, when the columns of $X$ are i.i.d. and each column has a (non-identity) population covariance matrix with a known (limiting) eigen-distribution, then the results in facilitate computation of the appropriate threshold for distinguishing the ‘noise portion’ of the eigen-spectrum from the ‘signal portion’.

If the form of population covariance matrix is misspecified then applying the tests based on this theory will lead to an overestimation of the rank of the signal matrix. Developing robust estimators of the signal rank that “work” without having to specify the symmetry structure (e.g. i.i.d. elements, i.i.d. columns, variance profile, etc.) of the noise random matrix remains an important open problem. Such estimators will have to exploit (symmetry-independent) ‘universal’ features of the spectrum in a way that present estimators do not.

This is where the algorithm we have developed really shines. Our algorithm takes as its input an estimate of the rank of the signal matrix and returns a (re)weighted approximation that largely mitigates the effect of rank overestimation in a manner that the EYM estimate cannot. Thus, advances in robust rank estimation when used with our algorithm will lead to improved signal matrix approximation. If the rank is correctly estimated, then the algorithm will better estimate weak subspace components of the signal matrix than the EYM algorithm.

4. Contributions

Characterizing the limiting solution of (4), computing the resulting limiting squared error, quantifying the improvement relative to the EYM estimator and developing an implementable algorithm that realizes these performance gains are the main contributions of this paper. Some of the ideas in this paper were initially presented in a conference paper by the author , in the context of the i.i.d Gaussian noise setting. This version goes beyond the Gaussian setting considered there. We also treat the setting where measurement matrix has missing entries, as considered in . In addition to rigorous results, we formulate some (empirically validated and theoretically justified) conjectures for the structure of the solution for various ‘rank-regularized’ variations of (4).

In related work, Hachem et al looked at the problem of structured subspace estimation arising in the context of parameter estimation in large arrays. They propose an oracle solution [36, Equation (13), pp. 435] and analyze its first and second order performance in the context of the MUSIC direction-of-arrival estimator.

If we were to apply the ideas and techniques developed in this paper to the problem

then, we would recover a solution that corresponds to their oracle solution. Here, we consider the problem of estimating the low-rank matrix; our results and our new algorithm can be analyzed using the techniques in to provide insights on the first and second order convergence properties. We leave the extension of our techniques to the estimation of projection matrices is relatively straightforward as an exercise to the reader.

The paper is organized as follows. The setup, the main theoretical results and a new algorithm based on the theoretical analysis are presented in Section 2. Simulation results to validate the theoretical predictions and a comparison of our method to other matrix regularization methods are contained in Section 3.

Main results and a new algorithm

Let $X_{n}$ be an $n\times m$ ($n\leq m$, without loss of generalityWe choose this convention to simplify the definition of the empirical singular value distribution.) random matrix whose ordered singular values we denote by $\sigma_{1}(X_{n})\geq\cdots\geq\sigma_{n}(X_{n})$. Let $\mu_{X_{n}}$ be the empirical singular value distribution, i.e., the probability measure defined as

Assume that the probability measure $\mu_{X_{n}}$ converges almost surely weakly, as $n,m\longrightarrow\infty$, to a non-random compactly supported probability measure $\mu_{X}$ that is supported on $[a,b]$. We assume that $\sigma_{1}\overset{\textrm{a.s.}}{\longrightarrow}b$, where $\overset{\textrm{a.s.}}{\longrightarrow}$ denotes almost sure convergence. These conditions are satisfied by the model where $X_{n}$ has i.i.d. entries mean zero entries with variance $1/m$ and bounded higher order moments.

For a given $r\geq 1$, let $\theta_{1}>\cdots>\theta_{r}>0$ be deterministic non-zero real numbers, chosen independently of $n$. For every $n$, let $S_{n}$ be an $n\times m$ signal matrix having rank $r$ with its $r$ non-zero distinct singular values equal to $\theta_{1},\ldots,\theta_{r}$.

We suppose that $X_{n}$ and $S_{n}$ are independent and that $X_{n}$, the noise-only matrix is bi-unitarily invariant while the low-rank signal matrix $S_{n}$ is deterministic. Recall that a random matrix is said to be bi-orthogonally invariant (or bi-unitarily invariant) if its distribution is invariant under multiplication on the left and right by orthogonal (or unitary) matrices. Alternately, if $S_{n}$ has isotropically random right (or left) singular vectors, then $X_{n}$ need not be unitarily invariant under multiplication on the right (or left, resp.) by orthogonal or unitary matrices. Equivalently, $X_{n}$ can have deterministic right and left singular vectors while $S_{n}$ can have isotropically random left and right singular vectors and we would get the same result stated shortly.

A matrix $X_{n}$ with i.i.d. Gaussian entries satisfies these assumption; our results extend well beyond the Gaussian setting. The main advantage of modeling the noise matrices as having isotropically random singular vectors is that it allows us to characterize the solution in terms of just the (marginal) singular value distribution of the noise-only matrix instead of having to model the full joint distribution of the elements of the noise-only matrix.

Since the singular value distribution of the noise-only part can be estimated from the singular value distribution of the signal-plus-noise matrix, we can develop a concrete, data-driven algorithm, presented in Section 2.5, that can applied to real-world datasets to improve low-rank signal matrix recovery.

We observe a signal-plus-noise matrix $\widetilde{X}_{n}$ modeled as,

is given by $w^{{\rm eym}}_{i}=\widehat{\sigma}_{i}$ for $i=1,\ldots,r$. This yields the rank $r$ signal matrix estimate $\sum_{i=1}^{r}w^{{\rm eym}}_{i}\widehat{u}_{i}\widehat{v}_{i}^{H}$ which, by the EYM theorem, is also the solution to the representation problem in (3).

Consider the denoising optimization problem

2. Theoretical results

The solution to (7) exhibits the following behavior in the asymptotic regime where $n,m\to\infty$ and $n/m\to c\in[0,\infty)$. We have that for every $1\leq i\leq r$,

The emergence of the $D$ transform in the limit characterization of the EYM and optimal coefficients follows from the results in . There it was shown that, in the large matrix limit, the principal singular values and singular vectors of $\widetilde{X}$ can be completely characterized in terms of the singular values of the signal matrix and the $D$-transform of the limiting noise-only singular value distribution. This is why, in Theorem 2.1, the limiting values of $w^{{\rm eym}}_{i}$ and $w^{{\rm opt}}_{i}$ only depend on the singular values $\theta_{i}$ (or $\rho_{i}$) of the signal matrix and the limiting noise-only singular value distribution $\mu_{X}$.

The $D$-transform is the analog of the log-Fourier transform in the sense that it describes how the distribution of the singular values of the sums of ‘freely’ independent matrices are related to the distribution of the singular values of the individual matrices . In that sense it is an asymptotically sufficient statistic and hence its appearance in Theorem 2.1 is rather natural. See Section 2.5 of for additional remarks.

We now characterize the limiting squared error for the optimal, EYM and other estimators with arbitrary weights.

Assuming that for $i=1,\ldots,r$, $\theta_{i}^{2}>1/D_{\mu_{X}}(b^{+})$. Then in the asymptotic regime considered in Theorem 2.1, the squared error, defined as in (6), exhibits the following limiting behavior:

Theorem 2.2 reveals that whenever $w_{i}-w^{{\rm opt}}_{i}$ is large, we can expect a significant increase in SE relative to the optimal estimator. The next result reveals the shrinkage-and-thresholding form of the optimal estimator.

When $r=1$, let the sole non-zero singular value of $S$ be denoted by $\theta$ and assume that $D_{\mu_{X}}^{\prime}(b^{+})=-\infty$. Then in the asymptotic regime considered, we have that

where $\rho=D_{\mu_{X}}^{-1}(1/\theta^{2})$.

Theorem 2.3 shows that when $b$ (the a.s. limit of the largest noise-only singular value) is $O(1)$, we can expect an $O(1)$ decrease in SE, relative to the EYM estimator, by thresholding whenever $\theta^{2}<1/D_{\mu_{X}}(b^{+})$. Note that when $X$ is i.i.d. Gaussian with mean zero, variance $1/m$ entries, then $b=(1+\sqrt{c})$ and $D^{\prime}_{\mu_{X}}(b^{+})=-\infty$ so that these results apply. More generally, whenever $\mu_{X}$ exhibits a square-root decay at $b$ then $D^{\prime}_{\mu_{X}}(b^{+})=-\infty$ will be satisfied. Silverstein and Choi show that a large class of (non i.i.d.) Gaussian noise models will satisfy this condition.

3. The missing data with i.i.d. noise setting

We now consider the setting where $\widetilde{X}$ has missing entries so that the signal-plus-noise matrix is modeled as

and $\odot$ denotes the Hadamard or element-wise product. Consider the optimization problem

Note that here we are approximating $pS$ instead of $S$ as in (6) (so that we can use the data-driven algorithm as-is). Setting $w^{{\rm opt}}_{i}\mapsto w^{{\rm opt}}_{i}/p$ will yield a solution to the denoising problem in (6). Let $||w||_{\infty}=\max_{i}|w_{i}|$ denote the element of the vector $w$ with the maximum absolute value.

Assume that the singular vectors $u_{i}$ and $v_{i}$ in (10) satisfy a ‘low-coherence’ condition in the following sense: we suppose that there exist non-negative constants $\eta_{u}$, $C_{u}$, $\eta_{v}$ and $C_{v}$, independent of $n$, such that for $i=1,\ldots,r$

Let the elements of $X_{ij}$ be i.i.d with mean zero, variance $1/m$ and bounded higher order moments. Then the solution to (11) exhibits the following limiting behavior. We have that for $p\in(0,1]$ and $i=1,\ldots,r$

Theorem 2.4 is a statement about the optimality of the shrinkage-and-thresholding form when there are missing entries in the signal-plus-noise matrix. Note that in this case, the equivalent noise-only matrix will not bi-unitarily invariant when $X$ is non-Gaussian. The proof (see Section 6), however, reveals that it asymptotically behaves as though it does so that the results of Theorems 2.1 and 2.3 still apply. Note that as a consequence, Theorem 2.2 can applied to compute the result asymptotic squared error. After the submission of this paper, we learned of recent work by Shabalin and Nobel for the $p=1$ setting of Theorem 2.4 with i.i.d. Gaussian noise; see .

4. The asymptotic equivalence of various rank-regularized estimators

Let us define the effective rank, $r_{\rm eff}$, of the signal matrix as

Thus, the effective rank quantifies the number of singular values in the signal-plus-noise matrix $\widetilde{X}$ that are ‘informative’, i.e., reveal the existence of a low-rank signal matrix. Clearly, $r_{\rm eff}\leq r$ but $r_{\rm eff}<r$ whenever the number of singular values that separate from the right edge $b$ of the spectrum is less than the latent signal matrix rank $r$. The following conjecture formalizes their relation to the number of ‘informative’ singular vectors in the signal-plus-noise matrix.

Assume that $D_{\mu_{X}}^{\prime}(b^{+})=-\infty$ and Note that these conditions are met when $X$ has i.i.d. entries of variance $1/m$. See Theorem 2.10 of . that for fixed $\widehat{r}$,

with very high probability. Then we have that for $r_{\rm eff}<i\leq\widehat{r}$ and $j=1,\ldots r$,

with high enough probability that we can establish their almost sure convergence to zero.

We now consider the principal rank-regularized optimization problem

We characterize the structure of the optimal estimator and the resulting MSE next.

Let $\widehat{r}$ be a fixed (with $n$) estimate of $r_{\rm eff}$ and $r_{\rm eff}$ be defined as in (13). Then, in the asymptotic regime considered, assuming Conjecture 2.5 holds, we have that

where $\rho_{i}=D_{\mu_{X}}^{-1}(1/\theta_{i})$ and we set $\theta_{i}=0$ for $i>r$. Consequently,

Corollary 2.6 reveals that the optimal estimator can realize a significant improvement in performance relative to the EYM estimator whenever $b=O(1)$ and $r_{\rm eff}<r$. The corollary highlights the importance of reliably estimating $r_{\rm eff}$ instead of $r$. Now, consider the rank regularized optimization problem

For arbitrary integer $1\leq\widehat{r}\leq q$, the solution to (15) is given by

We state a conjecture on the delocalization of the bulk singular vectors and characterize the asymptotic limit of (15) next.

Define $q=\min(m,n)$. Assume that $D_{\mu_{X}}^{\prime}(b^{+})=-\infty$ and that for all $1\leq i\leq q$, where $i$ depends on $n$

with very high probability. Then, we have that for large enough $n$ and every $i_{n}>r_{\rm eff}$ and $j=1,\ldots r$

with high enough probability that we can establish their almost sure convergence to zero.

Assuming Conjectures 2.5 and 2.8 hold, we have that

Consequently, even though, for finite $n$

Corollary 2.9 shows that when there is delocalization in the singular vectors then, in the large matrix limit, optimal performance is attained by estimating the effective rank $r_{\rm eff}$, applying shrinkage to the informative $r_{\rm eff}$ components and thresholding (to zero) the remaining components. In other words, there are vanishing (with $n$) performance losses when the coefficients given by $w^{{\rm opt}}(r_{\rm eff})$ are used in place of $\overline{w}^{{\rm opt}}(q)$. We believe that Conjectures 2.5 and 2.8 hold in the signal-plus-noise matrix with missing entries setting considered in Section 2.3 so that Corollaries 2.6 and 2.9 will apply there as well. This is pertinent because we now describe an algorithm for consistently estimating $w^{{\rm opt}}$ directly from data by exploiting the information in the singular value spectrum of the signal-plus-noise matrix.

5. A new algorithm for improved denoising

Equation (LABEL:eq:wopt_thm) shows that the optimal estimator in the large matrix limit is given by

where $\rho_{i}$ is the large matrix limit of the $i$-th largest singular value. In the finite $n,m$ setting, for $i=1,\ldots,r_{\rm eff}$, $\widehat{\rho}_{i}=\widehat{\sigma}_{i}$ is a biased, but asymptotically consistent estimator of $\rho_{i}$. We now describe an algorithm for estimating $w^{{\rm opt}}_{i}$ using a single signal-plus-noise matrix.

By construction (and the definition of the $D$-transform), $\widehat{D}(z;X)\overset{\textrm{a.s.}}{\longrightarrow}D_{\mu_{X}}(z)$ and $\widehat{D}^{\prime}(z;X)\overset{\textrm{a.s.}}{\longrightarrow}D^{\prime}_{\mu_{X}}(z)$ for $z$ outside the support of $\mu_{X}$. We now show how the spectrum of $\widetilde{X}$ can be used to estimate $\mu_{X}$. To that end, we establish a useful identify by first defining

Then, it is easy to see that for fixed (with $n$) $\widehat{r}$, $\mu_{X_{\widehat{r}}}\overset{\textrm{a.s.}}{\longrightarrow}\mu_{X,0}$. Thus, if

is a diagonal matrix containing the $q-\widehat{r}$ “noise” singular values of $\widetilde{X}$, then, by construction, and whenever $\widehat{\sigma}_{i}\overset{\textrm{a.s.}}{\longrightarrow}\rho_{i}>b$, then $\widehat{D}(\widehat{\sigma}_{i};\widehat{\Sigma}_{\widehat{r}})\overset{\textrm{a.s.}}{\longrightarrow}D_{\mu_{X}}(\rho_{i})$ and $\widehat{D}^{\prime}(\widehat{\sigma}_{i};\widehat{\Sigma}_{\widehat{r}})\overset{\textrm{a.s.}}{\longrightarrow}D^{\prime}_{\mu_{X}}(\rho_{i}).$ Hence, we form a consistent estimate of $w^{{\rm opt}}_{i}$ as described in Algorithm 1. The methods described in Section 1.3 can be used to form an estimate of $\widehat{r}$.

By Theorem 2.2, we can compute an estimate of the absolute and relative mean squared error (defined as $\textrm{MSE}/||S||_{F}^{2}$) as

respectively. A value for $\textrm{rel}\widehat{\textrm{MSE}}_{\widehat{r}}$ near indicates very good low-rank signal matrix approximation while a value near $1$ indicates a poor approximation. These metrics might be better proxies for the noisiness of a signal-plus-noise matrix than the condition number or the spectral gap. We conclude with a statement of the theoretical consistency of the $w^{{\rm opt}}_{i}$ produced by Algorithm 1.

Assume that $\widehat{r}=r_{\rm eff}$. Then for $1\leq i\leq\widehat{r}$, we have that

This is a straightforward consequence of Theorem 2.1-a) and the fact that the almost sure limit of (16a) leads (as described in the introduction of ) directly to the $D$-transform. ∎

Numerical Validation, Discussion and Extensions

We now numerically validate our predictions. In the experiments that follow, we consider the model in (1) with $r=1$, $n=m=400$ and select $X$ to be an $n\times m$ matrix with i.i.d. $\mathcal{N}(0,1/m)$ entries. For various values of $\theta$, Figure 1-a) compares empirically computed $w^{{\rm opt}}_{1}$ averaged over $100$ trials with the (limiting) theoretical prediction given by the $p=1$ result in Theorem 2.4. Figure 1-b) compares the realized normalized MSE and shows that the EYM solution is near-optimal for large values of $\theta$ but far from sub-optimal for small values of $\theta$. The simulations validate the shrinkage-and-thresholding form of the solution for $w^{{\rm opt}}_{1}$ given by Theorem 2.3 and show that Algorithm 1 realizes the predicted performance gains.

We now consider the optimization problem in (14) and evaluate the performance of the various algorithms for various values of $\widehat{r}$ for $\theta=10$ and $\theta=2$. Here, $r_{\rm eff}=1$ and Corollary 2.6 predicts that the optimal (oracle) algorithm should significantly outperform the EYM algorithm whenever $\widehat{r}>r_{\rm eff}$. Figure 2 shows the validity of this prediction and also shows that even though Algorithm 1 is suboptimal, relative to the oracle estimator, it is able to largely mitigate the effect of $r_{\rm eff}$ overestimation due to the shrinkage effect.

Figure 3 compares the normalized MSE estimates as a function of $\theta$, produced by Algorithm 1 to the empirical values for the setting where $\widehat{r}=r=1$ and $\theta_{1}=\theta$ and where $\widehat{r}=r=2$, $\theta_{1}=20$ and $\theta_{2}=\theta$. As expected the estimates, produced are accurate whenever $r_{\rm eff}=\widehat{r}$.

We now validate Theorem 2.4. We fix $r=1$ and $\theta_{1}=\theta=2$ in (10) and vary $p$, the proportion of entries with missing data. We sample $u_{1}$ and $v_{1}$ uniformly at random from the unit hypersphere so that the low-coherence conditions in Theorem 2.4 are met. Theorem 2.4 predicts that $w^{{\rm opt}}_{1}\to 0$ (asymptotically) when $p<\sqrt{n/m}/\theta^{2}=0.25$. Figure 4 shows the accuracy of the prediction and the significant improvement in performance of the oracle estimator and Algorithm 1 relative to the EYM estimator.

We now compare our algorithm to regularized matrix estimates obtained as the solution to the optimization problem

where $||\cdot||_{*}$ is the nuclear norm (or the sum of the singular values of the argument matrix). The optimization problem in (18) yields the closed-form solution

The resulting singular value thresholded (SVT) matrix corresponds to the weighting

Figure 5-a) and b) compare the resulting soft-thresholding operator associated with the SVT approximation with the optimal and the EYM solutions for $\lambda=1,2$ as a function of $\theta$ and $w^{{\rm eym}}$, respectively for the same $r=1$, $n=m$ setting in (1) with $X_{ij}$ i.i.d. $\mathcal{N}(0,1/m)$. Here $b=(1+\sqrt{c})=2$.

While SVT with $\lambda=2$ can yield comparable shrinkage (in the small $\theta$ regime) and thresholding (below $\theta=1$) as the optimal estimator, $w_{{\rm svt}}(2)-w^{{\rm opt}}$ will be large for moderate $\theta$ so that by Theorem 2.2-c) we expect SVT to be suboptimal for larger values of $\theta$. Figure 6 compares the performance of Algorithm 1 and the optimal estimator to the SVT algorithm with $\lambda=1$ and $\lambda=2$. SVT is significantly suboptimal as expected. Our results show that our algorithm would outperform SVT with convex shrinkage functions for any of the general family of noise models considered here.

2. Better singular value shrinkage with non-convex potential functions?

A closer examination of Figure 5-a) and b) reveals that the optimal estimator shrinks less for larger values of $\theta$ than the SVT possibly can. In fact, the optimal estimator will generically yield a non-convex shrinkage function which scales as

for large $\widehat{\sigma}_{i}$. Might singular value shrinkage with other non-convex potential functions generically outperform convex potential functions as well? These would be the non-convex analogs in the matrix setting of the non-negative Garrotte estimator in the vector setting. Fully understanding their benefits and shortfalls, relative to Algorithm 1, remains an open line of inquiry.

3. Role of informative components

We conclude by reexamining the role of the principal (or leading) $r_{\rm eff}$ singular vectors of $\widetilde{X}$ in the solution of the optimization problem (15). Theorem 2.7 shows that we should take the components $\widehat{u}_{i}$ and $\widehat{v}_{i}$ for which the inner product $(\widehat{u}_{i}^{H}u_{i})\,$ and $(v_{i}^{H}\widehat{v}_{i})$ is $O(1)$. The supposition in (7) is that the principal components are these components.

However, in an expository paper by the author , it is shown that if the (limiting) spectrum of the noise-only matrix is supported on two disconnected intervals, then the middle components can be more informative than the principal components. Thus, while this work (via Theorem 2.2) brings into focus the importance of accurately estimating $r_{\rm eff}$, it is equally important to be able to identify the most informative components. The development of fast, accurate algorithms for the same for large matrix-valued datasets remains an important open problem.

4. Extensions

We have initial numerical evidence that the algorithm presented here outperforms the EYM estimator for the variety of applications described in , even though they do not exactly fit the noise matrix models analyzed here. Extending the analysis of our algorithm to these models would shed further insight on the limits of low-rank signal matrix approximation.

We conclude by listing some directions of future research. These include 1) rigorously establishing the delocalization conjectures, 2) designing penalty functions that are robust to noise model mismatch, 3) clarifying the benefits, if any, of matrix regularization with convex or non-convex penalty functions relative to rank regularized solutions for the unstructured low-rank signal matrix setting, 4) extending the methods developed to problems involving estimation of signal matrices with an unstructured low-rank component and a sparse or diagonal component or low-rank structured component and 5) developing minimax estimators, along the lines of the work in , except for the more general class of noise models considered here.

Lastly, consider Theorem 2.3, where it is shown that for $\theta<1/D_{\mu_{X}}(b^{+})$, ${\rm SE}({w^{{\rm opt}}_{1}})\overset{\textrm{a.s.}}{\longrightarrow}\theta^{2}$. In this regime, is there another (non-SVD based) algorithm that can estimate the signal matrix with mean-squared-error $\theta^{2}-O(1)$? More generally, is there a non-SVD based algorithm that can (reliably) recover the (unstructured) low-rank signal matrix in the regime where the SVD based methods break down? This is a largely open question whose answer would better clarify the interplay between the limits of SVD-based estimation of the signal matrix singular vectors and the fundamental limits of estimation of the signal matrix itself. We leave these questions for future work.

Proof of Theorems 2.1, 2.3 and 2.7

We first prove Theorem 2.1 -b). Since $w^{{\rm eym}}_{i}=\widehat{\sigma}_{i}$, Theorem 2.1-b) follows immediately from Theorem 2.9 in . Next, we prove the first part of Theorem 2.1-a) by showing that

Theorem 2.7 follows by adopting the exact same approach, with some minor modifications so we shall omit its proof. We first establish some intermediate results.

where $\operatorname{diag}(\cdot)$ denotes a matrix with the arguments on the diagonal and zeros elsewhere (even for a rectangular matrix). Then

For fixed $r$, the solution to the optimization problem

Let $U_{r}=\begin{bmatrix}u_{1}&\ldots&u_{r}\end{bmatrix}$, $V_{r}=\begin{bmatrix}v_{1}&\ldots&v_{r}\end{bmatrix}$, $\Theta_{r}=\operatorname{diag}(\theta_{1},\ldots,\theta_{r})$, $\widehat{U}=\begin{bmatrix}\widehat{u}_{i}&\ldots\widehat{u}_{n}\end{bmatrix}$ and $\widehat{V}=\begin{bmatrix}\widehat{v}_{i}&\ldots\widehat{v}_{m}\end{bmatrix}$. Then for $W=\operatorname{diag}(w_{1},\ldots,w_{r},0,\ldots 0)$, the optimization problem can be rewritten as

By the unitary invariance of the Frobenius norm we have that

Let $K=\widehat{U}^{H}U_{r}\Theta_{r}V_{r}^{H}\widehat{V}$. Then,

Expanding out the diagonal entries of $K$ we get

follows immediately from (20) by the application of Corollary 4.1. We have thus proved the equality on the left-hand side of Theorem 2.1-a). It is easy to see how this approach yields Theorem 2.7.

We now prove the limit characterization portion of Theorem 2.1-a). In [5, Theorem 2.10 c)], it was proved that for $j=1,\ldots,r,$ and $i\neq j$ such that $\theta_{i}^{2}>1/D_{\mu_{X}}(b^{+})$, $\widehat{u}_{i}^{H}u_{j}\overset{\textrm{a.s.}}{\longrightarrow}0$ and $v_{j}^{H}\widehat{v}_{i}\overset{\textrm{a.s.}}{\longrightarrow}0$. Consequently,

Let $\rho_{i}=D^{-1}_{\mu_{X}}(1/\theta_{i})$. In [5, Theorem 2.10 c)] it was shown that

where $\widetilde{\mu}_{X}=c\mu_{X}+(1-c)\delta_{0}$ and for any probability measure $\mu$,

While there is ambiguity in the sign (or phase, when complex valued) of the individual singular vectors, the proof in shows that

However, $D_{\mu_{X}}(z)=\phi_{\mu}(z)\cdot\phi_{\widetilde{\mu}}(z)$ so that $\phi_{\mu}(\rho_{i})\cdot\phi_{\widetilde{\mu}}(\rho_{i})=D_{\mu_{X}}(\rho_{i})=D_{\mu_{X}}(D^{-1}_{\mu_{X}}(1/\theta_{i}^{2}))=1/\theta_{i}^{2}$, so that

This gives the limit on the right hand side of part a).

To prove part c), we note that $w^{{\rm eym}}_{i}>\theta_{i}$ (as a consequence of Horn’s interlacing inequalities ) while, for large enough $n$, $w^{{\rm opt}}_{i}=\theta_{i}\widehat{u}_{i}^{H}u_{i}v_{i}^{H}\widehat{v}_{i}+o(1)<\theta_{i}$. Thus $w^{{\rm eym}}_{i}>w^{{\rm opt}}_{i}$ for large enough $n$. Since $\widehat{u}_{i}^{H}u_{i}v_{i}^{H}\widehat{v}_{i}\to 1$ for $\theta_{i}\to\infty$, $w^{{\rm eym}}_{i}\overset{\textrm{a.s.}}{\longrightarrow}w^{{\rm opt}}_{i}$ as $\theta_{i}\to\infty$.

We now prove Theorem 2.3. Note that when $r=1$,

When $r=1$ and $\theta_{1}^{2}\leq 1/D_{\mu_{X}}(b^{+})$ and $D^{\prime}_{\mu_{X}}(b^{+})=-\infty$, then by Theorem 2.11 of , $\widehat{u}_{1}^{H}u_{1}\overset{\textrm{a.s.}}{\longrightarrow}0$ and $v_{1}^{H}\widehat{v}_{1}\overset{\textrm{a.s.}}{\longrightarrow}0$. Consequently, $w^{{\rm opt}}_{1}\overset{\textrm{a.s.}}{\longrightarrow}0$ and we have established the phase transition (or shrinkage-and-thresholding form) of $w^{{\rm opt}}_{1}$ in Theorem 2.3. The expressions for ${\rm SE}({w^{{\rm opt}}})$ and ${\rm SE}({w^{{\rm eym}}})$ are a straightforward consequence of Theorem 2.2.

Proof of Theorems 2.2 and Corollaries 2.6 and 2.9

In [5, Theorem 2.10 c)], it was proved that for $j=1,\ldots,r,$ and $i\neq j$ such that $\theta_{i}^{2}>1/D_{\mu_{X}}(b^{+})$, $\widehat{u}_{i}^{H}u_{j}\overset{\textrm{a.s.}}{\longrightarrow}0$ and $v_{j}^{H}\widehat{v}_{i}\overset{\textrm{a.s.}}{\longrightarrow}0$. Hence,

where we have substituted (23) to give us the final expression in the stated result.

Theorem 2.2-a) and b) follow from substituting the limiting values of $w^{{\rm opt}}_{i}$ and $w^{{\rm eym}}_{i}$ given by Theorem 2.1 in the derived expression. Theorem 2.2-c) follows easily by simple algebraic manipulation of the limiting expressions for ${\rm SE}({w})$ and ${\rm SE}({w^{{\rm opt}}})$. The portions of Corollaries 2.6 and 2.9 that characterize the structure of the limiting weights follows immediately from Conjecture 2.5 and Conjecture 2.8 via an application of Theorem 2.7.

We now consider the asymptotic squared error. Note that

Since we have just shown that $\overline{w}^{{\rm opt}}_{i}\overset{\textrm{a.s.}}{\longrightarrow}w^{{\rm opt}}_{i}$ for $i=1,\ldots,r_{\rm eff}$, we have

then we can conclude that ${\rm SE}({\overline{w}^{{\rm opt}}})-{\rm SE}({w^{{\rm opt}}(r_{\rm eff})})\overset{\textrm{a.s.}}{\longrightarrow}0$ and we are done. To that end, we shall utilize the claim from Conjecture 2.5 that the $o(n)$ leading coefficients of $\overline{w}^{{\rm opt}}_{i}$ corresponding to the edge (or principal) singular vectors will be bounded by $O(\log n~{}{\rm factors}/n^{1/3})$ and the claim from Conjecture 2.8 that $O(n)$ of $\overline{w}^{{\rm opt}}_{i}$ coefficients corresponding to the bulk singular vectors will be bounded by $O(\log n~{}{\rm factors}/n)$ with very high probability. This gives us

If the probability is high enough we will be able to conclude that ${\rm SE}({\overline{w}^{{\rm opt}}(r_{\rm eff})})\overset{\textrm{a.s.}}{\longrightarrow}{\rm SE}({w^{{\rm opt}}(r_{\rm eff})})$. Repeating this calculation with $w^{{\rm opt}}(\widehat{r})$ and utilizing Conjecture 2.5 gives us the expression for the asymptotic squared error in Corollary 2.6.

Proof of Theorem 2.4

where ${Z}$ is the noise-only random matrix with missing entries given by

so that, from (25), $\widetilde{X}=\overline{X}+\Delta_{S}$. Let $\overline{X}=\sum_{i}\overline{\sigma}_{i}\overline{u}_{i}\overline{v}_{i}^{H}$ be the SVD of $\overline{X}$. In lieu of (11), consider the slightly modified optimization problem

We will first show that $\overline{w}^{{\rm opt}}_{i}$ is characterized by the stated expression in Theorem 2.4. Then we will show that $\sigma_{1}(\Delta_{S})\overset{\textrm{a.s.}}{\longrightarrow}0$, which we will utilize to prove that $w^{{\rm opt}}_{i}\overset{\textrm{a.s.}}{\longrightarrow}\overline{w}^{{\rm opt}}_{i}$.

Comparing (7) to (28) reveals that the left hand side of Theorem 2.1-a) still holds except with $\theta_{i}\longmapsto p\,\theta_{i}$. Consequently,

We now establish the almost sure limit of the right hand side of (29).

where $a=\sqrt{p}(1-\sqrt{c})$ and $b=\sqrt{p}(1+\sqrt{c})$ are the end points of the support of $\mu_{Z}$. Here, $\mu_{Z}$ is the famous Marčenko-Pastur distribution . It is known , that $\sigma_{1}(Z)\overset{\textrm{a.s.}}{\longrightarrow}b=\sqrt{p}(1+\sqrt{c})$. Moreover, from the results of Bloemendal et al [8, Theorems 2.4 and 2.5], we have that for any $\{u_{i}\}_{i=1}^{r}$ and $\{v_{i}\}_{i=1}^{r}$, independent of $Z$,

where $\mu_{\widetilde{Z}}=c\mu_{Z}+(1-c)\delta_{0}$ (when $c<1$). An inspection of the proofs in reveals that the almost sure limits of these bilinear forms determine the almost sure limits of $\sigma_{i}(\overline{X})$ and $(\overline{u}_{i}^{H}u_{j})\,$ and $(v_{j}^{H}\overline{v}_{i})$ for $i=1,\ldots,r$. Equation (31) asserts that these limits are the same as the limits that we would have obtained if $Z$ were i.i.d. Gaussian (and hence bi-unitarily invariant) with matching mean and variance as the $Z$ in (26). Consequently, the almost sure limit of $\overline{w}^{{\rm opt}}_{i}$ in (29) will be the same as though $Z$ were i.i.d. Gaussian with mean zero and variance $p/m$ entries. Hence, by Theorem 2.1-b)

Computing the $D$-transform of $\mu_{Z}$ in (30) (see Example 3.1 in for the computation when $p=1$ from which the general $p$ answer can be easily deduced) gives us the pertinent expression for $\overline{w}^{{\rm opt}}_{i}$ and $\overline{w}_{i}^{\rm eym}$ which match the expressions in Theorem 2.4. The $r=1$ phase transition behavior for $\overline{w}^{{\rm opt}}_{1}$ follows from Theorem 2.3.

From the perturbation theory of singular values [37, Theorem 3.3.16-(c), pp. 178], we have that

so if we can show that $\sigma_{1}(\Delta_{S})\overset{\textrm{a.s.}}{\longrightarrow}0$ then we will have shown that $w^{{\rm eym}}_{i}\overset{\textrm{a.s.}}{\longrightarrow}\overline{w}_{i}^{\rm eym}$ and we have proved Theorem 2.4-a).

To prove that $w^{{\rm opt}}_{i}\overset{\textrm{a.s.}}{\longrightarrow}\overline{w}^{{\rm opt}}_{i}$ we need a more involved argument that requires showing that we get the same limiting behavior when $Z+\Delta_{S}$ is substituted for $Z$ in the bilinear forms on the left hand side of (31). We begin by noting that

as a consequence of the variational characterization of the largest singular value. To make further progress, we shall utilize the resolvent identityThis identity can be verified by multiplying by $(wI-B)$ on the left and $(wI-A)$ on the right of the expressions on either side of the equality. which states that

where $\Im w>0$ and $A$ and $B$ are Hermitian matrices. Applying this identity with $A=ZZ^{H}$ and $B=(Z+\Delta_{S})(Z+\Delta_{S})^{H}$ yields

Since $\sigma_{1}(AB)\leq\sigma_{1}(A)\cdot\sigma_{1}(B)$ [37, Theorem 3.3.16-(d), pp. 178] and $\sigma_{1}(A+B)\leq\sigma_{1}(A)+\sigma_{1}(B)$ [37, Theorem 3.3.16-(a), pp. 178], we have that

if $\sigma_{1}(\Delta_{S})\leq\sigma_{1}(Z)$ thus leading to the inequality

where $C$ is a universal constant (that does not depend on $n$ or $m$). This gives us

which implies that the largest singular value of a matrix is a $1$-Lipschitz function of the $nm$ entries of the matrix. Moreover, $\sigma_{1}(t\,A+(1-t)\,B)\leq t\sigma_{1}(A)+(1-t)\sigma_{1}(B)$, implying that the largest singular value is a convex, $1$-Lipschitz function. Since, by (36), the entries of the $\Delta_{S}$ are bounded, independent random variables, we can apply Talagrand’s concentration inequality (see [82, Theorem 2.1.13, pp. 73]) to obtain the tail bound

which implies, via the Borel-Cantelli lemma, that

Applying (38) to (32) yields the result that

This proves Theorem 2.4-a). Moreover, from (34), we have that

and by repeating the same argument we can show that

Using the same argument it can be shown that

Following the proofs in , the convergence of these bilinear forms implies that the almost sure limits of $\sigma_{i}(\widetilde{X})$ and $(\widehat{u}_{i}^{H}u_{j})\,$ and $(v_{j}^{H}\widehat{v}_{i})$ for $i,j=1,\ldots,r$ are identical to the almost sure limits of $\sigma_{i}(\overline{X})$ and $(\overline{u}_{i}^{H}u_{j})\,$ and $(v_{j}^{H}\overline{v}_{i})$ for $i,j=1,\ldots,r$. Consequently, $w^{{\rm opt}}_{i}\overset{\textrm{a.s.}}{\longrightarrow}\overline{w}^{{\rm opt}}_{i}$ and we have proved Theorem 2.4-b) and c).

Justification for assumptions in Conjectures 2.5 and 2.8

A key aspect (see [5, Lemma 4.1]) in rigorously proving Conjectures 2.5 and 2.8 is understanding the behavior of expressions of the form

where $z_{j}$ is a singular value of $\widetilde{X}$ but not of $X$. Let $X=U\Sigma V^{H}$ and $w=U^{H}u_{i}$. Then

When $X$ has isotropically random singular vectors, $w_{j}=O(1/n)$ with high probability so if $z_{j}\in[a,b]$ and $\max_{i}\sigma_{i}(XX^{H})-\sigma_{i+1}(XX^{H})$ is bounded with probability by $O(\log n/n)$ in the bulk and the right hand side of the above expression will get unbounded (with $n$) resulting delocalization of the associated singular vectors. When $\mu_{X}$ exhibits a square root decay at the edge, then we expect the singular values at the edge to be spaced $O(n^{-2/3})$ apart with high probability so we might delocalization via the same argument. See for an exposition of some of these issues and for recent results on the fine details of the spacing distribution of Wigner and Wishart random matrices.

References