Optimal Shrinkage of Singular Values

Introduction

For example, when choosing the square Frobenius loss, or mean square error (MSE)

where $X$ and $\hat{X}$ are $m$-by-$n$ matrices, we would like to find an estimator $\hat{X}$ with small mean square error (MSE). The default technique for estimating a low rank matrix in noise is the Truncated SVD (TSVD) : write

where $r=rank(X)$, assumed known, and $y_{1}\geq\ldots\geq y_{m}$. Being the best approximation of rank $r$ to the data in the least squares sense , and therefore the Maximum Likelihood estimator when $Z$ has Gaussian entries, the TSVD is arguably as ubiquitous in science and engineering as linear regression .

The TSVD estimator shrinks to zero some of the data singular values, while leaving others untouched. More generally, for any specific choice of scalar nonlinearity $\eta:[0,\infty)\to[0,\infty)$, also known as a shrinker, there is a corresponding singular value shrinkage estimator $\hat{X}_{\eta}$ given by

For scalar and vector denoising, univariate shrinkage rules have proved to be simple and practical denoising methods, with near-optimal performance guarantees under various performance measures . Shrinkage makes sense for singular values, too: presumably, the observed singular values $y_{1}\ldots y_{m}$ are “inflated” by the noise, and applying a carefully chosen shrinkage function, one can obtain a good estimate of the original signal $X$.

Indeed, there is a growing body of literature on matrix denoising by shrinkage of singular values, going back, to the best of our knowledge, to Owen and Perry and Shabalin and Nobel . Soft thresholding of singular values has been considered in , and hard thresholding in . In fact, and, very recently, considered shrinkers that are developed specifically for singular values, and measured their performance using Frobenius loss.

These developments suggest the following question: Is there a simple, natural shrinkage nonlinearity for singular values? If there is a simple answer to this question, surely it depends on the loss function $L$ and on specific assumptions on the signal matrix $X$.

In we have performed a narrow investigation that focused on hard and soft thresholding of singular values under the Frobenius loss (1). We adopted a simple asymptotic framework that models the situation where $X$ is low-rank, originally proposed in and inspired by Johnstone’s Spiked Covariance Model . In this framework, the signal matrix dimensions $m=m_{n}$ and $n$ both to infinity, such that their ratio converges to an asymptotic aspect ratio: $m_{n}/n\to\beta$, with $0<\beta\leq 1$, while the column span of the signal matrix remains fixed. Building on a recent probabilistic analysis of this framework we have discovered that, in this framework, there is an asymptotically unique admissible threshold for singular values, in the sense that it offers equal or better asymptotic MSE to that of any other threshold choice, no matter which specific low-rank model may be in force.

The main discovery reported here is that this phenomenon is in fact much more general: in this asymptotic framework, which models low-rank matrices observed in white noise, for each of a variety of loss functions, there exists a single asymptotically unique admissible shrinkage nonlinearity, in the sense that it offers equal or better asymptotic loss than any other shrinkage nonlinearity, at each specific low-rank model that can occur. In other words, once the loss function has been decided, in a definite asymptotic sense, there is a single rational choice of shrinkage nonlinearity.

In this paper, we develop a general method for finding the optimal shrinkage nonlinearity for a variety of loss functions. We explicitly work out the optimal shrinkage formula for the Frobenius norm loss, the nuclear norm loss, and the operator norm loss. Let us denote the Frobenius, Operator and Nuclear matrix norms by $\left|\left|\cdot\right|\right|_{F}$,$\left|\left|\cdot\right|\right|_{op}$ and $\left|\left|\cdot\right|\right|_{*}$, respectively. If the singular values of the matrix $X-\hat{X}$ are $\sigma_{1},\ldots,\sigma_{m}$, then these losses are given by

As we will see, the optimal nonlinearity for the Frobenius norm loss (4), in a natural noise scaling, is

In the asymptotically square case $\beta=1$ this reduces to

Optimal shrinker for Operator norm loss.

The operator norm loss (5) for matrix estimation has mostly been studied in the context of covariance estimation . Let us define

As we will see, the optimal nonlinearity for operator loss is just

Optimal Shrinkage for Nuclear norm loss.

The Nuclear norm loss (6) has also been proposed for matrix estimation. See and references within for discussion of the Nuclear norm and, more generally, of Schatten-$p$ norms as losses for matrix estimation.

As we will see, the optimal nonlinearity for nuclear norm loss is

where $x=x(y)$ is given in (8). Note that the formulas above are calibrated for the natural noise level $\sigma=1/\sqrt{n}$; see Section 8.1 below for usage in known noise level $\sigma$ or unknown noise level. In the code supplement for this paper we offer a Matlab implementation of each of these shrinkers in known or unknown noise.

Figure 2 shows the three nonlinearities (7), (9) and (10). As we will see, these nonlinearities, and many others that are not calculated explicitly in this paper, flow from a single general method for calculating optimal nonlinearities, developed here.

1 Optimal shrinkers vs. hard and soft thresholding

The optimal shrinkers presented have simple, closed-form formulas. Yet there are shrinkage rules that are simpler still, namely, hard and soft thresholding. These nonlinearities are extremely popular for scalar and vector denoising, due to their simplicity and various optimality properties . Recall that for $y\geq 0$,

It is worthwhile to ask how our optimal shrinkers differ, in shape and performance, from the popular hard and soft thresholding. To make a comparison, one should first decide how to tune the thresholds $\lambda$ and $s$. In our asymtotic framework, fortunately, there is a decisive answer to the tuning question: in previous work , we have restricted our attention to hard and soft thresholding under the Frobenius loss (4). It was shown that there exist optimal values $\lambda_{*}(\beta)$ and $s_{*}(\beta)$, which are unique admissible in the sense that they offer asymptotic performance equal to or better than the performance of any other thresold. The optimal thresholds are given by

where again $\beta$ is the limiting aspect ratio, $m_{n}/n\to\beta$.

Consider, for example, the square matrix case $\beta=1$. Under the MSE loss, the optimal hard threshold is then $\lambda_{*}=4/\sqrt{3}$, and the optimal soft threshold is $s_{*}=2$. Figure 2 shows the nonlinearities $\eta^{hard}_{\lambda_{*}}$ and $\eta^{soft}_{s_{*}}$ against our optimal shrinkers (7), (9) and (10). In high SNR ($y\gg 1$) the optimal shrinkers agree with hard thresholding and neither performs any shrinkage, while soft thresholding shrinks even strong signals. As shown in , the worst-case asymptotic MSE over a rank-$r$ matrix observed in noise level $1/\sqrt{n}$ is $2r$ for our optimal shrinker (7), $3r$ for the optimally tuned hard thresholding nonlinearity $\eta^{hard}_{\lambda_{*}}$ and $6r$ for the optimally tuned soft thresholding nonlinearity $\eta^{soft}_{s_{*}}$. Hard thresholding is worse in intermediate SNR levels; Soft thresholding is worse in strong SNR. For further discussion on this phenomenon, which stems from the random rotation of the data singular vectors due to noise, see . We conclude that optimal shrinkage, developed in this paper, offers significant performance improvement over hard and soft thresholding - even when they are optimally tuned.

Preliminaries

In the general model $Y=X+\sigma Z$, the noise level in the singular values of $Y$ is $\sqrt{n}\sigma$. Instead of specifying a different shrinkage rule that depends on the matrix size $n$, we calibrate our shrinkage rules to the “natural” model $Y=X+Z/\sqrt{n}$. In this convention, shrinkage rules stay the same for every value of $n$, and we conveniently abuse notation by writing $\hat{X}_{\eta}$ as in (3) for any $\hat{X}_{\eta}:M_{m\times n}\to M_{m\times n}$, keeping $m$ and $n$ implicit. To apply any denoiser $\hat{X}$ below to data from the general model $Y=X+\sigma Z$, use the denoiser

Throughout the text, we use $\hat{X}_{\eta}$ to denote singular value shrinker calibrated for noise level $1/\sqrt{n}$. In Section 8.1 below we provide a recipe for applying any denoiser $\hat{X}_{\eta}$ calibrated for noise level $\sigma=1/\sqrt{n}$ for data in the presence of unknown noise level.

2 Asymptotic framework and problem statement

In this paper, we consider a sequence of increasingly larger denoising problems

with $X_{n},Z_{n}\in M_{m_{n},n}$, satisfying the following assumptions:

Invariant white noise: The entries of $Z_{n}$ are i.i.d samples from a distribution with zero mean, unit variance and finite fourth moment. To simplify the formal statement of our results, we assume that this distribution is orthogonally invariant in the sense that $Z_{n}$ follows the same distribution as $AZ_{n}B$, for every orthogonal $A\in M_{m_{n},m_{n}}$ and $B\in M_{n,n}$. This is the case, for example, when the entries of $Z_{n}$ are Gaussian. In Section 8.2 we revisit this restriction and discuss general (not necessarily invariant) white noise.

Asymptotic aspect ratio $\beta$: The sequence $m_{n}$ is such that $m_{n}/n\to\beta$. To simplify our formulas, we assume that $0<\beta\leq 1$.

Note that while the signal rank $r$ and nonzero signal singular values $x_{1},\ldots,x_{r}$ are shared by all matrices $X_{n}$, the signal left and right singular vectors $U_{n}$ and $V_{n}$ are unknown and arbitrary. We also remark that the assumption, whereby the signal singular values are non-degenerate ($x_{i}>x_{i+1}$, $1\leq i<r$), is not necessary for our results to hold, yet it simplifies the analysis considerably.

Our results imply that the asymptotic loss $L_{\infty}$ exists and is well-defined, as a function of the signal singular values $\mathbf{x}$, for a large class of nonlinearities.

Optimal Shrinker. Let $L$ be a loss family. If a shrinker $\eta^{*}$ has an asymptotic loss that satisfies

3 Our contribution

At first glance, it seems too much to hope that optimal shrinkers in the sense of Definition 2 even exist. Indeed, existence of an optimal shrinker for a loss family $L$ implies that, asymptotically, the decision-theoretic picture is extremely simple and actionable: from the asymptotic loss perspective, there is a single rational choice for shrinker.

In our current terminology, Shabalin and Nobel have effectively shown that an optimal shrinker exists for Frobenius loss. The estimator they derive can be shown to be equivalent to the optimal shrinker (7), yet was given in a more complicated form. (In Section 4 we visit the special case of Frobenius loss in detail, and prove that (7) is the optimal shrinker.)

Our contribution in this paper is as follows.

We rigorously establish the existence of an optimal shrinker for a variety of loss families, including the popular Frobenius, operator and nuclear norm losses.

We provide a framework for finding the optimal shrinkers for a variety of loss families including these popular losses. As discussed in Section 8.1, our framework can be applied whether the noise level $\sigma$ is known or unknown.

We use our framework to find simple, explicit formulas for the optimal shrinkers for Frobenius, operator and nuclear norm losses, and show that it allows simple numerical evaluation of optimal shrinkers when a closed-form formula for the optimal shrinker is unavailable.

In the related problem of covariance estimation in the Spiked Covariance Model, in collaboration with I. Johnstone we identified a similar phenomenon, namely, existence of optimal eigenvalue shrinkers for covariance estimation .

The Asymptotic Picture

In the “null case” $X_{n}\equiv 0$, the empirical distribution of the singular values of $Y_{n}=Z_{n}/\sqrt{n}$ famously converges as $n\to\infty$ to the generalized quarter-circle distribution , whose density is

This distribution is compactly supported on $[\beta_{-},\beta_{+}]$, with

Moreover, in this null case we have $y_{n,1}\stackrel{{\scriptstyle a.s.}}{{\rightarrow}}1+\sqrt{\beta}$, see . We say that the singular values of $Y_{n}$ form a (generalized) quarter circle bulk and call $\beta_{+}$ the bulk edge.

Expanding seminal results of and many other authors, Benaych-Georges and Nadakuditi have provided a thorough analysis of a collection of models, which includes the model (12) as a special case. In this section we summarize some of their results regarding asymptotic behaviour of the model (12), which are relevant to singular value shrinkage.

Additional notation is required to state these facts formally. We rewrite the sequence of signal matrices in our asymptotic framework (13) as

Asymptotic location of the top $r$ data singular values. For $1\leq i\leq r$,

Asymptotic angle between signal and data singular vectors. Let $1\leq i\neq j\leq r$ and assume that $x_{i}\geq\beta^{1/4}$ is non-degenerate, namely, the value $x_{i}$ appears only once in $\mathbf{x}$. Then

If however $x_{i}<\beta^{1/4}$, then we have

We also note the following fact regarding the data singular values [25, proof of Theorem 2.9]:

Let $i>r$ be fixed. Then $y_{n,i}\stackrel{{\scriptstyle a.s.}}{{\rightarrow}}\beta_{+}$.

Optimal Shrinker for Frobenius Loss

As an introduction to the more general framework developed below, we first examine the Frobenius loss case, following the work of Shabalin and Nobel . Using Definition 1, let $L=\left\{L_{m,n}\right\}$ be the Frobenius loss family, namely $L_{m,n}$ is given by (4).

Directly expanding the Frobenius matrix norm, we obtain:

Frobenius loss of singular value shrinkage. For any shrinker $\eta:[0,\infty)\to[0,\infty)$, we have

This implies a lower bound on Frobenius loss of any singular value shrinker:

For any shrinker $\eta:[0,\infty)\to[0,\infty)$, we have

Combining Corollary 1, Lemma 1 and Lemma 2 we obtain a lower bound for the asymptotic Frobenius loss (see ):

For any continuous shrinker $\eta:[0,\infty)\to[0,\infty)$, we have

The notation $L_{2,2}$ in (26) will be made apparent below, see (48).

2 Optimal shrinker matching the lower bound

The singular value shrinker $\eta^{*}$, for which $\hat{X}_{\eta^{*}}$ minimizes the asymptotic lower bound, thus becomes a natural candidate for the optimal shrinker for Frobenius loss. Indeed, by definition, for $\hat{X}_{\eta^{*}}$ the limits of (23) and (24) are the smallest possible. It remains to show that the limit of (25) is the smallest possible.

It is clear from (25) that a necessary condition for a shrinker $\eta$ to be successful, let alone optimal, is that it must set to zero data eigenvalues that do not correspond to signal. With (25) in mind, we should only consider shrinkers $\eta$ for which $\eta(y)=0$ for any $y\leq\beta_{+}$. The following is a sufficient condition for a shrinker to achieve the lowest limit possible in the term (25), namely, for this term to converge to zero.

Assume that a continuous shrinker $\eta:[0,\infty)\to[0,\infty)$ satisfies $\eta(y)=0$ whenever $y\leq\beta_{+}+\varepsilon$ for some fixed $\varepsilon>0$. We say that $\eta$ is a Conservative shrinker.

By Lemma 1 and Lemma 2, it is clear that conservative shrinkers set to zero all data singular values $\left\{y_{i}\right\}$ which originate from pure noise ($x_{i}=0$), as well as all data singular values $\left\{y_{i}\right\}$ which are “engulfed” in the noise bulk, rendering their corresponding singular vectors useless ($x_{i}<\beta^{1/4}$). Conservative shrinkers are so called since they leave a (possibly infinitesimally small) safety margin $\varepsilon$. They enjoy the following key property:

Let $\eta:[0,\infty)\to[0,\infty)$ be a conservative shrinker. Then

By Lemma 3 we have $y_{n,r+1}\stackrel{{\scriptstyle a.s.}}{{\rightarrow}}\beta_{+}$. Let $N$ be the (random) index such that $y_{n,r+1}<\beta_{+}+\varepsilon$ for all $n>N$. Then for all $n>N$ and all $i>r$ we have $y_{n,i}<\beta_{+}+\varepsilon$, hence $\eta(y_{n,i})=0$. The desired almost sure convergence follows. ∎

Ironically, careful inspection of the candidate (7) reveals that it is continuous yet not strictly conservative: it only satisfies $\eta(y)=0$ for $y\leq\beta_{+}$, leaving no margin above the bulk edge $\beta_{+}$. In fact, building on Lipschitz continuity of the Frobenius loss itself, it can be shown that Lemma 5 remains true for the shrinker (7) as well ; this is however outside our present scope. Consequently, the asymptotic loss of (7) matches the lower bound from Corollary 2, and is lower than the asymptotic loss of any other continuous shrinker, for any low-rank model $(x_{1},\ldots,x_{r})$.

A Framework for Finding Optimal Shrinkers

With the previous section in mind, our main result may be summarized as follows: the basic ingredients that enabled us to find the optimal shrinker for Frobenius loss allow us to find the optimal shrinker for each of a variety of loss families. For these loss families, an optimal shrinker exists and is given by a simple formula. To avoid some technical nuisance, we focus on finding the optimal shrinker among conservative shrinkers.

To get started, let us describe the loss families to which our method applies.

Orthogonally invariant loss. A loss $L_{m,n}(\cdot,\cdot)$ is orthogonally invariant if for all $m,n$ we have $L_{m,n}(A,B)=L_{m,n}(UAV,UBV)$, for any orthogonal $U\in O_{m}$ and $V\in O_{n}$.

Decomposable loss family. Let $A,B\in M_{m\times n}$ and let $m=\sum_{i=1}^{k}m_{i}$ and $n=\sum_{i=1}^{k}n_{i}$. Assume that there are matrices $A_{i},B_{i}\in M_{m_{i},n_{i}}$, $1\leq i\leq k$, such that

in the sense that $A$ and $B$ are block-diagonal with blocks $\left\{A_{i}\right\}$ and $\left\{B_{i}\right\}$, respectively. A loss family $L=\{L_{m,n}\}$ is sum-decomposable if, for all $m,n$ and $A,B$ with block diagonal structure as above,

As primary examples, we consider loss families defined in Section 1: The Frobenius norm loss $L^{fro}$, the operator norm loss $L^{op}$ and the nuclear norm loss $L^{nuc}$. It is easy to check that (i) each of these losses are orthogonally invariant, and (ii) the families $L^{fro}$ and $L^{nuc}$ are sum-decomposable, while the family $L^{op}$ is max-decomposable. Our framework for finding optimal shrinkers can now be stated as follows.

Characterization of the optimal singular value shrinker. Let

and suppose that for any $x\geq\beta^{1/4}$ there exists a unique minimizer

such that $\eta^{**}$ is a conservative shrinker on $[\beta^{1/4},\infty)$. Further suppose that there exists a point $x_{0}\geq\beta^{1/4}$ such that

where $x(y)$ is defined in Eq. (8). Then for any conservative shrinker $\eta$, the asymptotic losses $L_{\infty}(\eta^{*}|\cdot)$ and $L_{\infty}(\eta|\cdot)$ exist, and

1 Discussion

Before we proceed to prove Theorem 1, we review the information it encodes about the problem at hand and its operational meaning. Theorem 1 is based on a few simple observations:

First, if $L$ is a sum– (resp. max–) decomposable family of orthogonally invariant losses, and if $\eta$ is a conservative shrinker, then the asymptotic loss $L_{\infty}(\eta|\mathbf{x})$ at $\mathbf{x}=(x_{1},\ldots,x_{r})$ can be written as a sum (resp. a maximum) over $r$ terms. These terms have identical functional form. When $x_{i}\geq\beta^{1/4}$, these terms have the form $L_{2,2}(A(x_{i}),B(\eta,x_{i}))$, and when $0\leq x_{i}<\beta^{1/4}$, these terms have the form $L_{1,1}(x_{i},0)+L_{1,1}(0,\eta)$ (resp. $\max\{L_{1,1}(x_{i},0)\,,\,L_{1,1}(0,\eta)\}$). As a result, one finds that the zero shrinker $\eta\equiv 0$ is necessarily optimal for $0\leq x\leq\beta^{1/4}$. For $x\geq\beta^{1/4}$, one just needs to minimize the loss of a specific $2$-by-$2$ matrix, namely the function $F$ from (29), to obtain the shrinker $\eta^{**}$ of (30).

Second, the asymptotic loss curve $L_{\infty}(\eta^{**}|x)$ necessarily crosses the asymptotic loss curve of the zero shrinker $L_{\infty}(\eta\equiv 0|x)$ at a point we will denote by $x_{0}$, with $x_{0}\geq\beta^{1/4}$.

Finally, by concatenating the zero shrinker and the shrinker $\eta^{**}$ precisely at the point $x_{0}$ where their asymptotic losses cross, one obtains a shrinker which is continuous ($x_{0}>\beta^{1/4}$) or possibly discontinuous ($x_{0}=\beta^{1/4}$). However, this shrinker always has a well-defined asymptotic loss. This loss dominates the asymptotic loss of any conservative shrinker.

For some loss families $L=\left\{L_{m,n}\right\}$, it is possible to find an explicit formula for the optimal shrinker using the following steps:

Write down an explicit expression for the function $F(\eta,x)$ from (29).

Explicitly solve for the minimizer $\eta^{**}(x)$ from (30).

Write down an explicit expression for the minimum $F(\eta^{**}(x),x)$.

Solve (31) for the crossing point $x_{0}$.

Compose $\eta^{**}(x)$ with the transformation $x(y)$ from (8) to obtain an explicit form of the optimal shrinker $\eta^{*}(y)$ from (32).

In Sections 6 and 7 we offer examples of this process: in Section 6 we follow it analytically and derive simple, explicit formulae of the optimal shrinkers for the Frobenius, operator and nuclear norm losses. In Section 7 we follow it numerically and compute the optimal shrinker for any Schatten-$p$ norm loss.

In the remainder of this section we describe a sequence of constructions and lemmas leading to the proof of Theorem 1.

2 Simultaneous Block Diagonalization

Let us start by considering a fixed signal matrix and noise matrix, without placing them in a sequence. To allow a gentle exposition of the main ideas, we initially make two simplifying assumptions: first, that $r=1$, namely that $X$ is rank-$1$, and second, that $\eta$ shrinks to zero all but the first singular values of $Y$, namely, $\eta(y_{i})=0$, $i>1$. Let $X\in M_{m\times n}$ be a signal matrix and let $Y=X+Z/\sqrt{n}\in M_{m\times n}$ be a corresponding data matrix. Denote their SVD by

Thus, if $L=\left\{L_{m,n}\right\}$ is a sum- or max-decomposable family of orthogonally invariant functions, we have

A similar argument gives a similar statement for rank-$r$ matrix $X$ with non-degenerate singular values:

where $R_{n}$ is a sequence of $2r$-by-$2r$ matrices such that

The lemma follows by permuting the coordinates, and then using the invariance and the decomposability properties of the loss family $L$.

3 Deterministic formula for the asymptotic loss

In Section 5.2 we analyzed a single matrix and shown that, for fixed $m$ and $n$, the loss $L_{m,n}(X\,,\,\hat{X}_{\eta}(Y))$ decomposes to “atomic” units of the form

Let us now return to the sequence model $Y_{n}=X_{n}+Z_{n}/\sqrt{n}$ and find the limiting value of these “atomic” units as $n\to\infty$. This will lead to a simple formula for the asymptotic loss $L_{\infty}(\eta|\mathbf{x})$.

for $i=1,\ldots,r$. Combining Lemma 7, Lemma 1 and Lemma 2 we obtain:

Let $Y_{n}=X_{n}+Z_{n}/\sqrt{n}$ be a matrix sequence in our asymptotic framework with signal singular values $\mathbf{x}=(x_{1},\ldots,x_{r})$. Assume that $\eta$ is continuous at $y(x_{i})$ for some fixed $1\leq i\leq r$. If $\beta^{1/4}\leq x_{i}$ then

where $B(\eta,x)$ is given by (28), while if $0\leq x_{i}<\beta^{1/4}$ then

As a result, we now obtain the asymptotic loss $L_{\infty}$ as a deterministic function of the nonzero signal singular values $x_{1},\ldots,x_{r}$. Observe that by Lemma 3, if $\eta$ is a conservative shrinker, then eventually $\eta(y_{n,i})=0$ for all $i>r$. Therefore the assumption $\eta(y_{i})=0$ for $i>r$, required for Lemma 7, is satisfied eventually. Combining Lemma 7 and Lemma 8, we obtain

A formula for the asymptotic loss of a conservative shrinker. Assume that $L=\left\{L_{m,n}\right\}$ is a sum- or max- decomposable family of orthogonally invariant losses. Extend the definition of $B(\eta,x)$ from (28) by setting $B(\eta,x)=diag(0,\eta)$ for $0\leq x<\beta^{1/4}$. If $\eta:[0,\infty)\to[0,\infty)$ is a conservative shrinker, then

The final step toward the proof of Theorem 1 involves the case when the shrinker $\eta$ is given as a special concatenation of two conservative shrinkers. Even if the two parts of $\eta$ do not match, forming a discontinuity point in which the limits from the left and from the right disagree, we may still have a formula for the asymptotic loss – provided that the loss functions match.

Assume that there exist a point $0<x^{*}$ and two shrinkers, $\eta_{1}:[0,x^{*})\to[0,\infty)$ and $\eta_{2}:[x^{*},\infty)\to[0,\infty)$, such that

We say that the asymptotic loss functions of $\eta_{1}$ and $\eta_{2}$ cross at $x^{*}$.

A formula for the asymptotic loss of a concatenation of two conservative shrinkers. Assume that $L=\left\{L_{m,n}\right\}$ is a sum- or max- decomposable family of orthogonally invariant losses. Extend the definition of $B(\eta,x)$ from (28) by setting $B(\eta,x)=diag(0,\eta)$ for $0\leq x<\beta^{1/4}$. Assume that there exist two shrinkers, $\eta_{1}:[0,x^{*})\to[0,\infty)$ and $\eta_{2}:[x^{*},\infty)\to[0,\infty)$, whose asymptotic loss functions cross at some point $0<x^{*}$. Define

Then $L_{\infty}(\eta|\cdot)$ exists and is given by (41) if $L$ is sum-decomposable, or (42) if $L$ is max-decomposable.

Consider the shrinker $\eta_{1}\equiv 0$. By Lemma 8, $\eta_{1}$ dominates any other conservative shrinker when $0\leq x<\beta^{1/4}$. By assumption, there exists a point $\beta^{1/4}\leq x_{0}$ such that $\eta_{1}$ also dominates any conservative shrinker on $[\beta^{1/4},x_{0})$, and such that $\eta^{**}$ dominates any other conservative shrinker on $[x_{0},\infty)$. Finally, by assumption, the asymptotic loss functions of $\eta_{1}$ and $\eta^{**}$ cross at $x_{0}$. By Lemma 10, the concatenated shrinker $\eta^{*}$ dominates any conservative shrinker on $[0,\infty)$. ∎

Finding Optimal Shrinkers Analytically: Frobenius, Operator & Nuclear Losses

Theorem 1 provides a general recipe for finding optimal singular value shrinkers, which was provided in Section 5.1. To see it in action, we turn to our three primary examples, namely, the Frobenius norm loss, the operator norm loss and the nuclear norm loss. In this section we find explicit formulas for the optimal singular value shrinkers in each of these losses.

We will need the following lemmas regarding $2$-by-$2$ matrices (see ):

The eigenvalues of any $2$-by-$2$ matrix $M$ with trace $trace(M)$ and determinant $det(M)$ are given by

These are the roots of the characteristic polynomial of $M$. ∎

Let $\Delta$ be a $2$-by-$2$ matrix with singular values $\sigma_{+}>\sigma_{-}>0$. Define $t=trace(\Delta\Delta^{\prime})=\left|\left|\Delta\right|\right|_{F}^{2}$, $d=det(\Delta)$ and $r^{2}=t^{2}-4d^{2}$. Assume that $\Delta$ depends on a parameter $\eta$ and let $\dot{\sigma}_{\pm}$, $\dot{t}$ and $\dot{d}$ denote the derivative of these quantities w.r.t the parameter $\eta$. Then

By Lemma 11 we have $2\sigma_{\pm}^{2}=t\pm r$ and therefore

Differentiating and expanding $\dot{\sigma}_{+}\pm\dot{\sigma}_{-}$ we obtain the relation

and the singular values $\sigma_{+}>\sigma_{-}$ of $\Delta$ are given by

Theorem 1 allows us to rediscover the optimal shrinker for Frobenius norm loss, which was derived from first principles in Section 4. To this end, observe that by (45) we have

2 Operator norm loss

The optimal shrinker for operator norm loss $\eta^{*}(y)=x(y)$ simply shrinks the data singular value back to the ”original” location of its corresponding signal singular value.

3 Nuclear norm loss

we find that only zero of $\partial(\sigma_{+}+\sigma_{-})/\partial\eta$ occurs when $\partial t/\partial\eta-\partial d/\partial\eta=0$, namely at

recovering the optimal shrinker (9). Inspection of (9) reveals that this optimal shrinker is in fact a conservative shrinker.

Finding Optimal Shrinkers Numerically: Schatten norm losses

In Section 6 we have followed the recipe discussed in Section 5.1 analytically, and explicitly solved for the optimal shrinkers of the Frobenius, Operator and Nuclear norm losses. In some cases, the optimization problem (30) does not admit a closed-form solution, and in other cases, the closed-form solution is unreasonably complicated. For such cases, we note that it is extremely easy to solve the problem (30) numerically, as it only involves minimization of a univariate function that depends on the two eigenvalues of a $2$-by-$2$ matrix. To demonstrate that our recipe for finding optimal shrinkers can be easily executed numerically, rather than analytically, in this section we find the optimal shrinker for any Schatten-$p$ norm loss numericallyWe thank the anonymous referee for this helpful suggestion., for any value $p>0$.

where the matrix size $m,n$ has been suppressed in the notation for simplicity.

Schatten-$p$ norms and quasi-norms have been considered in the literature for matrix estimation: see and references therein. (The case $0<p\leq 1$ is of special interest in matrix completion problems due to its low-rank inducing behavior.) So far in this paper we have carefully studied three special cases: $L_{m,n}^{fro}\equiv L_{m,n}^{S_{2}}$, $L_{m,n}^{nuc}\equiv L_{m,n}^{S_{1}}$ and $L_{m,n}^{op}\equiv L_{m,n}^{S_{\infty}}$. Observe that for any $0<p<\infty$, the Schatten-$p$ loss is orthogonally invariant and sum-decomposable, hence amenable to the our analysis.

While it is in principle possible to derive the optimal shrinker for the Schatten-$p$ loss analytically using Lemma 12 and Lemma 13, the result would be a very complicated expression. Instead, we follow the recipe of Section 5.1 numerically: We select points of interest $\left\{y_{i}\right\}$ in which we would like to evaluate the optimal shrinker $\eta^{*}(y)$. We define $x_{i}=x(y_{i})$ where $y\mapsto x(y)$ is the transformation from (8). For each of the values $\{x_{i}\}$ we form a symbolic expression for the function $F(\eta,x_{i})$ from (29), and minimize it numerically to obtain the minimizer $\eta^{**}(x_{i})$ from (30). The desired value of the optimal shrinker $\eta^{*}(y_{i})$ is then given by $\eta^{*}(y_{i})=\eta^{**}(x(y_{i}))=\eta^{**}(x_{i})$.

Figure 4 and Figure 5 show the optimal shrinker discovered numerically for the Schatten-$p$ loss, for a few values of $p$. Figure 4 focuses on the case $p\geq 1$, where the Schatten-$p$ loss is given by a norm. Note the familiar shapes for the values $p=1,2,10000$ (the latter is indistinguishable from the case $p=\infty$, namely the operator norm). It seems that the optimal shrinker for all cases $1\leq p<\infty$ are continuous, and that the discontinuity found analytically for the case $p=\infty$ forms only in the limit $p\to\infty$. Figure 5 focuses on the case $0<p<1$, where the Schatten-$p$ loss is given by a quasi-norm. The numerical findings are fascinating and prompt further research: for instance, while the optimal shrinker for $p=1$ is continuous, at an unknown value $p<1$ the shrinkers become discontinuous, with a discontinuity resembling that of the $p=\infty$ case. Furthermore, for small values of $p$, the optimal shrinkers are very similar to the hard thresholding nonlinearities, with a “hard threshold” that depends on $p$ and on the aspect ratio $\beta$. In other words, in these cases, the optimal shrinker and the optimal hard thresholding nonlinearity seem to approximately coincide. It also seems that as $p\to 0$, the optimal shrinkers tend to the zero shrinker. All these phenomena can be studied and evaluated precisely in further research using the framework developed in this paper.

Extensions

Our main results have been formulated and calibrated specifically for the model $Y=X+Z/\sqrt{n}\in M_{m\times n}$, where the distribution of the noise matrix $Z$ is orthogonally invariant. In this section we extend our main results to include the model $Y=X+\sigma Z\in M_{m\times n}$, and consider:

The setting where $\sigma$ is either known but does not necessarily equal $1/\sqrt{n}$, or is altogether unknown.

The setting where the noise matrix $Z$ has i.i.d entries, but its distribution is not necessarily orthogonally invariant.

Consider an asymptotic framework slightly more general than the one in Section 2.2, in which $Y_{n}=X_{n}+(\sigma/\sqrt{n})Z_{n}$, with $X_{n}$ and $Z_{n}$ as defined there. In this section we keep the loss family $L$ and the asymptotic aspect ratio $\beta$ fixed and implicit. We extend Definition 1 and write

When the noise level $\sigma$ is known, Eq. (11) allows us to re-calibrate any nonlinearity $\eta$, originally calibrated for noise level $1/\sqrt{n}$, to a different noise level. For a nonlinearity $\eta:[0,\infty)\to[0,\infty)$, write

If $\eta^{*}$ is an optimal shrinker for $Y_{n}=X_{n}+Z_{n}/\sqrt{n}$, namely,

When the noise level $\sigma$ is unknown, we are required to estimate it. See and references therein for existing literature on this estimation problem. The method below has been proposed in .

Consider the following robust estimator for the parameter $\sigma$ in the model $Y=X+\sigma Z$:

where $y_{med}$ is a median singular value of $Y$ and $\mu_{\beta}$ is the median of the Marcenko-Pastur distribution, namely, the unique solution in $\beta_{-}\leq x\leq\beta_{+}$ to the equation

where $\beta_{\pm}=1\pm\sqrt{\beta}$. Note that the median $\mu_{\beta}$ is not available analytically but can easily be obtained by numerical quadrature.

Let $\sigma>0$. For the sequence $Y_{n}=X_{n}+(\sigma/\sqrt{n})Z_{n}$ in our asymptotic framework,

Let $Y_{n}=X_{n}+(\sigma/\sqrt{n})Z_{n}$ be a sequence in our asymptotic framework and let $\eta^{*}$ be an optimal shrinker calibrated for $Y_{n}=X_{n}+Z_{n}/\sqrt{n}$. Then the random sequence of shrinkers $\eta^{*}_{\hat{\sigma}(Y_{n})}$ converges to the optimal shrinker $\eta^{*}_{\sigma}$:

Consequently, $\eta^{*}_{\hat{\sigma}(Y_{n})}$ asymptotically achieves optimal performance:

In practice, for denoising a matrix $Y\in M_{m\times n}$, assumed to satisfy $Y=X+\sigma Z$, where $X$ is low-rank and $Z$ has i.i.d entries, we have the following approximately optimal singular value shrinkage estimator:

when $\sigma$ is unknown. Here, $\eta^{*}$ is an optimal shrinker with respect to desired loss family $L$ in the natural scaling.

2 General white noise

Our results were formally stated for the sequence of models of the form $Y=X+\sigma Z$, where $X$ is a non-random matrix to be estimated, and the entries of $Z$ are i.i.d samples from a distribution that is orthogonally invariant (in the sense that the matrix $Z$ follows the same distribution as $AZB$, for any orthogonal $A\in M_{m,m}$ and $B\in M_{n,n}$). While Gaussian noise is orthogonally invariant, many common distributions, which one could consider to model white observation noise, are not.

The singular values of a signal matrix $X$ constitute a very widely used measure of the complexity, or information content, of $X$. In particular, they capture its rank. One attractive feature of the framework we adopt is that the loss $L_{m,n}(X,\hat{X})$ only depends on the signal matrix $X$ through its nonzero singular values $\mathbf{x}$. This allows the loss to be directly related to the complexity of the signal $X$. If the distribution of $Z$ is not orthogonally invariant, the loss no longer enjoys this property. This point is discussed extensively in .

In general white noise, which is not necessarily orthogonally invariant, one can still allow the loss to depend on $X$ only through its singular values by placing a prior distribution on $X$ and shifting to a model where it is a random, instead of a fixed, matrix. Specifically, consider an alternative asymptotic framework to the one in Section 2.2, in which the sequence denoising problems $Y_{n}=X_{n}+Z_{n}/\sqrt{n}$ satisfies the following assumptions:

General white noise: The entries of $Z_{n}$ are i.i.d samples from a distribution with zero mean, unit variance and finite fourth moment.

is a singular value decomposition of $X_{n}$, where $U_{n}$ and $V_{n}$ are uniformly distributed random orthogonal matrices. Formally, $U_{n}$ and $V_{n}$ are sampled from the Haar distribution on the $m$-by-$m$ and $n$-by-$n$ orthogonal group, respectively.

Asymptotic aspect ratio $\beta$: The sequence $m_{n}$ is such that $m_{n}/n\to\beta$.

The second assumption above implies that $X_{n}$ is a “generic” choice of matrix with nonzero singular values $\mathbf{x}$, or equivalently, a generic choice of coordinate systems in which the linear operator corresponding to $X$ is expressed.

The results of , which we have used, hold in this case as well. It follows that Lemma 1 and Lemma 2, and consequently all our main results, hold under this alternative framework. In short, in general white noise, all our results hold if one is willing to only specify the signal singular values, rather than the signal matrix, and consider a “generic” signal matrix with these singular values.

Simulation

Our results are exact only in the limit as the matrix size grows to infinity. To study the accuracy of the asymptotic loss on finite matrices, and to compare the optimal shrinker with optimally tuned hard and soft thresholding, we conducted two simulation studies.

We studied $n$-by-$n$ matrices of the form $Y=X+Z$. The signal matrix had exactly $r$ identical nonzero singular values. For brevity, we focused on the asymptotic Frobenius loss. Figure 6 compares the case $(n,r)=(20,1)$ with the case $(n,r)=(100,1)$. Figure 7 compares the case $(n,r)=(50,2)$ with the case $(n,r)=(50,4)$. In each case we show three different noise distributions: the entries of the noise matrix $Z$ are i.i.d draws from a Gaussian distribution (thin tails), uniform distribution (no tails) and Student-t with 6 degrees of freedom (fat tails). We overlay the predicted asymptotic loss from Eq. (41) and the observed loss for different values of the signal singular value $x$. The observed loss was obtained by averaging 50 Monte Carlo iterations. The shrinkers shown are the optimal shrinker for Frobenius loss from Eq. (7), and the optimally tuned hard and soft thresholds as described in Section 1.1. Simulations show qualitatively that our results are useful already for relatively small matrices, and that the low-rank assumption remains valid when $r/m\leq 0.1$, say.

Comparing optimal shrinkers with a brute-force calculation of the optimal shrinkage.

We studied $20$-by-$20$ matrices of the form $Y=X+Z$. The signal matrix was rank-$1$ and the noise matrix was i.i.d Gaussian. For each of the three losses $\{$ Frobenius, nuclear, operator $\}$, we calculated the optimal shrinkers using brute-force by scanning over a grid of possible values $\eta$ and finding the value that minimized the empirical loss as calculated by averaging over $10$ monte carlo draws. Figure 8 overlays the shrinkage calculated by brute-force over the asymptotically optimal shrinkers calculated for the three losses in Section 6. Note the agreement with the asymptotic formulae already for $n=20$ and rank fraction of $1/20=0.05$.

Conclusion

We have presented a general framework for finding optimal shrinkers, either analytically or numerically, for a variety of loss functions.

Note that our general method, summarized in Theorem 1, is guaranteed to find a shrinker that is asymptotically unique admissible, or optimal, among conservative shrinkers (in the sense of Definition 3). This is an artifact of our proof method, and it is best to think of Theorem 1 as a formal machine for finding “good” shrinkers, rather than a definite summary of their optimality properties. In fact, for all three loss functions considered in this paper, the optimal shrinkers we found dominate, in asymptotic loss, a much wider class of shrinkers. In particular, for all three losses, these optimal shrinkers dominate the class of continuous shrinkers with the property that $\eta(y)=0$ for all $y\leq\beta_{+}$, namely, shrinkers that truncate data singular values below the bulk edge $\beta_{+}$. In some sense, this is the class of “reasonable” shrinkers.

The challenging issue is how to control the manner in which “null” singular values $y_{n,i}$ ($i>r$) affect the loss function. When the noise distribution is Gaussian, is possible to prove an analogy of Lemma 5, showing that the cumulative effect of these “null” singular values is negligible. To formally appeal to this fact, we are required to consider only loss functions that enjoy a Lipschitz regularity property (on top of being decomposable and orthogonally invariant). Then one can show that the optimal shrinkers characterized in Theorem 1 dominate all “reasonable” shrinkers as above. See for more details.

Finally, we remark that closed-form solutions for the optimal shrinkers for Schatten-$p$ losses, and a generalization of our method to include Ky-Fan norms, both remain interesting problems for further study.

Reproducible Research

In the code supplement we offer a Matlab software library that includes:

A function that calculates the optimal singular value shrinkage w.r.t the Frobenius, operator and nuclear norm losses, both in known or unknown noise level.

Scripts that generate each of the figures in this paper.

Notably, the script which generates Figure 4 and Figure 5 includes an example of numerical evaluation of optimal shrinkers.

Acknowledgements

We thank Iain Johnstone for helpful comments. We also thank Amit Singer and Boaz Nadler for discussions stimulating this work, and Santiago Velasco-Forero for pointing out an error in an earlier version of the manuscript. We thank the anonymous referees for their helpful suggestions. This work was partially supported by NSF DMS-0906812 (ARRA). MG was partially supported by a William R. and Sara Hart Kimball Stanford Graduate Fellowship.

References