Optimal Shrinkage of Singular Values

arXiv:1405.7511v3 [math.ST] 15 May 2016

Matan Gavish †

David L. Donoho ∗

Abstract
We consider recovery of low-rank matrices from noisy data by shrinkage
of singular values, in which a single, univariate nonlinearity is applied to
each of the empirical singular values. We adopt an asymptotic framework,
in which the matrix size is much larger than the rank of the signal matrix to be
recovered, and the signal-to-noise ratio of the low-rank piece stays constant.
For a variety of loss functions, including Mean Square Error (MSE – square
Frobenius norm), the nuclear norm loss and the operator norm loss, we show
that in this framework there is a well-defined asymptotic loss that we evaluate precisely in each case. In fact, each of the loss functions we study admits a
unique admissible shrinkage nonlinearity dominating all other nonlinearities.
We provide a general method for evaluating these optimal nonlinearities,
and demonstrate our framework by working out simple, explicit formulas
for the optimal nonlinearities in the Frobenius, nuclear and operator norm
cases. For example, for a square low-rank n-by-n matrix observed in white
noise with level σ, the optimal
simply shrinks each
p nonlinearity for MSE loss√
2
2
data singular value y to y − 4nσ (or to 0 if y < 2 nσ). This optimal
nonlinearity guarantees an asymptotic MSE of 2nrσ 2 , which compares favorably with optimally tuned hard thresholding and optimally tuned soft
thresholding, providing guarantees of 3nrσ 2 and 6nrσ 2 , respectively. Our
general method also allows one to evaluate optimal shrinkers numerically to
arbitrary precision. As an example, we compute optimal shrinkers for the
Schatten-p norm loss, for any p > 0.

Keywords. Matrix denoising — singular value shrinkage — optimal shrinkage
— spiked model — low-rank matrix estimation — nuclear norm loss — unique
admissible — Schatten norm loss.

∗
†

Department of Statistics, Stanford University
School of Computer Science and Engineering, Hebrew University of Jerusalem

1

Contents

2

Contents
1 Introduction
1.1 Optimal shrinkers vs. hard and soft thresholding . . . . . . . . . . .

3
7

2 Preliminaries
2.1 Natural problem scaling . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Asymptotic framework and problem statement . . . . . . . . . . . .
2.3 Our contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
8
9
10

3 The Asymptotic Picture

10

4 Optimal Shrinker for Frobenius Loss
4.1 Lower bound on asymptotic loss . . . . . . . . . . . . . . . . . . . .
4.2 Optimal shrinker matching the lower bound . . . . . . . . . . . . .

13
13
14

5 A Framework for Finding Optimal Shrinkers
5.1 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Simultaneous Block Diagonalization . . . . . . . . . . . . . . . . . .
5.3 Deterministic formula for the asymptotic loss . . . . . . . . . . . . .

15
16
17
19

6 Finding Optimal Shrinkers Analytically:
clear Losses
6.1 Frobenius norm loss . . . . . . . . . .
6.2 Operator norm loss . . . . . . . . . . .
6.3 Nuclear norm loss . . . . . . . . . . . .

21
22
23
23

Frobenius, Operator & Nu. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .

7 Finding Optimal Shrinkers Numerically: Schatten norm losses

24

8 Extensions
8.1 Unknown Noise level . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.2 General white noise . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26
26
28

9 Simulation

29

10 Conclusion

29

1 Introduction

3

1 Introduction
Suppose that we are interested in an m-by-n matrix X, which is thought to be
either exactly or approximately of low rank, but we only observe a single noisy
m-by-n matrix Y , obeying Y = X + σZ; The noise matrix Z has independent,
identically distributed entries with zero mean, unit variance, and a finite fourth
moment. We choose a loss function Lm,n (·, ·) and wish to recover the matrix X
with some bound on the risk ELm,n (X, X̂), where X̂ is our estimated value of X.
For example, when choosing the square Frobenius loss, or mean square error
(MSE)
2

ro
Lfm,n
(X, X̂) = X − X̂

F

=

X
i,j

|Xi,j − X̂i,j |2 ,

(1)

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

m
X

yi vi ṽi′

(2)

i=1

for the Singular Value Decomposition of the data matrix Y , where vi ∈ Rm and
ṽi ∈ Rn (for i = 1, . . . , m) are the left and right singular vectors of Y corresponding to the singular value yi . The TSVD estimator is
X̂r =

r
X

yi vi ṽi′ ,

i=1

where r = rank(X), assumed known, and y1 ≥ . . . ≥ ym . Being the best approximation of rank r to the data in the least squares sense [2], and therefore
the Maximum Likelihood estimator when Z has Gaussian entries, the TSVD is
arguably as ubiquitous in science and engineering as linear regression [3–8].
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 η : [0, ∞) → [0, ∞), also known as a shrinker, there is a corresponding
singular value shrinkage estimator X̂η given by
X̂η =

m
X

η(yi )vi ṽi′ .

(3)

i=1

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 [9–13]. Shrinkage makes sense for
singular values, too: presumably, the observed singular values y1 . . . ym are “inflated” by the noise, and applying a carefully chosen shrinkage function, one can
obtain a good estimate of the original signal X.

1 Introduction

4

Singular value shrinkage arises when the estimator X̂ for the signal X has to
be bi-orthogonally invariant under rotations of the data matrix. The most general
form of an invariant estimator is X̂ = V D̂Ṽ ′ , where the matrices V and Ṽ contain
the left and right singular vectors of the data, and where D̂ is a diagonal matrix
that depends on the data singular values. In other words, the most general invariant estimator is equivalent to a vector map Rm → Rm acting on the vector of
data singular values. This is a wide and complicated class of estimators; focusing
on univariate singular value shrinkage (3) allows for a simpler discussion1 .
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
[15,16] and Shabalin and Nobel [17]. Soft thresholding of singular values has been
considered in [18–20], and hard thresholding in [20, 21]. In fact, [16, 17] and, very
recently, [22, 23] 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 [20] 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 [15–17] and inspired by Johnstone’s Spiked Covariance
Model [24]. In this framework, the signal matrix dimensions m = mn and n
both to infinity, such that their ratio converges to an asymptotic aspect ratio:
mn /n → β, with 0 < β ≤ 1, while the column span of the signal matrix remains fixed. Building on a recent probabilistic analysis of this framework [25]
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.
1

However, it is interesting to remark that, at least in the Frobenius loss case, and possibly
in other cases as well, the asymptotically optimal univariate shrinker, presented in this paper,
offers the same performance asymptotically as the best possible invariant estimator of the form
X̂ = V D̂Ṽ ′ - see [14].

1 Introduction

5

Some optimal shrinkers
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 ||·||F ,||·||op and ||·||∗ , respectively. If the singular values of the matrix
X − X̂ are σ1 , . . . , σm , then these losses are given by
ro
Lfm,n
(X, X̂)

=

X − X̂

Lop
m,n (X, X̂) =

X − X̂

Lnuc
m,n (X, X̂) =

X − X̂

2

=

F

σi2

(4)

i=1

op

∗

m
X

= max {σ1 , . . . , σm }

=

m
X

σi .

(5)
(6)

i=1

Optimal shrinker for the Frobenius norm loss. As we will see, the optimal
nonlinearity for the Frobenius norm loss (4), in a natural noise scaling, is
( p
√
1
(y 2 − β − 1)2 − 4β y ≥ 1 + β
∗
y
(7)
η (y) =
√ .
0
y ≤ 1+ β
In the asymptotically square case β = 1 this reduces to
p
η(y) = (y 2 − 4)+ .

Optimal shrinker for Operator norm loss. The operator norm loss (5) for matrix estimation has mostly been studied in the context of covariance estimation
[26–28]. Let us define
( q
p
√
2−β −1+
√1
y
(y 2 − β − 1)2 − 4β y ≥ 1 + β
2
(8)
x(y) =
√ .
0
y ≤1+ β
As we will see, the optimal nonlinearity for operator loss is just
η ∗ (y) = x(y) .

(9)

Optimal Shrinkage for Nuclear norm loss. The Nuclear norm loss (6) has also
been proposed for matrix estimation. See [29, 30] 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

∗

η (y) =

(

1
(x4
x2 y

0

−β−

√

βxy)

√
x4 ≥ β + βxy
.
√
x4 < β + βxy

(10)

1 Introduction

6

where x = x(y) is given in (8).
√
Note that the formulas above are calibrated for the natural noise level σ = 1/ n;
see Section 8.1 below for usage in known noise level σ or unknown noise level.
In the code supplement for this paper [31] 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.

β = 0.2
fro
op
nu c

3.5

fro
op
nu c

3.5

3

3

2.5

2.5

η ∗ (y)

η ∗ (y)

β = 0.5

2

2

1.5

1.5

1

1

0.5

0.5

0

0
0

0.5

1

1.5

2

y

2.5

3

3.5

4

0

0.5

1

1.5

β = 0.7

3

3.5

4

2.5

3

3.5

4

fro
op
nu c

3.5

3

3

2.5

2.5

η ∗ (y)

η ∗ (y)

2.5

β =1

fro
op
nu c

3.5

2

y

2

2

1.5

1.5

1

1

0.5

0.5

0

0
0

0.5

1

1.5

2

y

2.5

3

3.5

4

0

0.5

1

1.5

2

y

Figure 1: Optimal shrinkers for Frobenius, Operator and Nuclear norm losses for different values of β. All shrinkers asymptote to the identity η(y) = y as y → ∞. Curves
jittered in the vertical axis to avoid overlap. This figure can be reproduced using the code
supplement [31].

2 Preliminaries

7

1.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 [9–13]. Recall that for y ≥ 0,
ηssof t (y) = max (0, y − s)
ηλhard (y) = y · 1y≥λ .
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 λ and s. In our asymtotic framework, fortunately, there is a decisive answer to the tuning question: in previous
work [20], we have restricted our attention to hard and soft thresholding under
the Frobenius loss (4). It was shown that there exist optimal values λ∗ (β) and
s∗ (β), 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
s
8β
p
λ∗ (β) =
2(β + 1) +
(β + 1) + β 2 + 14β + 1
p
s∗ (β) = 1 + β ,

where again β is the limiting aspect ratio, mn /n → β.
Consider, for example, the square matrix√case β = 1. Under the MSE loss,
the optimal hard threshold is then λ∗ = 4/ 3, and the optimal soft threshold
t
against our optimal
and ηssof
is s∗ = 2. Figure 2 shows the nonlinearities ηλhard
∗
∗
shrinkers (7), (9) and (10). In high SNR (y ≫ 1) the optimal shrinkers agree with
hard thresholding and neither performs any shrinkage, while soft thresholding
shrinks even strong signals. As shown in [20],√the worst-case asymptotic MSE
over a rank-r matrix observed in noise level 1/ n is 2r for our optimal shrinker
and 6r for the
(7), 3r for the optimally tuned hard thresholding nonlinearity ηλhard
∗
sof t
optimally tuned soft thresholding nonlinearity η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 [20]. We conclude that optimal shrinkage,
developed in this paper, offers significant performance improvement over hard
and soft thresholding - even when they are optimally tuned.

2 Preliminaries
Column vectors are denoted by boldface lowercase letters, such as v, their transpose is v′ and their i-th coordinate is vi . The Euclidean inner product and norm
on vectors are denoted by hu , vi and ||u||2 , respectively. Matrices are denotes by

2 Preliminaries

8

β =1
9

β =1

op
h ard
soft

9

8

7

7

7

6

6

6

η ∗ (y)

8

5

5

5

4

4

4

3

3

3

2

2

2

1

1

1

0

0
0

2

4

y

6

8

10

nu c
h ard
soft

9

8

η ∗ (y)

η ∗ (y)

β =1

fro
h ard
soft

0
0

2

4

y

6

8

10

0

2

4

y

6

8

Figure 2: Optimal shrinkers for Frobenius (left), Operator (center) and Nuclear norm
(right) losses for square matrices (β = 1) plotted with optimally tuned hard and soft
thresholding nonlinearities. This figure can be reproduced using the code supplement
[31].

uppercase letters, such as X, their transpose is X ′ and theirP
i, j-th entry is Xi,j .
Mm×n denotes the space of real m-by-n matrices, hX , Y i = i,j Xi,j Yi,j denotes
the Hilbert-Schmidt inner product and ||X||F denotes the corresponding Frobenius norm on Mm×n . ||X||∗ and ||X||op denote the nuclear norm (sum of singular
values) and operator norm (maximal singular value) of X, respectively. For simplicity we only consider m ≤ n. We denote matrix denoisers, or estimators, by
a.s.
a.s.
X̂ : Mm×n → Mm×n . The symbols → and = denote almost sure convergence and
equality of a.s. limits, respectively. We use “fat” SVD of a matrix X ∈ Mm×n with
m ≤ n, that is, when writing X = UD Ũ we mean that U ∈ Mm×m , D ∈ Mm×n ,
and Ũ ∈ Mn×n . Symbols without tilde such as u are associated with left singular vectors, while symbols with tilde such as ũ are associated with right singular
vectors. By diag(x1 , . . . , xm ) we mean the m-by-n matrix whose main diagonal is
x1 , . . . , xm , with n implicit in the notation and inferred from context.

2.1 Natural problem scaling
In
√ the general model Y = X + σZ, the noise level in the singular values of Y is
nσ. 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/ n.
In this convention, shrinkage rules stay the same for every value of n, and we
conveniently abuse notation by writing X̂η as in (3) for any X̂η : Mm×n → Mm×n ,
keeping m and n implicit. To apply any denoiser X̂ below to data from the general model Y = X + σZ, use the denoiser
√
√
X̂η(n,σ) (Y ) = nσ · X̂η (Y / nσ) .
(11)
Throughout the
√ text, we use X̂η to denote singular value shrinker calibrated
for noise level 1/ n. In Section 8.1 below we provide a recipe for applying any
√
denoiser X̂η calibrated for noise level σ = 1/ n for data in the presence of unknown noise level.

10

2 Preliminaries

9

2.2 Asymptotic framework and problem statement
In this paper, we consider a sequence of increasingly larger denoising problems
√
(12)
Y n = Xn + Z n / n
with Xn , Zn ∈ Mmn ,n , satisfying the following assumptions:
1. Invariant white noise: The entries of Zn 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 Zn follows the same distribution as AZn B,
for every orthogonal A ∈ Mmn ,mn and B ∈ Mn,n . This is the case, for example, when the entries of Zn are Gaussian. In Section 8.2 we revisit this
restriction and discuss general (not necessarily invariant) white noise.
2. Fixed signal column span (x1 , . . . , xr ): Let the rank r > 0 be fixed and choose a
vector x ∈ Rr with coordinates x = (x1 , . . . , xr ) such that x1 > . . . > xr > 0.
Assume that for all n,
Xn = Un diag(x1 , . . . , xr , 0, . . . , 0) Ũn′

(13)

is an arbitrary singular value decomposition of Xn , where Un ∈ Mmn ,mn and
Ũn ∈ Mn,n .
3. Asymptotic aspect ratio β: The sequence mn is such that mn /n → β. To simplify our formulas, we assume that 0 < β ≤ 1.
Note that while the signal rank r and nonzero signal singular values x1 , . . . , xr
are shared by all matrices Xn , the signal left and right singular vectors Un and Vn
are unknown and arbitrary. We also remark that the assumption, whereby the
signal singular values are non-degenerate (xi > xi+1 , 1 ≤ i < r), is not necessary
for our results to hold, yet it simplifies the analysis considerably.
Definition 1. Asymptotic Loss. Let L = {Lm,n | (m, n) ∈ N×N} be a family of losses,
where each Lm,n : Mm×n × Mm×n → [0, ∞) is a loss function obeying Lm,n (X, X) = 0.
Let η : [0, ∞) → [0, ∞) be a nonlinearity and consider X̂η , √
the singular value shrinkage
denoiser (3) calibrated, as discussed above, for noise level 1/ n. Let mn be an increasing
sequence such that limn→∞ mn /n = β, implicit in our notation. Define the asymptotic
loss of the shrinker η (with respect to L) at the signal x = (x1 , . . . , xr ) by


a.s.
L∞ (η|x) = lim Lmn ,n Xn , X̂η (Xn + √1n Zn )
n→∞

when the limit exists.

Our results imply that the asymptotic loss L∞ exists and is well-defined, as a
function of the signal singular values x, for a large class of nonlinearities.

3 The Asymptotic Picture

10

Definition 2. Optimal Shrinker. Let L be a loss family. If a shrinker η ∗ has an asymptotic loss that satisfies
L∞ (η ∗ |x) ≤ L∞ (η|x)

for any other shrinker η in a certain class of shrinkers, any r ≥ 1 and any x ∈ Rr , then
we say that η ∗ is unique asymptotically admissible (of simply “optimal”) for the loss
sequence L and that class of shrinkers.

2.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 [17] 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.
1. We rigorously establish the existence of an optimal shrinker for a variety of
loss families, including the popular Frobenius, operator and nuclear norm
losses.
2. 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 σ is known or unknown.
3. 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 closedform 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 [14].

3 The Asymptotic Picture
In the “null
√ case” Xn ≡ 0, the empirical distribution of the singular values of
Yn = Zn / n famously converges as n → ∞ to the generalized quarter-circle
distribution [32], whose density is
p
4β − (x2 − 1 − β)2
(14)
1[1−√β,1+√β] (x) .
f (x) =
πβx

3 The Asymptotic Picture

11

This distribution is compactly supported on [β− , β+ ], with
p
β± = 1 ± β .

√
a.s.
Moreover, in this null case we have yn,1 → 1 + β, see [33, 34]. We say that the
singular values of Yn form a (generalized) quarter circle bulk and call β+ the bulk
edge.
Expanding seminal results of [35,36] and many other authors, Benaych-Georges
and Nadakuditi [25] 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.
3

14

12

2.5

10
2
8
1.5
6
1
4
0.5

0

2

0

20

40

60

80

0

100

0

0.5

1

1.5

2

2.5

3

Figure √
3: Singular values of a data matrix Y ∈ M100,100 drawn from the model Y =
X + Z/ 100, with r = 2 and x = (2.5, 1.7). Left: singular values in decreasing order.
Right: Histogram of the singular values (note the bulk edge close to 2). This figure can
be reproduced using the code supplement [31].
For x ≥ β 1/4 , define
y(x) =

s

1
x+
x



β
x+
,
x

c(x) =

s

x4 − β
x4 + βx2

c̃(x) =

r

x4 − β
.
x4 + x2

and

(15)
(16)
(17)

It turns out that y(x) from Eq. (15) is the asymptotic location of a data singular value corresponding to a signal singular value x, provided x ≥ β 1/4 . (Note
that the function x(y) from (8) is the inverse of y(x) when x ≥ β 1/4 , and that
y(β 1/4) = β+ .) Similarly, c(x) from Eq. (16) (resp. c̃(x) from Eq. (17)) is the cosine of the asymptotic angle between the signal left (resp. right) singular vector

3 The Asymptotic Picture

12

and the corresponding data left (resp. right) singular vector, provided that the
corresponding signal singular value x satisfies x ≥ β 1/4 .
Additional notation is required to state these facts formally. We rewrite the
sequence of signal matrices in our asymptotic framework (13) as
Xn =

r
X

xi un,i ũ′n,i ,

(18)

i=1

so that un,i ∈ Rmn (resp. ũn,i ∈ Rn ) is the left (resp. right) singular vector corresponding to the singular value xi , namely, i-th column of Un (resp. Ũn ) in (13).
Similarly, let Yn be a corresponding sequence of observed matrices in our framework, and write
Yn =

mn
X

′
yn,i vn,i ṽn,i

(19)

i=1

so that vn,i ∈ Rm (resp. ṽn,i ∈ Rn ) is the left (resp. right) singular vector corresponding to the singular value yn,i .
In our notation, Lemma 1 and Lemma 2 follow from Theorem 2.9 and Theorem
2.10 of [25]:
Lemma 1. Asymptotic location of the top r data singular values. For 1 ≤ i ≤ r,
(
y(xi ) xi ≥ β 1/4
a.s.
.
(20)
lim yn,i =
n→∞
β+
xi < β 1/4
Lemma 2. Asymptotic angle between signal and data singular vectors. Let 1 ≤
i 6= j ≤ r and assume that xi ≥ β 1/4 is non-degenerate, namely, the value xi appears
only once in x. Then
(
c(xi ) i = j
a.s.
,
(21)
lim hun,i , vn,j i =
n→∞
0
i 6= j
and
(
c̃(xi ) i = j
lim hũn,i , ṽn,j i =
.
n→∞
0
i=
6 j
a.s.

(22)

If however xi < β 1/4 , then we have
a.s.

a.s.

lim hun,i , vn,j i = lim hũn,i , ṽn,j i = 0 .

n→∞

n→∞

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

Lemma 3. Let i > r be fixed. Then yn,i → β+ .

4 Optimal Shrinker for Frobenius Loss

13

4 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 [17].
Using Definition 1, let L = {Lm,n } be the Frobenius loss family, namely Lm,n is
given by (4).

4.1 Lower bound on asymptotic loss
Directly expanding the Frobenius matrix norm, we obtain:
Lemma 4. Frobenius loss of singular value shrinkage. For any shrinker η : [0, ∞) →
[0, ∞), we have
r
X
2
 2

√
Xn − X̂η (Xn + Zn / n)
=
xi + (η(yn,i ))2
(23)
F

i=1

− 2

r
X

i,j=1
m
n
X

xi η(yn,i )hun,i , vn,j ihũn,i , ṽn,j i

(η(yn,i ))2

+

(24)
(25)

i=r+1

This implies a lower bound on Frobenius loss of any singular value shrinker:
Corollary 1. For any shrinker η : [0, ∞) → [0, ∞), we have
r
X
2
 2

√
Xn − X̂η (Xn + Zn / n)
≥
xi + (η(yn,i ))2
F

i=1

− 2

r
X

i,j=1

xi η(yn,i )hun,i , vn,j ihũn,i , ṽn,j i .

As n → ∞, this lower bound on the Frobenius loss is governed by three quantities: the asymptotic location of data singular value yn,i , the asymptotic angle
between the left signal singular vectors and left data singular vector hun,i , vn,i i,
and asymptotic angle between the right signal singular vectors and right data
singular vector hũn,i , ṽn,i i (see also [20]).
Combining Corollary 1, Lemma 1 and Lemma 2 we obtain a lower bound for
the asymptotic Frobenius loss (see [17]):
Corollary 2. For any continuous shrinker η : [0, ∞) → [0, ∞), we have
r
X
2
√
a.s.
L∞ (η|x) = lim Xn − X̂η (Xn + Zn / n)
≥
L2,2 (η|xi )
n→∞

F

i=1

where

and η = η(y).

L2,2 (η|x) = x2 + η 2 − 2xηc(x)c̃(x) .

The notation L2,2 in (26) will be made apparent below, see (48).

(26)

4 Optimal Shrinker for Frobenius Loss

14

4.2 Optimal shrinker matching the lower bound
By differentiating the asymptotic lower bound w.r.t η, we find that L∞ (η|x) ≥
P
r
∗
∗
i=1 L2,2 (η |xi ), where η (y(x)) = xc(x)c̃(x). Expanding c(x) and c̃(x) from Eqs.
(16) and (17), we find that η ∗ (y) is given by (7).
The singular value shrinker η ∗ , for which X̂η∗ minimizes the asymptotic lower
bound, thus becomes a natural candidate for the optimal shrinker for Frobenius
loss. Indeed, by definition, for X̂η∗ 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 η 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 η
for which η(y) = 0 for any y ≤ β+ . 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.
Definition 3. Assume that a continuous shrinker η : [0, ∞) → [0, ∞) satisfies η(y) = 0
whenever y ≤ β+ + ε for some fixed ε > 0. We say that η is a Conservative shrinker.
By Lemma 1 and Lemma 2, it is clear that conservative shrinkers set to zero
all data singular values {yi } which originate from pure noise (xi = 0), as well as
all data singular values {yi } which are “engulfed” in the noise bulk, rendering
their corresponding singular vectors useless (xi < β 1/4 ). Conservative shrinkers
are so called since they leave a (possibly infinitesimally small) safety margin ε.
They enjoy the following key property:
Lemma 5. Let η : [0, ∞) → [0, ∞) be a conservative shrinker. Then
mn
X

a.s.

(η(yn,i))2 → 0 .

i=r+1

a.s.

Proof. By Lemma 3 we have yn,r+1 → β+ . Let N be the (random) index such that
yn,r+1 < β+ +ε for all n > N. Then for all n > N and all i > r we have yn,i < β+ +ε,
hence η(yn,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 η(y) = 0 for y ≤ β+ , leaving no margin above the bulk edge β+ . 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 [14]; 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
(x1 , . . . , xr ).

5 A Framework for Finding Optimal Shrinkers

15

5 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.
Definition 4. Orthogonally invariant loss. A loss Lm,n (·, ·) is orthogonally invariant if for all m, n we have Lm,n (A, B) = Lm,n (UAV, UBV ), for any orthogonal U ∈ Om
and V ∈ On .
P
Definition 5. Decomposable loss family. Let A, B ∈ Mm×n and let m = ki=1 mi
P
and n = ki=1 ni . Assume that there are matrices Ai , Bi ∈ Mmi ,ni , 1 ≤ i ≤ k, such that
A = ⊕i Ai

B = ⊕i Bi

in the sense that A and B are block-diagonal with blocks {Ai } and {Bi }, respectively.
A loss family L = {Lm,n } is sum-decomposable if, for all m, n and A, B with block
diagonal structure as above,
X
Lm,n (A, B) =
Lmi ,ni (Ai , Bi ) .
i

Similarly, it is max-decomposable if

Lm,n (A, B) = max Lmi ,ni (Ai , Bi ) .
i

Examples. As primary examples, we consider loss families defined in Section 1:
The Frobenius norm loss Lf ro , the operator norm loss Lop and the nuclear norm
loss Lnuc . It is easy to check that (i) each of these losses are orthogonally invariant,
and (ii) the families Lf ro and Lnuc are sum-decomposable, while the family Lop is
max-decomposable.
Our framework for finding optimal shrinkers can now be stated as follows.
Theorem 1. Characterization of the optimal singular value shrinker. Let


x 0
(27)
A(x) =
0 0


c(x) c̃(x) c(x) s̃(x)
,
(28)
B(η, x) = η
c̃(x) s(x) s(x) s̃(x)
p
1 − c2 (x) and
where c(x)
and
c̃(x)
are
given
by
Eqs.
(16)
and
(17),
and
where
s(x)
=
p
s̃(x) = 1 − c̃2 (x). Assume that L = {Lm,n } is a sum- or max- decomposable family of
orthogonally invariant losses. Define
F (η, x) = L2,2 (A(x), B(η, x))

(29)

5 A Framework for Finding Optimal Shrinkers

16

and suppose that for any x ≥ β 1/4 there exists a unique minimizer
η ∗∗ (x) = argminη≥0 F (η, x) ,

(30)

such that η ∗∗ is a conservative shrinker on [β 1/4 , ∞). Further suppose that there exists a
point x0 ≥ β 1/4 such that
F (η ∗∗ (x), x) ≥ L1,1 (x, 0)

β 1/4 ≤ x ≤ x0 ,

with
F (η ∗∗ (x0 ), x0 ) = L1,1 (x0 , 0) .

(31)

(
η ∗∗ (x(y)) y(x0 ) ≤ y
,
η ∗ (y) =
0
0 ≤ y < y(x0 )

(32)

Define the shrinker

where x(y) is defined in Eq. (8). Then for any conservative shrinker η, the asymptotic
losses L∞ (η ∗ |·) and L∞ (η|·) exist, and
L∞ (η ∗ |x) ≤ L∞ (η|x)
for all r ≥ 1 and all x ∈ Rr .

5.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 η is a conservative shrinker, then the asymptotic loss
L∞ (η|x) at x = (x1 , . . . , xr ) can be written as a sum (resp. a maximum) over
r terms. These terms have identical functional form. When xi ≥ β 1/4 , these
terms have the form L2,2 (A(xi ), B(η, xi )), and when 0 ≤ xi < β 1/4 , these
terms have the form L1,1 (xi , 0) + L1,1(0, η) (resp. max{L1,1 (xi , 0) , L1,1 (0, η)}).
As a result, one finds that the zero shrinker η ≡ 0 is necessarily optimal for
0 ≤ x ≤ β 1/4 . For x ≥ β 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 η ∗∗ of
(30).
• Second, the asymptotic loss curve L∞ (η ∗∗ |x) necessarily crosses the asymptotic loss curve of the zero shrinker L∞ (η ≡ 0|x) at a point we will denote
by x0 , with x0 ≥ β 1/4 .
• Finally, by concatenating the zero shrinker and the shrinker η ∗∗ precisely
at the point x0 where their asymptotic losses cross, one obtains a shrinker
which is continuous (x0 > β 1/4 ) or possibly discontinuous (x0 = β 1/4 ). However, this shrinker always has a well-defined asymptotic loss. This loss
dominates the asymptotic loss of any conservative shrinker.

5 A Framework for Finding Optimal Shrinkers

17

For some loss families L = {Lm,n }, it is possible to find an explicit formula for
the optimal shrinker using the following steps:
1. Write down an explicit expression for the function F (η, x) from (29).
2. Explicitly solve for the minimizer η ∗∗ (x) from (30).
3. Write down an explicit expression for the minimum F (η ∗∗ (x), x).
4. Solve (31) for the crossing point x0 .
5. Compose η ∗∗ (x) with the transformation x(y) from (8) to obtain an explicit
form of the optimal shrinker η ∗ (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.

5.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 η shrinks to zero all but the first singular values of√Y , namely,
η(yi ) = 0, i > 1. Let X ∈ Mm×n be a signal matrix and let Y = X + Z/ n ∈ Mm×n
be a corresponding data matrix. Denote their SVD by
X = U · diag(x1 , 0, . . . , 0) · Ũ

Y = V · diag(y1 , . . . , ym ) · Ṽ .

Write 0m×n for the m-by-n matrix whose entries are all zeros. The basis pairs U, Ũ
and V, Ṽ diagonalize X and Y , respectively. Indeed, since η(yi ) = 0 for i ≥ 2, we
have
X = x1 u1 ũ′1

X̂η (Y ) = η(y1)v1 ṽ1′ .

Combining the basis pairs U, Ũ and V, Ṽ , we are lead to the following “common” basis pair, which we will denote by W, W̃ . Let w1 , . . . , wm denote the orthonormal basis constructed by applying the Gram-Schmidt process to the sequence u1 , v1 , v2 , . . . , vm−1 , where ui is the i-th column of U, namely the i-th left
singular vector of X, and similarly, vi is the i-th column of V . Denote by W the
matrix whose columns are w1 , . . . , wm . Repeating this construction for Ũ and Ṽ ,
let w̃1 , . . . , w̃n denote the orthonormal basis constructed by applying the GramSchmidt process to the sequence ũ1 , ṽ1 , ṽ2 , . . . , ṽm−1 , where ũi is the i-th column
of U, namely the i-th right singular vector of X, and similarly, ṽi is the i-th column
of Ṽ . Denote by W̃ the matrix whose columns are w̃1 , . . . , w̃m .

5 A Framework for Finding Optimal Shrinkers

18

Specifically, if v1 and u1 are not colinear, we let w1 = u1 and let w2 be the first
Gram-Schmidt step for the sequence u1 , v1 , vp
2 , . . . , vm−1 , namely, we let w2 =
−1
s1 (v1 − c1 u1 ), where c1 = hu1 , v1 i and s1 = 1 − c21 . If it happens that v1 and
u1 are colinear, we let w2 be any vector orthogonal to u1 , for example um . The
rest of the Gram-Schmidt process proceeds to add m − 2 additional unit vectors
orthogonal to both w1 and w2 . Observe that Span {u1 , v1 } ⊂ Span {w1 , w2 }, so
that W ′ u1 = (1, 0, . . . , 0)′ and W ′ v1 = (c1 , s1 , 0, . . . , 0)′ . Repeating the same construction for the right singular vector basis Ũ and Ṽ , we obtain the basis W̃ with
the property that W̃ ′ ũ1 = (1, 0, . . . , 0)′ and W̃ ′ ṽ1 = (c̃1 , s̃1 , 0, . . . , 0)′ .
Writing now X and X̂η in our new basis pair W, W̃ we get


x1 0
′
′
′
′
W X W̃ = x1 (W u1 )(W̃ ũ1 ) =
⊕ 0m−2×n−2
(33)
0 0


c1 c̃1 c1 s̃1
′
′
′
′
W X̂η (Y ) W̃ = η(y1) (W v1 )(W̃ ṽ1 ) = η(y1 )
⊕ 0m−2×n−2 . (34)
c̃1 s1 s1 s̃1
It is convenient to rewrite this as
W ′ X W̃ = A(x1 ) ⊕ 0m−2×n−2
W ′ X̂η (Y ) W̃ = B(η(y1 ), c1 , s1 , c̃1 , s̃1 ) ⊕ 0m−2×n−2

(35)
(36)

where


x 0
A(x) =
0 0


cc̃ cs̃
.
B(η, c, s, c̃, s̃) = η
c̃s ss̃

(37)
(38)

Thus, if L = {Lm,n } is a sum- or max-decomposable family of orthogonally
invariant functions, we have
Lm,n (X , X̂η (Y )) = Lm,n (W ′ X W̃ , W ′ X̂η (Y )W̃ )


= Lm,n A(x1 ) ⊕ 0m−2×n−2 , B(η(y1 ), c1 , s1 , c̃1 , s̃1 ) ⊕ 0m−2×n−2


= L2,2 A(x1 ) , B(η(y1), c1 , s1 , c̃1 , s̃1 ) .
We have proved:

P
′
Lemma 6. Let X = x1 u1 ũ′1 ∈ Mm×n be rank-1 and assume that Y = m
i=1 yi vi ṽi and
η are such that η(yi ) = 0, i > 1, where yi is the i-th largest singular value of Y . Let
L = {Lm,n } be a sum- or max-decomposable orthogonally invariant loss family. Then


Lm,n (X , X̂η (Y )) = L2,2 A(x1 ) , B(η(y1 ), c1 , s1 , c̃1 , s̃1 ) ,

where

c1 = hu1 , v1 i
c̃1 = hũ1 , ṽ1 i

q
s1 = 1 − c21
q
s̃1 = 1 − c̃21 .

5 A Framework for Finding Optimal Shrinkers

19

A similar argument gives a similar statement for rank-r matrix X with nondegenerate singular values:
P
Lemma 7. A decomposition for the loss. Let X =P ri=1 xi ui ũ′i ∈ Mm×n be rankm
′
r with x1 > . . . > xr > 0, and assume that Y =
i=1 yi vi ṽi and η are such that
η(yi ) = 0, i > r, where yi is the i-th largest singular value of Y . Let L = {Lm,n } be a
sum- or max-decomposable family of orthogonally invariant functions. Then
Lm,n (X , X̂η (Y )) =

r
X
i=1



L2,2 A(xi ) , B(η(yi ), ci , si , c̃i , s̃i ) ,

if L is sum-decomposable, or
r





Lm,n (X , X̂η (Y )) = max L2,2 A(xi ) , B(η(yi ), ci , si , c̃i , s̃i ) ,
i=1

if L is max-decomposable. Here,
ci = hui , vi i

si =

c̃i = hũi , ṽi i

s̃i =

for i = 1, . . . , r.

q

q

1 − c2i
1 − c̃2i ,

See [14] for the complete proof. In short, we apply the Gram-Schmidt process
to the vector sequences u1 , . . . , ur , v1 , . . . , vm−r and ũ1 , . . . , ũr , ṽ1 , . . . , ṽm−r to obtain the orthogonal matrices W and W̃ , whose columns constitute a basis pair
similar to the one above. In this basis pair,
i
h
(39)
W ′ X W = diag(x1 , . . . , xr ) ⊕ Ir ⊕ Im−2r
W ′ X̂η (Y ) W = Rn ⊕ In−2r

(40)

where Rn is a sequence of 2r-by-2r matrices such that

a.s.
Rn → ⊕ri=1 B η(yi ), c1 , s1 , c̃1 , s̃1 .

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

5.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 Lm,n (X , X̂η (Y )) decomposes to “atomic” units of the form


L2,2 A(xi ) , B(η(yi ), ci , si , c̃i , s̃i ) .

√
Let us now return to the sequence model Yn = Xn + Zn / n and find the
limiting value of these “atomic” units as n → ∞. This will lead to a simple
formula for the asymptotic loss L∞ (η|x).

5 A Framework for Finding Optimal Shrinkers

20

Each of these “atomic” units only depend on yi , the i-th data singular value,
and on ci (resp. c̃i ), the angle between the i-th left (resp. right) signal and data
singular vectors. In the special case of Frobenius norm loss, we have already
encountered this phenomenon (Lemma 4), where we have seen that these quantities converge to deterministic functions that depend on xi , the i-th signal singular
value alone.
√
For the sequence Yn = Xn + Zn / n, recall our notation yn,i , un,i , ũn,i , vn,i , ṽn,i
from (18) and (19), and define
q
cn,i = hun,i , vn,i i
si = 1 − c2n,i
q
c̃n,i = hũn,i , ṽn,i i
s̃i = 1 − c̃2n,i ,

for i = 1, . . . , r. Combining Lemma 7, Lemma 1 and Lemma 2 we obtain:
√
Lemma 8. Let Yn = Xn + Zn / n be a matrix sequence in our asymptotic framework
with signal singular values x = (x1 , . . . , xr ). Assume that η is continuous at y(xi ) for
some fixed 1 ≤ i ≤ r. If β 1/4 ≤ xi then


a.s.
lim L2,2 (A(xi ) , B(η(yn,i ), cn,i , sn,i , c̃n,i , s̃n,i ) = L2,2 A(xi ) , B(η(y(xi )), xi ) ,
n→∞

where B(η, x) is given by (28), while if 0 ≤ xi < β 1/4 then
a.s.

lim L2,2 (A(xi ) , B(η(yn,i ), cn,i , sn,i , c̃n,i , s̃n,i ) =

n→∞

=



L2,2 A(xi ) , diag(0, η(y(xi)))
L1,1 (xi , 0) + L1,1 (0, η(y(xi ))) .

As a result, we now obtain the asymptotic loss L∞ as a deterministic function
of the nonzero signal singular values x1 , . . . , xr . Observe that by Lemma 3, if η
is a conservative shrinker, then eventually η(yn,i ) = 0 for all i > r. Therefore
the assumption η(yi ) = 0 for i > r, required for Lemma 7, is satisfied eventually.
Combining Lemma 7 and Lemma 8, we obtain
Lemma 9. A formula for the asymptotic loss of a conservative shrinker. Assume
that L = {Lm,n } is a sum- or max- decomposable family of orthogonally invariant losses.
Extend the definition of B(η, x) from (28) by setting B(η, x) = diag(0, η) for 0 ≤ x <
β 1/4 . If η : [0, ∞) → [0, ∞) is a conservative shrinker, then
r
X



L2,2 A(xi ), B(η(y(xi )), xi )

(41)



r
L∞ (η|x) = max L2,2 A(xi ), B(η(y(xi )), xi )

(42)

L∞ (η|x) =

i=1

if L is sum-decomposable, or
i=1

if L is max-decomposable.

6 Finding Optimal Shrinkers Analytically: Frobenius, Operator & Nuclear
Losses

21

The final step toward the proof of Theorem 1 involves the case when the
shrinker η is given as a special concatenation of two conservative shrinkers. Even
if the two parts of η 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.
Definition 6. Assume that there exist a point 0 < x∗ and two shrinkers, η1 : [0, x∗ ) →
[0, ∞) and η2 : [x∗ , ∞) → [0, ∞), such that




∗
∗
∗
∗
∗
∗
L2,2 A(x ), B(η1 (y(x )), x ) = L2,2 A(x ), B(η2 (y(x )), x ) .
We say that the asymptotic loss functions of η1 and η2 cross at x∗ .

Lemma 10. A formula for the asymptotic loss of a concatenation of two conservative shrinkers. Assume that L = {Lm,n } is a sum- or max- decomposable family of orthogonally invariant losses. Extend the definition of B(η, x) from (28) by setting B(η, x) = diag(0, η) for 0 ≤ x < β 1/4 . Assume that there exist two shrinkers,
η1 : [0, x∗ ) → [0, ∞) and η2 : [x∗ , ∞) → [0, ∞), whose asymptotic loss functions cross
at some point 0 < x∗ . Define
(
η1 (x) 0 ≤ x ≤ x∗
.
η(x) =
η2 (x) x∗ < x
Then L∞ (η|·) exists and is given by (41) if L is sum-decomposable, or (42) if L is maxdecomposable.
Proof of Theorem 1. Consider the shrinker η1 ≡ 0. By Lemma 8, η1 dominates
any other conservative shrinker when 0 ≤ x < β 1/4 . By assumption, there exists a point β 1/4 ≤ x0 such that η1 also dominates any conservative shrinker on
[β 1/4 , x0 ), and such that η ∗∗ dominates any other conservative shrinker on [x0 , ∞).
Finally, by assumption, the asymptotic loss functions of η1 and η ∗∗ cross at x0 . By
Lemma 10, the concatenated shrinker η ∗ dominates any conservative shrinker on
[0, ∞).

6 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 [14]):

6 Finding Optimal Shrinkers Analytically: Frobenius, Operator & Nuclear
Losses

22

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


p
1
2
λ± (M) = 2 trace(M) ± trace(M) − 4det(M) .
(43)
Proof. These are the roots of the characteristic polynomial of M.

Lemma 12. Let ∆ be a 2-by-2 matrix with singular values σ+ > σ− > 0. Define
t = trace(∆∆′ ) = ||∆||2F , d = det(∆) and r 2 = t2 − 4d2 . Assume that ∆ depends
on a parameter η and let σ̇± , ṫ and d˙ denote the derivative of these quantities w.r.t the
parameter η. Then
˙ ṫ − 2d)
˙ .
r(σ̇+ + σ̇− )(σ̇+ − σ̇− ) = 2(ṫ + 2d)(
2
Proof. By Lemma 11 we have 2σ±
= t ± r and therefore

√

ṫ ± ṙ
2σ̇± = √
.
2 t±r

Differentiating and expanding σ̇+ ± σ̇− we obtain the relation
(σ̇+ + σ̇− ) =

˙ ṫ − 2d)
˙
(8d/r)(ṫ + 2d)(
(t2 − r 2 )(σ̇+ − σ̇− )

(44)

and the result follows.
√
1 − c2 and s̃ = 1 − c̃2 . Define


ηcc̃ − x ηcs̃
∆ = ∆(η, c, c̃, s, s̃) =
ηc̃s
ηss̃ .

Lemma 13. Let η, c, c̃ ≥ 0 and set s =

√

Then
||∆||2F = η 2 + x2 − 2xηcc̃
det(∆) = −xηss̃ ,
and the singular values σ+ > σ− of ∆ are given by
r
q
σ± =

√1
2

||∆||2F ±

||∆||4F − 4det(∆)2 .

(45)
(46)

(47)

Proof. Apply Lemma 11 to the matrix ∆.

6.1 Frobenius norm loss
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

 
 
x 0
cc̃ cs̃
f ro
(48)
= ||∆||2F = η 2 + x2 − 2xηcc̃ .
,
L2,2 η
0 0
c̃s ss̃

6 Finding Optimal Shrinkers Analytically: Frobenius, Operator & Nuclear
Losses

23

To find the optimal shrinker, we solve ∂ ||∆||2F /∂η = 0 for η and use the fact that
c2 c̃2 + c2 s̃2 + s2 c̃2 + s2 s̃2 = 1. We find that the minimizer of ||∆||2F is η ∗∗ (x) = x cc̃.
Defining η ∗∗ (x) = x c(x)c̃(x) for x ≥ β 1/4 , we find that the asymptotic loss of
η ∗∗ (x) and of η ≡ 0 cross at x0 = β 1/4 . Simplifying (32), where x(y) is given by (8),
we find that η ∗ (y) is given by (7). By Theorem 1, this is the optimal shrinker. (As
mentioned in Section 4, the optimal shrinker is not itself strictly conservative, yet
can be shown to have the asymptotic loss predicted by our framework.)

6.2 Operator norm loss
By (47),
Lop
2,2


 
 
x 0
cc̃ cs̃
= ||∆||op = σ+ .
,
η
0 0
c̃s ss̃

(49)

To find the optimal shrinker, we solve ∂ ||∆||op /∂η = 0 for η on β 1/4 ≤ x. We
find that the minimizer of ||∆||op is at η ∗∗ (x) = x. By (47), the asymptotic loss of
p
η ∗∗ is given by x 1 − c(x)c̃(x) + |c(x) − c̃(x)|, so that the asymptotic loss of η ∗∗ (x)
and of η ≡ 0 cross at x0 = β 1/4 . Simplifying (32), we recover the shrinker η ∗ (y)
mentioned above√in (9). Observe that although this shrinker is discontinuous
at y(β 1/4 ) = 1 + β, its asymptotic loss L∞ (η ∗ |·) exists by Lemma 10, and, by
Theorem 1, dominates the asymptotic loss any conservative shrinker. We note
that as in the Frobenius case, this optimal shrinker is not strictly conservative,
as it only satisfies η(y) = 0 for y ≤ β+ . Here too a delicate argument beyond
our present scope shows that the asymptotic loss exists and matches the formula
predicted by our framework [14].
Remark. The optimal shrinker for operator norm loss η ∗ (y) = x(y) simply shrinks
the data singular value back to the ”original” location of its corresponding signal
singular value.

6.3 Nuclear norm loss
Again by (47)
Lnuc
2,2

 
 

cc̃ cs̃
x 0
η
,
= ||∆||∗ = σ+ + σ− .
c̃s ss̃
0 0

(50)

To find the optimal shrinker, assume first that x is such that c(x)c̃(x) ≥ s(x)s̃(x).
By Lemma 12, with
t = η 2 + x2 − 2xηcc̃

and

d = −xηss̃ ,

we find that only zero of ∂(σ+ + σ− )/∂η occurs when ∂t/∂η − ∂d/∂η = 0, namely
at


∗∗
η (x) = x c(x)c̃(x) − s(x)s̃(x) .

7 Finding Optimal Shrinkers Numerically: Schatten norm losses

24

Direct calculation using (47), (46) and (45) shows
 that the square of the
 asymptotic

loss of η ∗∗ is simply x2 + (η ∗∗ (x))2 − 2xη ∗∗ (x) c(x)c̃(x) − s(x)s̃(x) , so that the
asymptotic loss of η ∗∗ (x) and of η ≡ 0 cross at the unique x0 satisfying c(x0 )c̃(x0 ) =
s(x
c(x) from (16), c̃ = c̃(x) from (17) and also s = s(x) =
p 0 )s̃(x0 ). Substituting c =p
2
1 − c(x) and s̃ = s̃(x) = 1 − c̃(x)2 , we find
√ 
 4
β
x −β
∗
,
+
η (y) =
2
xy
x +
recovering the optimal shrinker (9). Inspection of (9) reveals that this optimal
shrinker is in fact a conservative shrinker.

7 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 numerically2 , for any value p > 0.
Let 0 < p ≤ ∞. Recall that the Schatten-p norm (resp. quasi-norm if 0 < p < 1)
of a matrix is the ℓp norm (resp. quasi-norm) of its singular values vector. If the
singular values of the m-by-n matrix X − X̂ are σ1 , . . . , σm , define the Schatten-p
loss
v
u m
uX p
p
p
LSm,n
(X, X̂) = t
σi
0<p<∞
i=1

p
(X, X̂) = max {σ1 , . . . , σm }
LSm,n

p = ∞,

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 [29, 30] and references therein. (The case 0 < p ≤ 1 is of
special interest in matrix completion problems due to its low-rank inducing bero
havior.) So far in this paper we have carefully studied three special cases: Lfm,n
≡
S2
nuc
S1
op
S∞
Lm,n , Lm,n ≡ Lm,n and Lm,n ≡ Lm,n . Observe that for any 0 < p < ∞, the Schattenp loss is orthogonally invariant and sum-decomposable, hence amenable to the
our analysis.
2

We thank the anonymous referee for this helpful suggestion.

7 Finding Optimal Shrinkers Numerically: Schatten norm losses

25

While it is in principle possible to derive the optimal shrinker for the Schattenp 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 {yi } in which we would like to evaluate the
optimal shrinker η ∗ (y). We define xi = x(yi ) where y 7→ x(y) is the transformation from (8). For each of the values {xi } we form a symbolic expression for the
function F (η, xi ) from (29), and minimize it numerically to obtain the minimizer
η ∗∗ (xi ) from (30). The desired value of the optimal shrinker η ∗ (yi ) is then given
by η ∗ (yi ) = η ∗∗ (x(yi )) = η ∗∗ (xi ).
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 ≥ 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 = ∞, namely
the operator norm). It seems that the optimal shrinker for all cases 1 ≤ p < ∞
are continuous, and that the discontinuity found analytically for the case p = ∞
forms only in the limit p → ∞. 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 = ∞ 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 β.
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 → 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.
β =1
4
3.5
3

β = 0.25
4

1
1. 1
1. 5
2
5
10000

3.5
3
2.5

η ∗ (y)

η ∗ (y)

2.5
2
1.5

2
1.5

1

1

0.5

0.5

0
1

1
1. 1
1. 5
2
5
10000

0
1.5

2

2.5

y

3

3.5

4

1

1.5

2

2.5

y

3

3.5

Figure 4: Numerically computed optimal shrinkers for the Schatten-p loss, for different
values of p ≥ 1 shown in the legend. Left, β = 1 (the case of square matrix); Right,
β = 0.25 (the case of four times as many columns as rows.) Each shrinker was evaluated
on a grid of 200 points. This figure can be reproduced using the code supplement [31].

4

8 Extensions

26
β =1

15

β = 0.25
15

0. 15
0. 25
0. 5
0. 75
1

η ∗ (y)

10

η ∗ (y)

10

0. 15
0. 25
0. 5
0. 75
1

5

5

0

0
2

4

6

8

y

10

12

14

2

4

6

8

y

10

12

14

Figure 5: Numerically computed optimal shrinkers for the Schatten-p loss, for different
values of 0 < p ≤ 1 shown in the legend. Left, β = 1 (the case of square matrix); Right,
β = 0.25 (the case of four times as many columns as rows.) Each shrinker was evaluated
on a grid of 200 points. This figure can be reproduced using the code supplement [31].

8 Extensions
Our main results
√ have been formulated and calibrated specifically for the model
Y = X + Z/ n ∈ Mm×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 + σZ ∈ Mm×n , and consider:
√
1. The setting where σ is either known but does not necessarily equal 1/ n,
or is altogether unknown.
2. The setting where the noise matrix Z has i.i.d entries, but its distribution is
not necessarily orthogonally invariant.
The results below follow [20].

8.1 Unknown Noise level
Consider an asymptotic framework
slightly more general than the one in Section
√
2.2, in which Yn = Xn +(σ/ n)Zn , with Xn and Zn as defined there. In this section
we keep the loss family L and the asymptotic aspect ratio β fixed and implicit. We
extend Definition 1 and write


a.s.
L∞ (η|x, σ) = lim Lmn ,n Xn , X̂η (Xn + √σn Zn ) .
n→∞

When the noise level σ is known, Eq. (11) allows
us to re-calibrate any nonlin√
earity η, originally calibrated for noise level 1/ n, to a different noise level. For a
nonlinearity η : [0, ∞) → [0, ∞), write
ηc (y) = c · η(y/c) .
We clearly have:

8 Extensions

27

√
Lemma 14. If η ∗ is an optimal shrinker for Yn = Xn + Zn / n, namely,
L∞ (η ∗ |x) ≤ L∞ (η|x)

for any r ≥ 1, any x ∈ Rr and any continuous
nonlinearity η, and σ > 0, then ησ∗ is an
√
optimal shrinker for Yn = Xn + (σ/ n)Zn , namely
L∞ (ησ∗ |x, σ) ≤ L∞ (η|x, σ)

for any r ≥ 1, any x ∈ Rr and any continuous nonlinearity η.

When the noise level σ is unknown, we are required to estimate it. See [15–17,
22,23,37] and references therein for existing literature on this estimation problem.
The method below has been proposed in [20].
Consider the following robust estimator for the parameter σ in the model Y =
X + σZ:
ymed
,
(51)
σ̂(Y ) = √
n · µβ

where ymed is a median singular value of Y and µβ is the median of the MarcenkoPastur distribution, namely, the unique solution in β− ≤ x ≤ β+ to the equation
Zx p
(β+ − t)(t − β− )
1
dt = ,
2πt
2
β−

√

where β± = 1 ± β. Note that the median µβ is not available analytically but can
easily be obtained by numerical quadrature.
√
Lemma 15. Let σ > 0. For the sequence Yn = Xn + (σ/ n)Zn in our asymptotic
framework,
σ̂(Yn ) a.s.
√ = σ.
lim
n→∞ 1/ n
√
our asymptotic framework and
Corollary 3. Let Yn = Xn + (σ/ n)Zn be a sequence in√
∗
let η be an optimal shrinker calibrated for Yn = Xn +Zn / n. Then the random sequence
∗
of shrinkers ησ̂(Y
converges to the optimal shrinker ησ∗ :
n)
a.s.

lim ησ̂(Yn ) (y) = ησ (y) ,

n→∞

y > 0.

∗
Consequently, ησ̂(Y
asymptotically achieves optimal performance:
n)


σ
√
∗
Z
)
= L∞ (ησ∗ |x, σ) .
lim Lmn ,n Xn , X̂ησ̂(Y
(X
+
n
n n
)
n→∞

n

In practice, for denoising a matrix Y ∈ Mm×n , assumed to satisfy Y = X + σZ,
where X is low-rank and Z has i.i.d entries, we have the following approximately
optimal singular value shrinkage estimator:
√
√
X̂(Y ) = nσ X̂η∗ (Y /( nσ))
(52)
when σ is known, and
X̂(Y ) =

√
∗

√
nσ̂(Y ) · X̂η∗ (Y /( nσ̂(Y )))

(53)

when σ is unknown. Here, η is an optimal shrinker with respect to desired loss
family L in the natural scaling.

8 Extensions

28

8.2 General white noise
Our results were formally stated for the sequence of models of the form Y =
X + σ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 ∈ Mm,m
and B ∈ Mn,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 Lm,n (X, X̂)
only depends on the signal matrix X through its nonzero singular values 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 [17].
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 Yn = Xn + Zn / n
satisfies the following assumptions:
1. General white noise: The entries of Zn are i.i.d samples from a distribution
with zero mean, unit variance and finite fourth moment.
2. Fixed signal column span and uniformly distributed signal singular vectors: Let
the rank r > 0 be fixed and choose a vector x ∈ Rr with coordinates x =
(x1 , . . . , xr ). Assume that for all n,
Xn = Un diag(x1 , . . . , xr , 0, . . . , 0) Vn′

(54)

is a singular value decomposition of Xn , where Un and Vn are uniformly
distributed random orthogonal matrices. Formally, Un and Vn are sampled
from the Haar distribution on the m-by-m and n-by-n orthogonal group,
respectively.
3. Asymptotic aspect ratio β: The sequence mn is such that mn /n → β.
The second assumption above implies that Xn is a “generic” choice of matrix
with nonzero singular values x, or equivalently, a generic choice of coordinate
systems in which the linear operator corresponding to X is expressed.
The results of [25], 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.

10 Conclusion

29

9 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.
Comparing the asymptotic loss with the empirical loss. 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 ≤ 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 η 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.

10 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

10 Conclusion

30

(20, 20, 1) unifor m

(20, 20, 1) gaus s ian
4.5

4.5

X̂ λ ∗
X̂ s ∗
X̂ o p t

4

4

4

3

2.5

2.5

2

2

1.5

1.5

1

1
1

1.5

2

2.5

3

x

3.5

4

4.5

3.5

3

5

MSE

MSE

MSE

4.5

3.5

3.5

2

1.5

1

1.5

2

2.5

3

x

3.5

4

4.5

1

5

3

2

2

2

1.5

1.5

2.5

3

x

3.5

4

4.5

5

1

3.5

4

4.5

5

3

1.5

2

3

x

2.5

2.5

1.5

2.5

3.5

MSE

MSE

2.5

2

X̂ λ ∗
X̂ s ∗
X̂ o p t

4

3.5

3

1

1.5

(100, 100, 1) s tudent-t(6)
4.5

X̂ λ ∗
X̂ s ∗
X̂ o p t

4

3.5

1

1

(100, 100, 1) unifor m
4.5

X̂ λ ∗
X̂ s ∗
X̂ o p t

4

3

2.5

(100, 100, 1) gaus s ian
4.5

MSE

(20, 20, 1) s tudent-t(6)

X̂ λ ∗
X̂ s ∗
X̂ o p t

1

1.5

2

2.5

3

x

3.5

4

4.5

5

1

1

1.5

2

2.5

3

x

3.5

4

4.5

5

Figure 6: Comparison of asymptotic (solid line) and empirically observed (dots)
Frobenius loss for n = 20, 100 and r = 1. Horizontal axis is the singular value
of the signal matrix X. Shown are optimally tuned soft threshold X̂λ∗ , optimally
tuned hard threshold X̂s∗ and optimal shrinker X̂opt from (7). This figure can be
reproduced using the code supplement [31].
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 η(y) = 0 for all y ≤ β+ , namely, shrinkers that
truncate data singular values below the bulk edge β+ . In some sense, this is the
class of “reasonable” shrinkers.
The challenging issue is how to control the manner in which “null” singular
values yn,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 [14] 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.

10 Conclusion

31

(50, 50, 2) unifor m

(50, 50, 2) gaus s ian
9

8

7

6

6

MSE

7

7

6

5

5

4

4

4

3

3

3

1

1.5

2

2.5

3

x

3.5

4

4.5

2

5

2
1

1.5

(50, 50, 4) gaus s ian

2

2.5

3

3.5

x

4

4.5

5

1

18

X̂ λ ∗
X̂ s ∗
X̂ o p t

16

12

MSE

14

12

MSE

14

12

10

10

10

8

8

8

6

6

6

4
1.5

2

2.5

3

x

3.5

4

4.5

5

2.5

3

x

3.5

4

4.5

5

X̂ λ ∗
X̂ s ∗
X̂ o p t

16

14

4

2

(50, 50, 4) s tudent-t(6)

18

X̂ λ ∗
X̂ s ∗
X̂ o p t

1

1.5

(50, 50, 4) unifor m

18

16

X̂ λ ∗
X̂ s ∗
X̂ o p t

8

5

2

MSE

(50, 50, 2) s tudent-t(6)
9

X̂ λ ∗
X̂ s ∗
X̂ o p t
MSE

8

MSE

9

X̂ λ ∗
X̂ s ∗
X̂ o p t

4
1

1.5

2

2.5

3

3.5

x

4

4.5

5

1

1.5

2

2.5

3

x

3.5

4

4.5

5

Figure 7: Comparison of asymptotic (solid line) and empirically observed (dots)
Frobenius loss for n = 50 and r = 2, 4. Horizontal axis is the singular value of
the signal matrix X. Shown are optimally tuned soft threshold X̂λ∗ , optimally
tuned hard threshold X̂s∗ and optimal shrinker X̂opt from (7). This figure can be
reproduced using the code supplement [31].
fro

op

nuc

3.5

4

3

3.5

2.5

2.5

3

2

1.5

η ∗ (y)

3

η ∗ (y)

η ∗ (y)

3.5

2

1.5

2.5

2

1

1

0.5

0.5

1

0

0

0.5

2

2.5

3

y

3.5

4

1.5

2

2.5

3

y

3.5

4

2

2.5

3

y

3.5

Figure 8: Optimal shrinkage computed by brute-force against the asymptotically optimal shrinkers from Section 6, for β = 1. This figure can be reproduced using the code
supplement [31].

Reproducible Research
In the code supplement [31] we offer a Matlab software library that includes:

4

References

32

1. 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.
2. Scripts that generate each of the figures in this paper.
3. 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
[1] Gene H. Golub and William Kahan. Calculating the Singular Values and
Pseudo-Inverse of a Matrix. Journal of the Society for Industrial & Applied Mathematics: Series B, 2(2):205–224, 1965.
[2] Carl Eckart and Gale Young. The approximation of one matrix by another of
lower rank. Psychometrika, 1(3), 1936.
[3] Orly Alter, Patrick O. Brown, and David Botstein. Singular value decomposition for genome-wide expression data processing and modeling. Proceedings
of the National Academy of Sciences, 97(18):10101–10106, aug 2000.
[4] Raymond B. Cattell. The scree test for the number of factors. Multivariate
Behavioral Research, 1(2):245–276, 1966.
[5] Donald A. Jackson. Stopping rules in principal components analysis: a comparison of heuristical and statistical approaches. Ecology, 1993.
[6] Terrence D. Lagerlund, Frank W. Sharbrough, and Neil E. Busacker. Spatial
filtering of multichannel electroencephalographic recordings through principal component analysis by singular value decomposition. Journal of clinical
neurophysiology : official publication of the American Electroencephalographic Society, 14(1):73–82, jan 1997.
[7] Alkes L. Price, Nick J. Patterson, Robert M. Plenge, Michael E. Weinblatt,
Nancy A. Shadick, and David Reich. Principal components analysis corrects for stratification in genome-wide association studies. Nature genetics,
38(8):904–9, aug 2006.

References

33

[8] Ove Edfors and Magnus Sandell. OFDM channel estimation by singular
value decomposition. IEEE Transactions on Communications, 46(7):931–939,
1998.
[9] David L. Donoho and Iain M. Johnstone. Ideal Spatial Adaptation by
Wavelet Shrinkage. Biometrika, 81(3):425–455, 1994.
[10] David L. Donoho and Iain M. Johnstone. Adapting to Unkown Smoothness via Wavelet Shrinkage. Journal of the American Statistical Association,
90(432):1200–1224, 1995.
[11] David L. Donoho, Iain M. Johnstone, Gerard Kerkyacharian, and Dominique
Picard. Wavelet shrinkage: Asymptopia? Journal of the Royal Statistical Society. Series B, 57(2):301–369, 1995.
[12] David L. Donoho. De-Noising by Soft-Thresholding. IEEE Transactions on
Information Theory, 41(3):613–627, may 1995.
[13] David L. Donoho and Iain M. Johnstone. Minimax estimation via wavelet
shrinkage. Annals of statistics, 26(3):879–921, 1998.
[14] David L. Donoho, Matan Gavish, and Iain M. Johnstone. Optimal shrinkage
of eigenvalues in the Spiked Covariance Model. arXiv:1311.0851, 2013.
[15] Art B. Owen and Patrick O. Perry. Bi-cross-validation of the SVD and the
nonnegative matrix factorization. The Annals of Applied Statistics, 3(2):564–
594, jun 2009.
[16] Patrick O. Perry. Cross validation for unsupervised learning. PhD thesis, Stanford University, 2009.
[17] Andrey A. Shabalin and Andrew B. Nobel. Reconstruction of a low-rank
matrix in the presence of Gaussian noise. Journal of Multivariate Analysis,
118:67–76, jul 2013.
[18] Jian-Feng Cai, Emmanuel J. Candès, and Zuowei Shen. A Singular Value
Thresholding Algorithm for Matrix Completion. SIAM Journal on Optimization, 20(4):1956, 2008.
[19] David L. Donoho and Matan Gavish. Minimax Risk of Matrix Denoising by
Singular Value Thresholding. Annals of Statistics, 42(6):2413–2440, 2014.
[20] Matan Gavish and David L. Donoho. The Optimal Hard Threshold for Singular Values is 4/sqrt(3). IEEE Transactions on Information Theory, 60(8):5040–
5053, 2014.
[21] Sourav Chatterjee. Matrix estimation by universal singular value thresholding. Annals of Statistics, 43(1):177–214, 2015.

References

34

[22] Raj Rao Nadakuditi. Optshrink: an algorithm for improved low-rank signal
matrix denoising by optimal, data-driven singular value shrinkage. IEEE
Transactions on Information Theory, 60(5):3002–3018, 2014.
[23] Julie Josse and Sylvain Sardy. Adaptive shrinkage of singular values. Statistics and Computing, 26:715–724, 2016.
[24] Iain M. Johnstone. On the distribution of the largest eigenvalue in principal
components analysis. Annals of Statistics, 29(2):295–327, 2001.
[25] Florent Benaych-Georges and Raj Rao Nadakuditi. The singular values and
vectors of low rank perturbations of large rectangular random matrices.
Journal of Multivariate Analysis, 111:120–135, oct 2012.
[26] Peter J. Bickel and Elizaveta Levina. Covariance regularization by thresholding. The Annals of Statistics, 36(6):2577–2604, dec 2008.
[27] Clifford Lam and Jianqing Fan. Sparsistency and Rates of Convergence in
Large Covariance Matrix Estimation. Annals of statistics, 37(6B):4254–4278,
jan 2009.
[28] Tony Cai, Cun-Hui Zhang, and Harrison H. Zhou. Optimal rates of convergence for covariance matrix estimation. The Annals of Statistics, 38(4):2118–
2144, aug 2010.
[29] Angelika Rohde and Alexandre B. Tsybakov.
Estimation of highdimensional low-rank matrices. The Annals of Statistics, 39(2):887–930, 2011.
[30] Vladimir Koltchinskii, Karim Lounici, and Alexandre B. Tsybakov. Nuclearnorm penalization and optimal rates for noisy low-rank matrix completion.
The Annals of Statistics, 39(5):2302–2329, oct 2011.
[31] Matan Gavish and David Donoho. Code supplement to ”Optimal Shrinkage
of Singular Values” http://purl.stanford.edu/kv623gt2817, 2015.
[32] Greg W. Anderson, Alice Guionnet, and Ofer Zeitouni. An Introduction to
Random Matrices. Cambridge university press, Cambridge, 2010.
[33] Y Q Yin, Zhidong Bai, and P R Krishnaiah. On the limit of the largest eigenvalue of the large dimensional sample covariance matrix. Probability Theory
and Related Fields, 78:509–521, 1988.
[34] Zhidong Bai and YQ Yin. Limit of the smallest eigenvalue of a large dimensional sample covariance matrix. The annals of Probability, 21(3):1275–1294,
1993.
[35] Debashis Paul. Asymptotics of Sample Eigenstructure for a Large Dimensional Spiked Covariance Model. Statistica Sinica, 17(4):1617–1642, 2007.

References

35

[36] R. Brent Dozier and Jack W. Silverstein. On the empirical distribution
of eigenvalues of large dimensional information-plus-noise-type matrices.
Journal of Multivariate Analysis, 98(4):678–694, apr 2007.
[37] Shira Kritchman and Boaz Nadler. Non-parametric detection of the number
of signals: Hypothesis testing and random matrix theory. Signal Processing,
IEEE Transactions, 57(10):3930–3941, 2009.

