arXiv:1703.05038v2 [math.NA] 27 Feb 2018

Harmonic Mean Iteratively Reweighted Least Squares for
Low-Rank Matrix Recovery
Christian Kümmerle∗, Juliane Sigl†
February 28, 2018

Abstract
We propose a new iteratively reweighted least squares (IRLS) algorithm for the recovery of
a matrix X ∈ Cd1 ×d2 of rank r  min(d 1 , d 2 ) from incomplete linear observations, solving a
sequence of low complexity linear problems. The easily implementable algorithm, which we
call harmonic mean iteratively reweighted least squares (HM-IRLS), optimizes a non-convex
Schatten-p quasi-norm penalization to promote low-rankness and carries three major strengths,
in particular for the matrix completion setting. First, we observe a remarkable global convergence behavior of the algorithm’s iterates to the low-rank matrix for relevant, interesting cases,
for which any other state-of-the-art optimization approach fails the recovery. Secondly, HMIRLS exhibits an empirical recovery probability close to 1 even for a number of measurements
very close to the theoretical lower bound r (d 1 + d 2 − r ), i.e., already for significantly fewer linear
observations than any other tractable approach in the literature. Thirdly, HM-IRLS exhibits a
locally superlinear rate of convergence (of order 2 − p) if the linear observations fulfill a suitable
null space property. While for the first two properties we have so far only strong empirical
evidence, we prove the third property as our main theoretical result.

1

Introduction

The problem of recovering a low-rank matrix from incomplete linear measurements or observations
has gained considerable attention in the last few years due to the omnipresence of low-rank models
in different areas of science and applied mathematics. Low-rank models arise in a variety of areas such as system identification [LHV13, LV10], signal processing [AR15], quantum tomography
[GLF+ 10, Gro11] and phase retrieval [CSV13, CESV13, GKK15]. An instance of this problem of particular importance, e.g., in recommender systems [SRJ05, GNOT92, CR09], is the matrix completion
problem, where the measurements correspond to entries of the matrix to be recovered.
Although the low-rank matrix recovery problem is NP-hard in general, several tractable algorithms have been proposed that allow for provable recovery in many important cases. The nuclear
norm minimization (NNM) approach [Faz02, CR09], which solves a surrogate semidefinite program,
is particularly well-understood. For NNM, recovery guarantees have been shown for a number of
∗ Department of Mathematics, Technische Universität München, Boltzmannstr. 3, 85748 Garching/Munich, Germany
E-mail: christian.kuemmerle@ma.tum.de
† Department of Mathematics, Technische Universität München, Boltzmannstr. 3, 85748 Garching/Munich, Germany
E-mail: juliane.sigl@ma.tum.de

1

measurements on the order of the information theoretical lower bound r (d 1 + d 2 − r ), if r denotes
the rank of a d 1 × d 2 -matrix [RFP10, CR09]; i.e., for a number of measurements m ≥ ρr (d 1 + d 2 − r )
with some oversampling constant ρ ≥ 1. Even though NNM is solvable in polynomial time, it can
be computationally very demanding if the problem dimensions are large, which is the case in many
potential applications. Another issue is that although the number of measurements necessary for
successful recovery by nuclear norm minimization is of optimal order, it is not optimal. More precisely, it turns out that the oversampling factor ρ of nuclear norm minimization has to be much
larger than the oversampling factor of some other, non-convex algorithmic approaches [ZL15, TW13].
These limitations of convex relaxation approaches have led to a rapidly growing line of research
discussing the advantages of non-convex optimization for the low-rank matrix recovery problem
[JMD10, TW13, HH09, JNS13, WYZ12, TW16, Van13, WCCL16, TBS+ 15]. For several of these nonconvex algorithmic approaches, recovery guarantees comparable to those of NNM have been derived
[CLS15, TBS+ 15, ZL15, SL16]. Their advantage is a higher empirical recovery rate and an often
more efficient implementation. While there are some results about global convergence of firstorder methods minimizing a non-convex objective [GLM16, BNS16] so that a success of the method
might not depend on a particular initialization, the assumptions of these results are not always
optimal, e.g., in the scaling of the numbers of measurements m in the rank r [GLM16, Theorem 5.3].
In general, the success of many non-convex optimization approaches relies on a distinct, possibly
expensive initialization step.

1.1

Contribution of this paper

In this spirit, we propose a new iteratively reweighted least squares (IRLS) algorithm for the lowrank matrix recovery problem1 that strives to minimize a non-convex objective function based on
the Schatten-p quasi-norm
p
min kX kSp subject to Φ(X ) = Y ,
(1)
X

for 0 < p < 1, where Φ : Cd1 ×d2 → Cm is the linear measurement operator and Y ∈ Cm is the
data vector that define the problem. The overall strategy of the proposed IRLS algorithm is to mimic
this minimization by a sequence of weighted least squares problems. This strategy is shared by
the related previous algorithms of Fornasier, Rauhut & Ward [FRW11] and Mohan & Fazel [MF12]
which minimize (1) by defining iterates as
1

X (n+1) = min kWL(n) 2 X kF2 subject to Φ(X ) = Y ,
X

(2)

p−2

where WL(n) ≈ (X (n)X (n)∗ ) 2 is a so-called weight matrix which reweights the quadratic penalty by
operating on the column space of the matrix variable. Thus, we call this column-reweighting type
of IRLS algorithms IRLS-col. Due to the inherent symmetry, it is evident to conceive, still in the
spirit of [FRW11, MF12], the algorithm IRLS-row
1

X (n+1) = min kWR(n) 2 X ∗ kF2 subject to Φ(X ) = Y
X

1 The

(3)

algorithm and partial results were presented at the 12th International Conference on Sampling Theory and
Applications in Tallinn, Estonia, July 3–7, 2017. The corresponding conference paper has been published in its proceedings
[KS17].

2

p−2

with WR(n) ≈ (X (n)∗X (n) ) 2 , which reweights the quadratic penalty by acting on the row space of the
matrix variable. We note that even for square dimensions d 1 = d 2 , IRLS-col and IRLS-row do
not coincide.
In this paper, as an important innovation, we propose the use of a different type of weight matrices, so-called harmonic mean weight matrices, which can be interpreted as the harmonic mean of the
matrices WL(n) and WR(n) above. This motivates the name harmonic mean iteratively reweighted least
squares (HM-IRLS) for the corresponding algorithm. The harmonic mean of the weight matrices
of IRLS-col and of IRLS-row in HM-IRLS is able to use the information in both the column
and the row space of the iterates, and it also gives rise to a qualitatively better behavior than the use
of more obvious symmetrizations as, e.g., the arithmetic mean of weight matrices would allow for,
both in theory and in practice.
We argue that the choice of harmonic mean weight matrices as in HM-IRLS leads to an efficient algorithm for the low-rank matrix recovery problem with fast convergence and superior performance in terms of sample complexity, also compared to algorithms based on strategies different
from IRLS.
On the one hand, we show that the accumulation points of the iterates of HM-IRLS converge
to stationary points of a smoothed Schatten-p functional under the linear constraint, as it is known
for, e.g., IRLS-col, c.f. [FRW11, MF12]. On the other hand, we extend the theoretical guarantees
which are based on a Schatten-p null space property (NSP) of the measurement operator [OMFH11,
FR13], to HM-IRLS.
Our main theoretical result is that HM-IRLS exhibits a locally superlinear convergence rate of
order 2 − p in the neighborhood of a low-rank matrix for the non-convexity parameter 0 < p < 1
connected to the Schatten-p quasinorm, if the measurement operator fulfills the mentioned NSP of
sufficient order. For p  1, this means that the convergence rate is almost quadratic.
Although parts of our theoretical results, as in the case of the IRLS algorithms algorithms of
Fornasier, Rauhut & Ward [FRW11] and Mohan & Fazel [MF12], do not apply to the matrix completion setting, due to the popularity of the problem and for reasons of comparability with other
algorithms, we conduct numerical experiments to explore the empirical performance of HM-IRLS
also for this setting. Surprisingly enough we observe that the theoretical results comply with our
numerical experiments also for matrix completion. In particular, the theoretically predicted local
convergence rate of order 2 − p can be observed very precisely for this important measurement
model as well (see Figures 3 to 5).
This local superlinear convergence rate is unprecedented for IRLS variants such as IRLS-col
and as those that use the arithmetic mean of the one-sided weight matrices: this means that neither
can a superlinear rate be verified numerically, nor is it possible to show such a rate by our proof
techniques for any other IRLS variant.
To the best of our knowledge, HM-IRLS is the first algorithm for low-rank matrix recovery
which achieves superlinear rate of convergence for low complexity measurements as well as for
larger problems.
Additionally, we conduct extensive numerical experiments comparing the efficiency of HMIRLS with previous IRLS algorithms as IRLS-col, Riemannian optimization techniques [Van13],
alternating minimization approaches [HH09, TW16], algorithms based on iterative hard thresholding [KC14, BTW15], and others [PKCS16], in terms of sample complexity, again for the important
case of matrix completion.
The experiments lead to the following observation: HM-IRLS recovers low-rank matrices sys3

tematically with an optimal number of measurements that is very close to the theoretical lower
bound on the number of measurements that is necessary for recovery with high empirical probability. We consider this result to be remarkable, as it means that for problems of moderate dimensionality (matrices of ≈ 107 variables, e.g. (d 1 × d 2 )-matrices with d 1 ≈ d 2 ≈ 5 · 103 ) the proposed
algorithm needs fewer measurements for the recovery of a low rank matrix than all the state-of-the-art
algorithms we included in our experiments (see Figure 6).
An important practical observation of HM-IRLS is that its performance is very robust to the
choice of the initialization and can be used as a stand-alone algorithm to recover low-rank matrices
also starting from a trivial initialization. This is suggested by our numerical experiments since
even for random or adversary initializations, HM-IRLS converges to the low-rank matrix, even
though it is based on an objective function which is highly non-convex. While a complete theoretical
understanding of this behavior is not yet achieved, we regard the empirical evidence in a variety of
interesting cases as strong. In this context, we consider a proof of the global convergence of HMIRLS for non-convex penalizations under appropriate assumptions as an interesting open problem.

1.2

Outline

We proceed in the paper as follows. In the next section, we provide some background on Kronecker
and Hadamard products of matrices as these concepts are used in the analysis of the algorithm to be
discussed. Moreover, we explain different reformulations of the Schatten-p quasi-norm in terms of
weighted `2 -norms, which lead to the derivation of the harmonic mean iteratively reweighted least
squares (HM-IRLS) algorithm in Section 3. We present our main theoretical results, the convergence guarantees and the locally superlinear convergence rate for the algorithm in Section 4. Numerical experiments and comparisons to state-of-the-art methods for low-rank matrix recovery are
carried out in Section 5. In Section 6, we interpret the algorithm’s different steps as minimizations
of an auxililary functional with respect to its arguments and show theoretical guarantees for HMIRLS extending similar guarantees for IRLS-col. After this, we detail the proof of the locally
superlinear convergence rate under appropriate assumptions on the null space of the measurement
operator.

2
2.1

Notation and background
General notation, Schatten-p and weighted norms

In this section, we explain some of the notation we use in the course of this paper.
The set of matrices X ∈ Cd1 ×d2 is denoted by Md1 ×d2 . Unless stated otherwise, vectors x ∈ Cd are
h
iT
considered as column vectors. We also use the vectorized form X vec = X 1T , . . . , X jT , . . . , X dT2 ∈

Cd1d2 of a matrix X ∈ Md1 ×d2 with columns X j , j ∈ {1, . . . , d 2 }. The
 reverse recast of a vector
d
d
1
2
x ∈C
into a matrix of dimension d 1 × d 2 is denoted by x mat(d1,d2 ) = X 1 , . . . , X j , . . . , X d2 , where
X j = [x (d1 −1)·j+1 , . . . , x (d1 −1)·j+d1 ]T , j = 1, . . . , d 2 are column vectors, or X mat if the dimensions are
clear from the context. Obviously, it holds that X = (X vec )mat .
The identity matrix in dimension d×d is denoted by Id . With 0d1 ×d2 ∈ Md1 ×d2 and 1d1 ×d2 ∈ Md1 ×d2
we denote the matrices with only 0- or 1-entries respectively. The set of Hermitian matrices is
denoted by Hd ×d := {X ∈ Md ×d | X = X ∗ }. We write X + ∈ Md1 ×d2 for the Moore-Penrose inverse
of the matrix X ∈ Md1 ×d2 .
4

Let Ud = {U ∈ Cd ×d ; U U ∗ = Id } denote the set of unitary matrices. Then the singular value
decomposition of a matrix X ∈ Md1 ×d2 can be written as X = U ΣV ∗ with U ∈ Ud1 , V ∈ Ud2 and
Σ ∈ Md1 ×d2 , where Σ is diagonal and contains the singular values of X such that Σii = σi (X ) ≥ 0 for
i ∈ {1, . . . , min(d 1 , d 2 )}. We define the Schatten-p (quasi-)norm of X ∈ Md1 ×d2 as

kX kSp :=



rank(X ),


 hÍ


i 1/p
min(d 1,d 2 ) p
σ
(X
)
j
j=1




σmax (X ),


for p = 0,
,

for 0 < p < ∞,

(4)

for p = ∞.

Note that for p = 1, the Schatten-p norm is also called nuclear norm, written as kX k∗ := kX kS1 .
The trace tr[X ] of a matrix X ∈ Md1 ×d2 is defined by the sum of its diagonal elements, tr[X ] =
Ímin(d1,d2 )
p
X j j . It can be seen that the p-th power of the Schatten-p norm coincides with kX kSp =
j=1


tr (X ∗X )p/2 . The Schatten-2 norm is also called Frobenius norm and has the property that it is
p
induced by the Frobenius scalar product hX , Y iF = tr [X ∗Y ], i.e., kX kF = kX kS2 = hX , X iF . We
define the weighted Frobenius scalar product of two matrices X , Y ∈ Md1 ×d2 weighted by the the
positive definite weight matrix W ∈ Hd1 ×d1 as hX , Y iF (W ) := hW X , Y iF = hX ,W Y iF . This scalar
p
p
product induces the weighted Frobenius norm kX kF (W ) = hX , X iF (W ) = tr[(W X )∗X ]. It is clear
that the Frobenius norm of a matrix X coincides with the `2 -norm of its vectorization X vec , i.e.,
kX kF = kX vec k`2 .
Similar to weighted Frobenius norms, we define the weighted `2 -scalar product of vectors x, y ∈
Cd weighted by the positive definite weight matrix W ∈ Hd ×d as hx, yi`2 (W ) = x ∗W y = y ∗W x and
√
its induced weighted `2 -norm as kx k`2 (W ) = x ∗W x. We use the notation X  0 for a positive
definite matrix
X ∈ Hd ×d . Furthermore, we denote the range of a linear map Φ : Md1 ×d2 → Cm by

Ran(Φ) = Y ∈ Cm ; there is X ∈ Md1 ×d2 such that Y = Φ(X ) and its null space by N (Φ) = X ∈
Md1 ×d2 ; Φ(X ) = 0 .

2.2

Problem setting and characterization of Sp - and reweighted Frobenius norm
minimizers

Given a linear map Φ : Md1 ×d2 → Cm such that m  d 1d 2 , we want to uniquely identify and
reconstruct an unknown matrix X 0 from its linear image Y := Φ(X 0 ) ∈ Cm . However, basic linear
algebra tells us that this is not possible without further assumptions, since Φ is not injective if
m < d 1d 2 . Indeed, there is a (d 1d 2 − m)-dimensional affine space {X 0 } + N (Φ) fulfilling the linear
constraint
Φ(X ) = Y .
Nevertheless, under the additional assumption that the matrix X 0 ∈ Md1 ×d2 has rank r < min(d 1 , d 2 )
and under appropriate assumptions on the map Φ, the recovery of X 0 is possible by solving the affine
rank minimization problem
min rank(X ) subject to Φ(X ) = Y .

(5)

The unique solvability of (5) is given with high probability if, for example, Φ is a linear map whose
matrix representation has i.i.d. Gaussian entries [ENP12] and m = Ω(r (d 1 + d 2 )). Unfortunately,
solving (5) is intractable in general, but the works [CR09, RFP10, CP11] suggest solving the tractable
5

convex optimization program
min kX kS1 subject to Φ(X ) = Y ,

(6)

also called nuclear norm minimization (NNM), as a proxy.
As discussed in the introduction, there are empirical as well as theoretical results (e.g., in [DDFG10,
Cha07]) coming from the related sparse vector recovery problem that suggest alternative relaxation
approaches. These results indicate that it might be even more advantageous to solve the non-convex
problem
p
min F p (X ) := kX kSp subject to Φ(X ) = Y ,
(7)
for 0 < p < 1, i.e., minimizing the p-th power of the Schatten-p quasi-norms under the affine
constraint. Heuristically, the choice of p < 1 relatively small can be motivated by the observation
that by the definition (4) of the Schatten-p quasi-norm
p→0

p

kX kSp −−−→ rank(X ) =: kX kS0 .
The above consideration suggests that the solution of (7) might be closer to (5) than (6) for small
p. On the other hand, again, it is in general computationally intractable to find a global minimum
of the non-convex optimization problem (7) if p < 1. Therefore it is a natural and very relevant
question to ask which optimization algorithm to use to find global minimizers of (7).
In this paper, we discuss an algorithm striving to solve (7) that is based on the following observations: Assume for the moment that we are given a square matrix X ∈ Md1 ×d2 with d 1 = d 2 of full
rank. Then, we can rewrite the p-th power of its Schatten-p quasi-norm as a square of a weighted
Frobenius norm, or, using Kronecker product notiation as explained in Appendix A, as a square of
a weighted `2 -norm (if we use the vectorized notation X vec ): Iit turns out that
p

p

(i) kX kSp = tr[(XX ∗ ) 2 ] = tr[(XX ∗ )

p−2
2

1

(XX ∗ )] = tr(WL XX ∗ ) = kWL2 X kF2

1

= kX kF2 (WL ) = k(Id2 ⊗ WL ) 2 X vec k`22 = kX vec k`22 (Id
p−2
2

2

⊗WL ) ,

whereWL is the symmetric weight matrix (XX ∗ ) in Md1 ×d1 and Id2 ⊗WL is the block diagonal
weight matrix in Md1d2 ×d1d2 with d 2 instances of WL on the diagonal blocks, but also that
p

p

(ii) kX kSp = tr[(X ∗X ) 2 ] = tr[(X ∗X )(X ∗X )

p−2
2

1

] = tr(X ∗XWR ) = kXWR2 kF2

1

= kX ∗ kF2 (WR ) = k(WR ⊗ Id1 ) 2 X vec k`22 = kX vec k`22 (WR ⊗Id ) ,
p−2

1

where WR is the symmetric weight matrix (X ∗X ) 2 in Md2 ×d2 . It follows from the definition
of the Kronecker product that the weight matrix WR ⊗ Id1 ∈ Md1d2 ×d1d2 is a block matrix of
diagonal blocks of the type diag (WR )i j , . . . , (WR )i j ∈ Md1 ×d1 , i, j ∈ [d 2 ].

6

(a) Id2 ⊗ WL

(b) WR ⊗ Id1

Figure 1: Sparsity structure of the weight matrices ∈ Md1d2 ×d1d2
The sparsity structures of Id2 ⊗ WL and WR ⊗ Id1 are illustrated in Fig. 1. Note that a representation
1

p

of kX kSp by squares of Frobenius norms can be achieved by multiplying X by WL2 from the left in
1

(i), or by WR2 from the right in (ii).
The above calculations are not well-defined if X is not of full rank or if d 1 , d 2 , since in these
cases at least one of the matrices XX ∗ ∈ Md1 ×d1 or X ∗X ∈ Md2 ×d2 is singular, prohibiting the
p−2
p−2
definition of the matrices WR = (X ∗X ) 2 or WL = (XX ∗ ) 2 for p < 2. However, these issues can be
overcome by introducing a smoothing parameter ϵ > 0 and smoothed weight matrices WL (X , ϵ) ∈
Md1 ×d1 and WR (X , ϵ) ∈ Md2 ×d2 defined by
WL (X , ϵ) := (XX ∗ + ϵ 2 Id1 )

p−2
2

,

(8)

WR (X , ϵ) := (X ∗X + ϵ 2 Id2 )

p−2
2

.

(9)

Remark 1. The weight matrices WL (X , ϵ) and WR (X , ϵ) are symmetric and positive definite.
The possibility to rewrite the p-th power of the Schatten-p of a matrix as a weighted Frobenius
norm gives rise to the general strategy of IRLS algorithms for low-rank matrix recovery: Weighted
least squares problems of the type
min

X ∈Md1 ×d2
Φ(X )=Y

kX kF2 (WL )

or

min

X ∈Md1 ×d2
Φ(X )=Y

kX ∗ kF2 (WR )

are solved and weight matrices WL are updated alternatingly, leading to the algorithms columnreweighting IRLS-col and row-reweighting IRLS-row, respectively [MF12, FRW11].

2.3

Averaging of weight matrices

While the algorithms IRLS-col and IRLS-row provide a tractable local minimization strategy
of smoothed Schatten-p functionals under the linear constraint, we argue that it is suboptimal to
follow either one of the two approaches as they do not exploit the symmetry of the problem in an
optimal way: They either use low-rank information in the column space or in the row space.
A first intuitive approach towards a symmetric exploitation of the low-rank structure is inspired
by the following identity, by combing the calculations (i) and (ii) carried out in Section 2.2.
7

Lemma 2. Let 0 < p ≤ 2 and X ∈ Md1 ×d2 with d = d 1 = d 2 be a full rank matrix. Then
p
kX kSp


1

1
1
1
WL ⊕ WR 2
2
2
2 2
=
X vec
kWL X kF + kXWR kF =
2
2

2

`2

= kX vec k`22 (W(arith) ) ,

where

 WL ⊕ WR
1
=: W(arith)
Id2 ⊗ WL + WR ⊗ Id1 =
2
2
is the arithmetic mean matrix of the symmetric and positive definite weight matrices Id2 ⊗ WL and
p−2
p−2
WR ⊗ Id1 , WL := (XX ∗ ) 2 , and WR := (X ∗X ) 2 .
Unfortunately, the introduction of arithmetic mean weight matrices does not prove to be particularly advantageous compared to one-sided reweighting strategies. No convincing improvements
can be noted neither in numerical experiments nor in the theoretical investigations for the convergence rate of IRLS for low-rank matrix recovery, cf. also Section 5.2 and Remark 22.
In contrast, we want to promote the usage of the harmonic mean of the weight matrices Id2 ⊗ WL
 −1
 −1
and WR ⊗ Id1 , i.e., weight matrices of the type 2 WR−1 ⊗ Id1 + Id2 ⊗ WL−1
= 2 WL−1 ⊕ WR−1
=:
W(harm) . In the remaining parts of the paper, we explain why W(harm) is able to significantly outperform other weighting variants both theoretically and practically.
The following lemma verifies that also the harmonic mean of the weight matrices Id2 ⊗ WL and
WR ⊗ Id1 leads to a legitimate reformulation of the Schatten-p quasi-norm power.
Lemma 3. Let 0 < p ≤ 2 and X ∈ Cd1 ×d2 with d = d 1 = d 2 be a full rank matrix. Then
p

kX kSp = 2 WL−1 ⊕ WR−1
where
2 WR−1 ⊗ Id1 + Id2 ⊗ WL−1

 − 21

 −1

2

X vec

`2

= kX vec k`22 (W(harm) ) ,

= 2 WL−1 ⊕ WR−1

 −1

=: W(harm)

is the harmonic mean matrix of the symmetric and positive definite weight matrices Id2 ⊗ WL and
p−2
p−2
WR ⊗ Id2 , WL := (XX ∗ ) 2 and WR := (X ∗X ) 2 .
Í
Proof. Let X = U ΣV ∗ = di=1 σi ui vi∗ ∈ Md ×d be the singular value decomposition of X . Therefore
for the vectorized version, X vec = (V ⊗ U )Σvec holds true. By the definitions of WL and WR , we can
Í
Í
2−p
2−p
write WL−1 = di=1 σi ui ui∗ and WR−1 = di=1 σi vi vi∗ . Using the Kronecker sum inversion formula
of Lemma 23 in Appendix A, we obtain
1

2
X vec k`22 = 2 WL−1 ⊕ WR−1
kX vec k`22 (W(harm) ) = kW(harm)

∗

 −1
= 2 tr WL−1 ⊕ WR−1 X vec
X

 − 21

X vec

2
`2

mat

=

d Õ
d
d Õ
Õ

2σk

2−p
i=1 j=1 k =1 σi

=2

d
Õ
σi2

2−p
i=1 2σi

!

v v ∗v u ∗ u u ∗
2−p j i k k i i

+ σj

d
Õ
l =1

σl ul vl∗

p

= kX kSp ,

which finishes the proof.


8

3

Harmonic mean iteratively reweighted least squares algorithm

In this section, we use this idea to formulate a new iteratively reweighted least squares algorithm
for low-rank matrix recovery. The so-called harmonic mean iteratively reweighted least squares algorithm (HM-IRLS) solves a sequence of weighted least squares problems to recover a low-rank
matrix X 0 ∈ Md1 ×d2 from few linear measurements Φ(X 0 ) ∈ Cm . The weight matrices appearing in
the least squares problems can be seen as the harmonic mean of the weight matrices in (8) and (9),
i.e., the ones used by IRLS-col and IRLS-row.
More precisely, for 0 < p ≤ 1 and d = min(d 1 , d 2 ), D = max(d 1 , d 2 ), given a non-increasing
∞ and the sequence of iterates (X (n) )∞ produced by
sequence of non-negative real numbers (ϵ (n) )n=1
n=1
the algorithm, we update our weight matrices such that
h
i −1
e (n) = 2 U (n) (Σd(n) )2−p U (n)∗ ⊕ V (n) (Σd(n) )2−p V (n)∗
W
,
1
2

(10)

(n)

with the diagonal matrices Σdt ∈ Mdt ×dt for dt = {d 1 , d 2 } such that
(n)
(Σdt )ii

1

(σi (X (n) )2 + ϵ (n)2 ) 2
=
0
(

if i ≤ d,
if d < i ≤ D,

(11)

and the matrices U (n) ∈ Ud1 and V (n) ∈ Ud2 , containing the left and right singular vectors of X (n)
in its columns, respectively.
e (n) can be seen as a stabilized version of the harmonic mean
We note that this definition of W
e (n) becomes very ill-conditioned
weight matrixW(harm) of Lemma 3. This stabilization is necessary as W
2−p
as soon as some of the singular values of X (n) approach zero and, related to that, (X (n)X (n)∗ ) 2 ⊕
2−p
(X (n)∗X (n) ) 2 would even be singular as soon as X (n) is not of full rank.
Additionally, for the formualtion of the algorithm and any n ∈ N, it is convenient to define the
f(n) )−1 : Md1 ×d2 → Md1 ×d2 as
linear operator (W
f(n) )−1 (X ) :=
(W

i
1 h (n) (n) 2−p (n)∗
(n)
U (Σd1 ) U X + XV (n) (Σd2 )2−p V (n)∗ ,
2

e (n) on Md1 ×d2 .
describing the operation of the inverse of W
Finally, HM-IRLS can be formulated in pseudo code as follows.

9

(12)

Algorithm 1 Harmonic Mean IRLS for low-rank matrix recovery (HM-IRLS)
Input: A linear map Φ : Md1 ×d2 → Cm , image Y = Φ(X 0 ) of the ground truth matrix X 0 ∈ Md1 ×d2 ,
rank estimate e
r , non-convexity parameter 0 < p ≤ 1.
0
Output: Sequence (X (n) )nn=1
⊂ Md1 ×d2 .
(0)
e
Initialize n = 0, ϵ = 1 and W (0) = Id1d2 ∈ Md1d2 ×d1d2 .
repeat

f(n) −1 Φ∗ (Φ ◦ (W
f(n) )−1 ◦ Φ∗ )−1 (Y ) ,
X (n+1) = arg min kX vec k`2 (W
f(n) )= (W )

(13)



ϵ (n+1) = min ϵ (n) , σer +1 (X (n+1) ) ,
h
i −1
e (n+1) = 2 U (n+1) (Σd(n+1) )2−p U (n+1)∗ ⊕ V (n+1) (Σd(n+1) )2−p V (n+1)∗
W
,
1
2

(14)

Φ(X )=Y

2

(15)

where U (n+1) ∈ Ud1 and V (n+1) ∈ Ud2 are matrices containing the left and right singular vectors
(n+1)

of X (n+1) in its columns, and the Σdt

are defined for t ∈ {1, 2} according to (11).

n = n + 1,
until stopping criterion is met.
Set n 0 = n.
From a practical point of view, it is beneficial that the explicit calculation of the very large
e (n+1) ∈ Hd1d2 ×d1d2 from (15) is not necessary in implementations of Algorithm 1.
weight matrices W
e (n+1) )−1
As suggested by formulas (12) and (13), it can be seen that just the operation of its inverse (W
(n)
−1
e ) is needed, which can be implemented by matrix-matrix multiplications on the space
resp. (W
e ∈ Md1 ×d2 , we have that W
e (n)X vec = X
evec if and only if X vec = (W
e (n) )−1X
evec ,
Md1 ×d2 : For matrices X , X
which can be written in matrix variables as
1 h (n) (n) 2−p (n)∗ e e (n) (n) 2−p (n)∗ i
X =
U (Σd1 ) U X + XV (Σd2 ) V
.
2
e (n) and the Kronecker sum, cf. (15) and Appendix A.
The last equivalence is due to the definitions of W
Note that the smoothing parameters ϵ (n) are chosen in dependence on a rank estimate r˜ here,
which will be an important ingredient for the theoretical analysis of the algorithm. In practice, how∞ are possible
ever, other choices of non-increasing sequences of non-negative real numbers (ϵ (n) )n=1
and can as well lead to (a maybe even faster) convergence when tuned appropriately.
We refer to Section 5.4 for a further discussion of implementation details.
Example With a simple example, we illustrate the versatility of HM-IRLS: Let d 1 = d 2 = 4, and
assume that we want to reconstruct the rank-1 matrix
1
2
3
4
1
© ª
©
ª

­ 10 20 30 40 ®
­ 10 ®
X 0 = uv ∗ = ­ ® 1 2 3 4 = ­
®
­ −2 −4 −6 −8 ®
­ −2 ®
«0.1 0.2 0.3 0.4¬
«0.1¬
from m = d f = r (d 1 + d 2 − r ) = 7 sampled entries Φ(X 0 ), where Φ is the linear map Φ : M 4×4 →
C7 , Φ(X ) = X 2,1 , X 4,1 , X 3,2 , X 4,2 , X 4,3 , X 1,4 , X 2,4 . Since the linear map Φ samples some
10

entries of matrices in M 4×4 and does not see the others, this is an instance of the problem that is
called matrix completion.
In general, reconstructing a (d 1 × d 2 ) rank-r matrix from m = r (d 1 + d 2 − r ) entries is a hard
problem, as it is known that if m < r (d 1 + d 2 − r ), there is always more than one matrix X such
that Φ(X ) = Φ(X 0 ), and even for equality, the property that Φ is invertible on (most) rank-r matrices
might be hard to verify [KTT15].
It can be argued that the specific matrix completion problem we consider is in some sense a
hard one, since, e.g., the deterministic sufficient condition for unique completability of [PABN16,
Theorem 2] is not fulfilled (less then 2 observed entries in the third column), and since the classical
kuu ∗ e k 2
kvv ∗ e k 2
coherence parameters µ(u) = d 1 max ku ki4 2 ≈ 3.81 and µ(v) = d 2 max kv k 4i 2 ≈ 2.13 that are
1≤i ≤4

1≤i ≤4

2

2

used to analyze the behavior of many matrix completion algorithms [CR09, JNS13] are quite large,
with µ(u) being quite close to the maximal value of 4.
On the other hand, as the problem is small and X 0 has rank r = 1, it is possible to impute the
missing values of
∗
∗
∗
4
©
ª
∗ 40®
­ 10 ∗
­
®
∗®
­ ∗ −4 ∗
«0.1 0.2 0.3 ∗ ¬

by solving very simple linear equations, since, for example, X 4,4 = u 4v 4 , X 2,1 = u 2v 1 , X 2,4 = u 2v 4 , and
X X
X 4,1 = u 4v 1 , and therefore X 4,4 = 4,X12, 12, 4 = 0.4. This shows that the only rank-1 matrix compatible
with Φ(X 0 ) is X 0 .
It turns out that—without using the combinatorial simplicity of the problem—the classical NNM
does not solve the problem, as the nuclear norm minimizer (solution of (6) for Y = Φ(X 0 )) produced
by the semidefinite program of the convex optimization package CVX [GB14] converges to

X nuclear

1
0.023 0.041
4
©
ª
0.232 0.411
40 ®
­ 10
≈­
®,
­−0.056 −4 −0.200 −0.226®
0.2
0.3
0.400 ¬
« 0.1

a matrix with 45.74 ≈ kX nuclear kS1 < kX 0 kS1 = σ1 (X 0 ) ≈ 56.13 and a relative Frobenius error of
kX nuclear −X 0 kF
kX 0 kF

= 0.661.
Interestingly, HM-IRLS is able to solve the problem, if p is chosen small enough, with very
high precision already after few iterations, for example, up to a relative error of 4.18 · 10−13 after
24 iterations if p = 0.1. This is in contrast to the behavior of IRLS-col, IRLS-row and also
to AM-IRLS, the IRLS variant that uses weight matrices derived from the arithmetic mean the
weights of IRLS-col and IRLS-row, cf. Lemma 2. The iterates X (n) for iteration n = 2000
of these algorithms exhibit relative errors of 0.240, 0.489 and 0.401, respectively, for the choice of
p = 0.1—furthermore, there is no choice of p that would lead to a convergence to X 0 .
To understand this very different behavior, we note that the n-th iterate of any of the four IRLS
variants can be written, using Appendix A, in a concise way as
X (n+1) = arg min hX vec ,W (n)X vec i,
Φ(X )=Y

11

(16)

(a) HM-IRLS

(b) IRLS-col

(c) IRLS-row

(d) AM-IRLS

Figure 2: Values of the matrix H (1) of ”weight coefficients” corresponding to the orthonormal basis
(ui(1)v j(1)∗ )i,4 j=1 after the first iteration in the example
where
4
Õ


(n) 2
(n)
hX vec ,W (n)X vec i = hX , U (n) H (n) ◦ (U (n)∗XV (n) ) V (n)∗ iF =
Hi(n)
j |hu i , Xv j i|

(17)

i, j=1

with X (n) = U (n) Σ(n)V (n)∗ =

Hi(n)
j

Í4

(n) (n) (n)
i=1 σi u i v i

being the SVD of X (n) , and


 2−p  −1
2−p


2 (σi(n) )2 + (ϵ (n) )2 ) 2 + (σ j(n) )2 + (ϵ (n) )2 ) 2



 (n) 2
 p−2

 (σ ) + (ϵ (n) )2 2

i
=
 p−2
(n) 2


) + (ϵ (n) )2 2
(σ

j

p−2
p−2 



 0.5 · (σ (n) )2 + (ϵ (n) )2  2 + (σ (n) )2 + (ϵ (n) )2  2
i
i


for HM-IRLS,
for IRLS-col,
for IRLS-row, and
for AM-IRLS,

for i, j ∈ {1, 2, 3, 4} and ϵ (n) = min(σ2(n) , ϵ (n−1) ).
The values of the matrix H (1) of weight coefficients after the first iteration in the above example
are visualized in Figure 2, for each of the four IRLS versions above.
The intuition for the superior behavior of HM-IRLS is now the following: Since large entries
of H (n) penalize the corresponding parts of the space Md1 ×d2 = span{ui(n)v j(n)∗ , i ∈ [d 1 ], j ∈ [d 2 ]}
in the minimization problem (16), large areas of blue and dark blue in Figure 2 indicate a benign
optimization landscape where the minimizer X (n+1) of (16) is able to improve considerably on the
previous iterate X (n) .
In particular, it can be seen that in the case of HM-IRLS, the penalties on the whole direct sum
of column and row space of the best rank-r approximation of X (n)


∗

T (n) := u 1(n) , . . . , ur(n) Z 1∗ + Z 2 v 1(n) , . . . , vr(n) : Z 1 ∈ Md1 ×r , Z 2 ∈ Md2 ×r ,
are small compared to the other penalites, since the coefficients of H (1) corresponding to T (1) are
exactly the ones in the first row and first column of the (4 × 4) matrices in Figure 2—a contrast that
becomes more and more pronounced as X (n) approaches the rank-r ground truth X 0 (with r = 1 in
the example).
On the other hand, IRLS-col, IRLS-row and AM-IRLS only have small coefficients on
smaller parts of T (n) , which, from a global perspective, explains why their usage might lead to nonglobal minima of the Schatten-p objective.
12

We note that the space T (n) plays also an important role in Riemannian optimization approaches
for matrix recovery problems [Van13], since it is also the tangent space of the smooth manifold of
rank-r matrices at the best rank-r approximation of X (n) .

4

Convergence results

In the following part, we state our main theoretical results about convergence properties of the
algorithm HM-IRLS. Furthermore, their relation to existing results for IRLS-col and IRLSrow is discussed.
It cannot be expected that a low-rank matrix recovery algorithm like HM-IRLS succeeds to
converge to a low-rank matrix without any assumptions on the measurement operator Φ that defines
the recovery problem (5). For the purpose of the convergence analysis of HM-IRLS, we introduce
the following strong Schatten-p null space property [FRW11, OMFH11, FR13].
Definition 4 (Strong Schatten-p null space property). Let 0 < p ≤ 1. We say that a linear map
Φ : Md1 ×d2 → Cm fulfills the strong Schatten-p null space property (Schatten-p NSP) of order r with
constant 0 < γr ≤ 1 if
Õ
 p/2

 Õ
r
d
γr
p
2
σi (X )
<
σi (X )
(18)
p
r 1− 2 i=r +1
i=1
for all X ∈ N (Φ) \ {0}.
Intuitively explained, if a map Φ fulfills the strong Schatten-p null space property of order r ,
there are no rank-r matrices in the null space and all the elements of the null space must not have
a quickly decaying spectrum.
Null space properties have already been used to guarantee the success of nuclear norm minimization (6), or Schatten-1 minimization in our terminology, for solving the low-rank matrix recovery problem [RXH11].
We note that the definitions of Schatten-p null space properties are quite analogous to the `p null space property in classical compressed sensing [FR13, Theorem 4.9], applied to the vector of
singular values. In particular, (18) implies that
r
Õ
i=1

p
σi (X )

<

d
Õ
i=r +1

p

σi (X )

for all X ∈ N (Φ) \ {0},

(19)

since kX kSp ≤ r 1/p−1/2 kX kS2 for X that is rank-r . This, in turn, ensures the existence of unique
solutions to (7) if Y = Φ(X 0 ) are the measurements of a low-rank matrix X 0 .
Proposition 5 ([Fou18]). Let Φ : Md1 ×d2 → Cm be a linear map, let 0 < p ≤ 1 and r ∈ N. Then every
matrix X 0 ∈ Md1 ×d2 such that rank(X 0 ) ≤ r and Φ(X 0 ) = Y ∈ Cm is the unique solution of Schatten-p
minimization (7) if and only if Φ fulfills (19).
Remark 6. The sufficiency of the Schatten-p NSP (19) in Proposition 5 already been pointed out by Oymak et al. [OMFH11]. The necessity as stated in the theorem, however, is due to a recent generalization
of Mirsky’s singular value inequalities to concave functions [Aud14, Fou18].

13

It can be seen that the (weak) Schatten-p NSP of (19) is a stronger property for larger p in the
sense that if 0 < p 0 ≤ p ≤ 1, the Schatten-p property implies the Schatten-p 0 property. Very related
to that, it can be seen that for any 0 < p ≤ 1, the strong Schatten-p null space property is implied
by a sufficiently small rank restricted isometry constant δr , which is a classical tool in the analysis of
low-rank matrix recovery algorithms [RFP10, CP11].
Definition 7 (Restricted isometry property (RIP)). The restricted isometry constant δr > 0 of order
r of the linear map Φ : Md1 ×d2 → Cm is defined as the smallest number such that
(1 − δr )kX kF2 ≤ kΦ(X )k`22 ≤ (1 + δr )kX kF2
for all matrices X ∈ Md1 ×d2 of rank at most r .
Indeed, it follows from the proof of [CDK15, Theorem 4.1] that a restricted isometry constant of
order 2r such that δ 2r < √ 2 ≈ 0.4531 implies the strong Schatten-p NSP of order r with a constant
2+3

√2
2+3
√
p
δ 2r
( 2+1)p
2p
(1−δ 2r )p .

γr < 1 for any 0 < p ≤ 1. More precisely, it can be seen that δ 2r <

implies that the strong

Schatten-p NSP (18) of order r holds with the constant γr =
Linear maps that are instances drawn from certain random models are known to fulfill the
restricted isometry property with high probability if the number of measurements is sufficiently
large [DR16], and, a fortiori, the Schatten-p null space property. In particular, this is true for (sub)Gaussian linear measurement maps Φ : Md1 ×d2 → Cm whose matrix representation is such that
1
Φ ∈ Cm×d1d2 , where e
Φ has i.i.d. standard (sub-)Gaussian entries,
√ e
m

(20)

as it is summarized in the following lemma.
Lemma 8. For any 0 < p ≤ 1, 0 < γ < 1 and any (sub-)Gaussian random operator Φ : Md1 ×d2 → Cm
(e.g. as defined in (20)), there exist constants C 1 > 1, C 2 > 0 such that if m ≥ C 1r (d 1 + d 2 ), the strong
Schatten-p null space property (18) of order r with constant γr < γ is fulfilled with probability at least
1 − e −C2m .

4.1

Local convergence for p < 1

In this section, we provide a convergence analysis for HM-IRLS covering several aspects. We are
able to show that the algorithm converges to stationary points of a smoothed Schatten-p functional
p
дϵ as in (21) without any additional assumptions on the measurement map Φ. Such guarantees have
already been obtained for IRLS algorithms with one-sided reweighting as IRLS-col and IRLSrow, in particular for p = 1 by Fornasier, Rauhut & Ward [FRW11] and for 0 < p ≤ 1 by Mohan &
Fazel [MF12].
Beyond that, assuming the measurement operator fulfills an appropriate Schatten-p null space
property as defined in Definition 4, we show the a-posteriori exact recovery statement that HMIRLS converges to the low-rank matrix X 0 if lim ϵn = 0, which only was shown for one-sided
n→∞

IRLS for the case p = 1 by [FRW11].
Moreover, we provide a local convergence guarantee stating that HM-IRLS recovers the lowrank matrix X 0 if we obtain an iterate X (n) that is close enough to X 0 , which is novel for IRLS
algorithms.
14

Let 0 < p ≤ 1 and ϵ > 0. To state the theorem, we introduce the ϵ-perturbed Schatten-p functional
p
дϵ : Md1 ×d2 → R ≥0 such that
d
Õ
p
p
дϵ (X ) =
(σi (X )2 + ϵ 2 ) 2
(21)
i=1

where σ (X ) ∈ Rd denotes the vector of singular values of X ∈ Md1 ×d2 .
Theorem 9. Let Φ : Md1 ×d2 → Cm be a linear operator and Y ∈ Ran(Φ) a vector in its range. Let
(X (n) )n ≥1 and (ϵ (n) )n ≥1 be the sequences produced by Algorithm 1 for input parameters Φ, Y , r and
0 < p ≤ 1, let ϵ = limn→∞ ϵ (n) .
(i) If ϵ = 0 and if Φ fulfills the strong Schatten-p NSP (18) of order r with constant 0 < γr < 1,
then the sequence (X (n) )n ≥1 converges to a matrix X ∈ Md1 ×d2 of rank at most r that is the
unique minimizer of the Schatten-p minimization problem (7). Moreover, there exists an absolute
constant Ĉ > 0 such that for any X with Φ(X ) = Y and any e
r ≤ r , it holds that
p

kX − X kF ≤
where Ĉ =

1−p/2

2p+1γ r
1−γ r

Ĉ

βer (X )Sp ,
r 1−p/2

and βer (X )Sp is the best rank-e
r Schatten-p approximation error of X , i.e.,

ek p , X
e ∈ Md1 ×d2 has rank e
βer (X )Sp := inf kX − X
r .
Sp

(22)

(ii) If ϵ > 0, then each accumulation point X of (X (n) )n ≥1 is a stationary point of the ϵ-perturbed
p
Schatten-p functional дϵ of (21) under the linear constraint Φ(X ) = Y . If additionally p = 1,
p
then X is the unique global minimizer of дϵ .
(iii) Assume that there exists a matrix X 0 ∈ Md1 ×d2 with Φ(X 0 ) = Y such that rank(X 0 ) = r ≤
min(d 1,d 2 )
, a constant 0 < ζ < 1 and an iteration n ∈ N such that
2
kX (n) − X 0 kS∞ ≤ ζ σer (X 0 )
and ϵ n = σr +1 (X n ). If Φ fulfills the strong Schatten-p NSP of order 2r with γ 2r < 1 and if the
0)
condition number κ = σσr1 (X
(X 0 ) of X 0 and ζ are sufficiently small (see condition (25) and formula
(26)), then
X (n) → X 0 for n → ∞.
It is important to note that by using Lemma 8, it follows that the assertions of Theorem 9(i) and
(iii) hold for (sub-)Gaussian operators (20) with high probability in the regime of measurements of
optimal sample complexity order. In particular, there exist constant oversampling factors ρ 1 , ρ 2 ≥ 1
such that the assertions of (i) and (iii) hold with high probability if m > ρ k r (d 1 + d 2 ), k ∈ {1, 2},
respectively.
Remark 10. However, if m < d 1d 2 , null space property-type assumptions as (18) or (19) do not hold for
the important case of matrix completion-type measurements [CR09], where Φ(X ) is given as m sample
entries
Φ(X )` = X i ` , j ` ,
` = 1, . . . , m,
(23)
15

and (i ` , j ` ) ∈ [d 1 ] × [d 2 ] for all ` ∈ [m], of the matrix X ∈ Md1 ×d2 , which also were considered in the
example of Section 3.
This means that parts (i) and (iii) of Theorem 9 do, unfortunately, not apply for matrix completion
measurements, which define a very relevant class of low-rank matrix recovery problems. This problem
is shared by any existing theory for IRLS algorithms for low-rank matrix recovery [FRW11, MF12].
However, in Section 5, we provide strong numerical evidence that HM-IRLS exhibits properties as
predicted by (i) and (iii) of Theorem 9 even for the matrix completion setting. We leave the extension
of the theory of HM-IRLS to matrix completion measurements as an open problem to be tackled by
techniques different from uniform null space properties [DR16, Section V].

4.2

Locally superlinear convergence rate for p < 1

Next, we state the second main theoretical result of this paper, Theorem 11. It shows that in a
neighborhood of a low-rank matrix X 0 that is compatible with the measurement vector Y , the algorithm HM-IRLS converges to X 0 with a convergence rate that is superlinear of the order 2 − p, if the
operator Φ fulfills an appropriate Schatten-p null space property.
Theorem 11 (Locally Superlinear Convergence Rate). Assume that the linear map Φ : Md1 ×d2 → Cm
fulfills the strong Schatten-p NSP of order 2r with constant γ 2r < 1 and that there exists a matrix
X 0 ∈ Md1 ×d2 with rank(X 0 ) = r ≤ min(d21,d2 ) such that Φ(X 0 ) = Y , let Φ, Y , r and 0 < p ≤ 1 be the input
0)
(n) := X (n) −X
parameters of Algorithm 1. Moreover, let κ = σσr1 (X
0
(X 0 ) be the condition number of X 0 and η
be the error matrices of the n-th output of Algorithm 1 for n ∈ N.
Assume that there exists an iteration n ∈ N and a constant 0 < ζ < 1 such that
kη (n) kS∞ ≤ ζ σr (X 0 )

(24)

and ϵ (n) = σr +1 (X (n) ). If additionally the condition number κ and ζ are small enough, or more precisely,
if
p(1−p)
<1
(25)
µ kη (n) kS∞
with the constant
 γ (3 + γ )(1 + γ )  2−p  d − r  2− p2 σ (X )p(p−1)
2r
2r
2r
r
0
µ := 2 (1 + γ 2r )
rp
κp
(1 − γ 2r )
r
(1 − ζ )2p
5p

then

p


 2−p
kη (n+1) kS∞ ≤ µ 1/p kη (n) kS∞

and

(26)


 2−p
kη (n+1) kSp ≤ µ 1/p kη (n) kSp

for all n ≥ n.
We think that the result of Theorem 11 is remarkable, since there are only few low-rank recovery
algorithms which exhibit either theoretically or practically verifiable superlinear convergence rates.
In particular, although the algorithms of [MMBS13] and NewtonSLRA of [SS16] do show superlinear convergence rates, the first are not competitive to HM-IRLS in terms of sample complexity
and the second has neither applicable theoretical guarantees for most of the interesting problems
nor the ability of solving medium size problems.

16

Remark 12. We observe that while the statement describes the observed rates of convergence very
accurately (cf. Section 5.2), the assumption (25) on the neighborhood that enables convergence of a
rate 2 − p is more pessimistic than our numerical experiments suggest. Our experiments confirm that
the local convergence rate of order 2 − p also holds for matrix completion measurements, where the
assumption of a Schatten-p null space property fails to hold, cf. Section 5.

4.3

Discussion and comparison with existing IRLS algorithms

Optimally, we would like to have a statement in Theorem 9 about the accumulation points X being
p
global minimizers of дϵ , instead of mere stationary points [FRW11, Theorem 6.11], [DDFG10, Theorem 5.3]. A statement that strong is, unfortunately, difficult to achieve due to the non-convexity
p
of the Schatten-p quasinorm and of the ϵ-perturbed version дϵ . Nevertheless, our theorems can
be seen as analogues of [DDFG10, Theorem 7.7], which discusses the convergence properties of an
IRLS algorithm for sparse recovery based on `p -minimization with p < 1.
As already mentioned in previous sections, Fornasier, Rauhut & Ward [FRW11] and Mohan &
Fazel [MF12] proposed IRLS algorithms for low-rank matrix recovery and analysed their convergence properties. The algorithm of [FRW11] corresponds (almost) to IRLS-col with p = 1 as
explained in Section 3. In this context, Theorem 9 recovers the results [FRW11, Theorem 6.11(i-ii)]
for p = 1 and generalize them, with weaker conclusions due to the non-convexity, to the cases
0 < p < 1. The algorithm IRLS-p of [MF12] is similar to the former, but differs in the choice
of the ϵ-smoothing and also covers non-convex choices 0 < p < 1. However, we note that in the
non-convex case, its convergence result [MF12, Theorem 5.1] corresponds to Theorem 9(ii), but does
not provide statements similar to (i) and (iii) of Theorem 9.
Theorem 11 with its analysis of the convergence rate is new in the sense that to the best of
our knowledge, there are no convergence rate proofs for IRLS algorithms for the low-rank matrix
recovery problem in the literature. Indeed, we refer to Remark 22 in Section 6.3 for an explanation
why the variants of [FRW11] and [MF12] cannot exhibit superlinear convergence rates, unlike HMIRLS.
We also note that there is a close connection between the statements of Theorems 9 and 11 and
results that were obtained by Daubechies, DeVore, Fornasier and Güntürk [DDFG10, Theorems 7.7
and 7.9] for an IRLS algorithm dedicated to the sparse vector recovery problem.

5

Numerical experiments

In this section, we demonstrate first that the superlinear convergence rate that was proven theoretically for Algorithm 1 (HM-IRLS) in Theorem 11 can indeed be accurately verified in numerical
experiments, even beyond measurement operators fulfilling the strong null space property, and
compare its performance to other variants of IRLS.
In Section 5.3, we then examine the recovery performance of HM-IRLS for the matrix completion setting with the performance of other state-of-the-art algorithms comparing the measurement
complexities that are needed for successful recovery for many random instances.
The numerical experiments are conducted on Linux and Mac systems with MATLAB R2017b. An
implementation of the HM-IRLS algorithm and a minimal test example are available at https://wwwm15.ma.tum.de/Allgemeines/SoftwareSite.

17

5.1

Experimental setup

In the experiments, we sample (d 1 × d 2 ) dimensional ground truth matrices X 0 of rank r such that
X 0 = U ΣV ∗ , where U ∈ Rd1 ×r and V ∈ Rd2 ×r are independent matrices with i.i.d. standard Gaussian
entries and Σ ∈ Rr ×r is a diagonal matrix with i.i.d. standard Gaussian diagonal entries, independent
from U and V .
We recall that a rank-r matrix X ∈ Md1 ×d2 has d f = r (d 1 +d 2 −r ) degrees of freedom, which is the
theoretical lower bound on the number of measurements that are necessary for exact reconstruction
[CP11]. The random measurement setting we use in the experiments can be described as follows:
We take measurements of matrix completion type, sampling m = bρd f c entries of X 0 uniformly
over its d 1d 2 indices to obtain Y = Φ(X 0 ). Here, ρ is such that dd1df 2 ≥ ρ ≥ 1 and parametrizes the
difficulty of the reconstruction problem, from very hard problems for ρ ≈ 1 to easier problems for
larger ρ.
However, this uniform sampling of Φ could yield instances of measurement operators whose
information content is not large enough to ensure well-posedness of the corresponding low-rank
matrix recovery problem, even if ρ > 1. More precisely, it is impossible to recover a matrix exactly
if the number of revealed entries in any row or column is smaller than its rank r , which is explained
and shown in the context of the proof of [PABN16, Theorem 1].
Thus, in order to provide for a sensible measurement model for small ρ, we exclude operators Φ
that sample fewer than r entries in any row or column. Therefore, we adapt the uniform sampling
model such that operators Φ are discarded and sampled again until the requirement of at least r
entries per column and row is met and recovery can be achieved from a theoretical point of view.
We note that the described phenomenon is very related to the fact that matrix completion recovery guarantees for the uniform sampling model require at least one additional log factor, i.e.,
they require at least m ≥ log(max(d 1 , d 2 ))d f sampled entries [DR16, Section V].
While we detail the experiments for the matrix completion measurement setting just described
in the remaining section, we add that Gaussian measurement models also lead to very similar results
in experiments.

5.2

Convergence rate comparison with other IRLS algorithms

In this subsection, we vary the Schatten-p parameter between 0 and 1 and compare the corresponding convergence behavior of HM-IRLS with the IRLS variant IRLS-col, which performs the
reweighting just in the column space, and with the arithmetic mean variant AM-IRLS. The latter
two coincide with Algorithm 1 except that the weight matrices are chosen as described in Equation (17) in Section 3.
We note that IRLS-col is very similar to the IRLS algorithms of [FRW11] and [MF12] and
differs from them basically just in the choice of the ϵ-smoothing. We present the experiments with
IRLS-col to isolate the influence of the weight matrix type, but very similar results can be observed for the algorithms of [FRW11] and [MF12].2
In the matrix completion setup of Section 5.1, we choose d 1 = d 2 = 40, r = 10 and distinguish
easy, hard and very hard problems corresponding to oversampling factors ρ of 2.0, 1.2 and 1.0,
respectively. The algorithms are provided with the ground truth rank r and are stopped whenever
2 Implementations of the mentioned authors’ algorithms were downloaded from https://faculty.
washington.edu/mfazel/ and https://github.com/rward314/IRLSM, respectively.

18

100

10-5

p=1, HM-IRLS
p=1, AM-IRLS
p=1.0, IRLS-col
p=0.8, HM-IRLS
p=0.8, AM-IRLS
p=0.8, IRLS-col
p=0.5, HM-IRLS
p=0.5, AM-IRLS
p=0.5, IRLS-col
p=0.1, HM-IRLS
p=0.1, AM-IRLS
p=0.1, IRLS-col

10-10

10-15

0

5

10

15

20

25

30

35

40

45

50

Iteration n

Figure 3: Relative Frobenius errors as a function of the iteration n for oversampling factor ρ = 2
(easy problem).
the relative change of Frobenius norm kX (n) − X (n−1) kF /kX (n−1) kF drops below the threshold of
10−10 or a maximal iteration of iterations n max is reached.
5.2.1

Convergence rates

First, we study the behavior of the three IRLS algorithms for the easy setting of an oversampling
2 −r )
factor of ρ = 2, which means that 2r (dd1 1+d
= 0.875 of the entries are sampled, and parameters
d2
p ∈ {0.1, 0.5, 0.8, 1}.
In Figure 3, we observe that for p = 1, HM-IRLS, AM-IRLS and IRLS-col have a quite
similar behavior, as the relative Frobenius errors kX (n) − X 0 kF /kX 0 kF decrease only slowly, i.e.,
even a linear rate is hardly identifiable. For choices p < 1 that correspond to non-convex objectives,
we observe a very fast, superlinear convergence of HM-IRLS, as the iterates X (n) converge up to
a relative error of less than 10−12 within fewer than 20 iterations for p ∈ {0.8, 0.5, 0.1}. Precise
calculations verify that the rate of convergences are indeed of order 2 − p, the order predicted by
Theorem 11. We note that this fast convergence rate kicks in not only locally, but starting from the
very first iteration.
On the other hand, it is easy to see that AM-IRLS and IRLS-col converge linearly, but
not superlinearly to the ground truth X 0 for p ∈ {0.8, 0.5, 0.1}. The linear rate of AM-IRLS is
slightly better than the one of IRLS-col, but the numerical stability of AM-IRLS deteriorates
for p = 0.1 close to the ground truth (after iteration 43). This is due to a bad conditioning of the
quadratic problems as the X (n) are close to rank-r matrices. In contrast, no numerical instability
issues can be observed for HM-IRLS.
For the hard matrix completion problems with oversampling factor of ρ = 1.2, we observe that
for p = 0.8, the three algorithms typically do not converge to ground truth. This can be seen in
19

10

10

0

p=0.8, HM-IRLS

p=0.5, HM-IRLS

p=0.25, HM-IRLS

p=0.01, HM-IRLS

p=0.8, AM-IRLS

p=0.5, AM-IRLS

p=0.25, AM-IRLS

p=0.01, AM-IRLS

p=0.8, IRLS-col

p=0.5, IRLS-col

p=0.25, IRLS-col

p=0.01, IRLS-col

-1

10 -2

0

10

20

30

40

50

60

70

80

90

100

Iteration n

Figure 4: Relative Frobenius errors as a function of the iteration n for oversampling factor ρ = 1.2
(hard problem). Left column: y-range [10−10 ; 100 ]. Right column: Enlarged section of left column
corresponding to y-range of [10−2 ; 100 ].
the example that is shown in Figure 4, where HM-IRLS, AM-IRLS and IRLS-col all exhibit
a relative error of 0.27 after 100 iterations. We do not visualize the result for p = 1, as the iterates
of the three algorithms do not converge to the ground truth either, which is to be expected: In
some sense, they implement nuclear norm minimization, which is typically not able to recover a
low-rank matrix from measurements with an oversampling factor as small as ρ = 1.2 [DGM13].
The dramatically different behavior between HM-IRLS and the other approaches becomes very
apparent for more non-convex choices of p ∈ {0.01, 0.25, 0.5}, where the former converges up to a
relative Frobenius error of less than 10−10 within 15 to 35 iterations, while the others do not reach a
relative error of 10−2 even after 100 iterations. For HM-IRLS, the convergence of order 2 −p can be
very well locally observed also here, it just takes some iterations until the superlinear convergence
begins, which is due to the increased difficulty of the recovery problem.
Finally, we see in the example shown in Figure 5 that even for the very hard problems where
ρ = 1, which means that the number of sampled entries corresponds exactly to the degrees of
freedom r (d 1 + d 2 − r ), HM-IRLS can be successful to recover the rank-r matrix if the parameter
p is chosen small enough (here: p ≤ 0.25). This is not the case for the algorithms AM-IRLS and
IRLS-col.
5.2.2 HM-IRLS as the best extension of IRLS for sparse recovery
We summarize that among the three variants HM-IRLS, AM-IRLS and IRLS-col, only HMIRLS is able to solve the low-rank matrix recovery problem for very low sample complexities
corresponding to ρ ≈ 1. Furthermore, it is the only IRLS algorithm for low-rank matrix recovery
that exhibits a superlinear rate of convergence at all.
It is worthwhile to compare the properties of HM-IRLS with the behavior of the IRLS algorithm
of [DDFG10] designed to solve the sparse vector recovery problem by mimicking `p -minimization
for 0 < p ≤ 1. While neither IRLS-col nor AM-IRLS are able to generalize the superlinear
convergence behavior of [DDFG10] (which is illustrated in Figure 8.3 of the same paper) to the
20

10 0
10 0

10

p=0.5, HM-IRLS

p=0.25, HM-IRLS

p=0.01, HM-IRLS

p=0.8, AM-IRLS

p=0.5, AM-IRLS

p=0.25, AM-IRLS

p=0.01, AM-IRLS

p=0.8, IRLS-col

p=0.5, IRLS-col

p=0.25, IRLS-col

p=0.01, IRLS-col

-2

p=0.8, HM-IRLS
p=0.8, AM-IRLS
p=0.8, IRLS-col
p=0.5, HM-IRLS
p=0.5, AM-IRLS
p=0.5, IRLS-col
p=0.25, HM-IRLS
p=0.25, AM-IRLS
p=0.25, IRLS-col
p=0.01, HM-IRLS
p=0.01, AM-IRLS
p=0.01, IRLS-col

10 -4

10 -6

10

p=0.8, HM-IRLS

10 -1

-8

10 -2

10 -10

0

10

20

30

40

50

60

70

80

90

100

0

10

20

30

40

50

60

70

80

90

100

Iteration n

Iteration n

Figure 5: Relative Frobenius errors as a function of the iteration n for oversampling factor ρ = 1.0
(very hard problem). Left column: y-range [10−10 ; 100 ]. Right column: Enlarged section of left
column corresponding to y-range of [10−2 ; 100 ].
low-rank matrix recovery problem, HM-IRLS is, as can be seen in Figures 3 to 5.
Taking the theoretical guarantees as well as the numerical evidence into account, we claim that
HM-IRLS is the presently best extension of IRLS for vector recovery [DDFG10] to the low-rank matrix
recovery setting, providing a substantial improvement over the reweighting strategies of [FRW11]
and [MF12].
Moreover, we mention two observations which suggest that HM-IRLS has in some sense even
more favorable properties than the algorithm of [DDFG10]: First, the discussion of [DDFG10, Section 8] states that a superlinear convergence can only be observed locally after a considerable
amount of iterations with just a linear error decay. In contrast to that, HM-IRLS exhibits a superlinear error decay quite early (i.e., for example as early as after two iterations), at least if the
sample complexity is large enough, cf. Figure 3.
Secondly, it can be observed that the convergence of the algorithm of [DDFG10] to a sparse
vector often breaks down if p is smaller than 0.5 [DDFG10, Section 8]. In contrast to that, we
observe that HM-IRLS does not suffer from this loss of global convergence for p  0.5. Thus, a
choice of very small parameters p ≈ 0.1 or smaller is suggested as such a choice is accompanied by
a very fast convergence.

5.3

Recovery performance compared to state-of-the-art algorithms

After comparing the performance of HM-IRLS with other IRLS variants, we now conduct experiments to compare the empirical performance of HM-IRLS also to that of low-rank matrix recovery
algorithms different from IRLS.
To obtain a comprehensive picture, we consider not only the IRLS variants AM-IRLS and
IRLS-col, but a variety of state-of-the-art methods in the experiments, as Riemannian optimization technique Riemann Opt [Van13], the alternating minimization approaches AltMin
[HH09], ASD [TW16] and BFGD [PKCS16], and finally the algorithms Matrix ALPS II [KC14]
21

1

Empirical probability of successful recovery

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

Riemann_Opt
0.1

0
1.0

1.2

1.5

1.7

2.0

2.3

2.6

Figure 6: Comparison of empirical success rates of state-of-the-art algorithms, as a function of the
oversampling factor ρ
and CGIHT Matrix [BTW15], which are based on iterative hard thresholding. As the IRLS variants we consider, all these algorthms use knowledge about the true ground truth rank r .
In the experiments, we examine the empirical recovery probabilities of the different algorithms
systematically for varying oversampling factors ρ, determining the difficulty of the low-rank recovery problem as the sample complexity fulfills m = bρd f c. We recall that a large parameter ρ
corresponds to an easy reconstruction problem, while a small ρ, e.g., ρ ≈ 1, defines a very hard
problem.
We choose d 1 = d 2 = 100 and the r = 8 as parameter of the experimental setting, conducting
the experiments to recover rank-8 matrices X 0 ∈ R100×100 . We remain in the matrix completion
measurement setting described in Section 5.1, but sample now 150 random instances of X 0 and Φ
for different numbers of measurements varying between m min = 1500 to m max = 4000. This means
that the oversampling factor ρ increases from ρ min = 0.975 to ρ max = 2.60. For each algorithm, a
successful recovery of X 0 is defined as a relative Frobenius error kX out − X 0 kF /kX 0 kF of the matrix
X out returned by the algorithm of smaller than 10−3 . The algorithms are run until stagnation of
the iterates or until the maximal number of iterations n max = 3000 is reached. The number n max is
chosen large enough to ensure that a recovery failure is not due to a lack of iterations.
In the experiments, except for AltMin, for which we used our own implementation, we used
implementations provided by the authors of the corresponding papers for the respective algorithms,
using default input parameters provided by the authors. The respective code sources can be found
in the references.

22

5.3.1

Beyond the state-of-the-art performance of HM-IRLS

The results of the experiment can be seen in Figure 6. We observe that HM-IRLS exhibits a very
high empirical recovery probability for p = 0.1 and p = 0.5 as soon as the sample complexity
parameter ρ is slightly larger than 1.0, which means that m = bρr (d 1 + d 2 − r )c measurements
suffice to recover (d 1 × d 2 )-dimensional rank-r matrices with ρ close to 1. This is very close to the
information theoretical lower bound of d f = r (d 1 + d 2 − r ). Very interestingly, it can be observed
that the empirical recovery percentage reaches almost 100% already for an oversampling factor of
ρ ≈ 1.1, and remains at exactly 100% starting from ρ ≈ 1.2.
Quite good success rates can also be observed for the algorithms AM-IRLS and IRLS-col
for non-convex parameter choices p ∈ {0.1, 0.5}, reaching an empirical success probability of almost
100% at around ρ = 1.5. AM-IRLS performs only marginally better than the classical IRLS strategy
IRLS-col, which are both outperformed considerably by HM-IRLS. It is important to note that
in accordance to what was observed in Section 5.2, in the successful instances, the error threshold
that defines successful recovery is achieved already after a few dozen iterations for HM-IRLS, while
typically only after several or many hundreds for AM-IRLS and IRLS-col. Furthermore, it is
interesting to observe that the algorithm IRLS-MF, which corresponds to the variant studied and
implemented by [MF12] and differs from IRLS-col mainly only in the choice of the ϵ-smoothing
(14), has a considerably worse performance than the other IRLS methods. This is plausible since the
smoothing influences severely the optimization landscape of the objective to be minimized.
The strong performance of HM-IRLS is in stark contrast to the behavior of all the algorithms
that are based on different approaches than IRLS and that we considered in our experiments. They
basically never recover any rank-r matrix if ρ < 1.2, and most of the algorithms need a sample
complexity parameter of ρ > 1.7 to exceed a empirical recovery probability of a mere 50%. A
success rate of close to 80% is reached not before raising ρ above 2.0 in our experimental setting,
and also only for a subset of the comparison algorithms, in particular for Matrix ALPS II,
ASD, AltMin. The empirical probability of 100% is only reached for some of the IRLS methods,
and not for any competing method in our experimental setting, even for quite large oversampling
factors such as ρ = 2.5. While we do not rule out that a possible parameter tuning could improve
the performance of any of the algorithms slightly, we conclude that for hard matrix completion
problems, the experimental evidence for the vast differences in the recovery performance of HMIRLS compared to other methods is very apparent.
Thus, our observation is that the proposed HM-IRLS algorithm recovers low-rank matrices systematically with nearly the optimal number of measurements and needs fewer measurements than all
the state-of-the-art algorithms we included in our experiments, if the non-convexity parameter p is
chosen such that p  1.
We also note that the very sharp phase transition between failure and success that can be observed in Figure 6 for HM-IRLS indicates that the sample complexity parameter ρ is indeed the
major variable determining the success of HM-IRLS. In contrast, the wider phase transitions for
the other algorithms suggest that they might depend more on other factors, as the realizations of
the random sampling model and the interplay of measurement operator Φ and ground truth matrix
X0.
Another conclusion that can be drawn from the empirical recovery probability of 1 is that, despite the severe non-convexity of the underlying Schatten-p quasinorm for, e.g., p = 0.1, HM-IRLS
with the initialization of X (1) as the Frobenius norm minimizer does not get stuck in stationary
points if the oversampling factor is large enough. Further experiments conducted with random
23

initializations as well as severely adversary initializations, e.g., with starting points chosen in the
orthogonal complement of the spaces spanned by the singular vectors of the ground truth matrix
X 0 , lead to comparable results. Therefore, we claim that HM-IRLS exhibits a global convergence
behavior in interesting application cases and for oversampling factor ranges for which competing
non-convex low-rank matrix recovery algorithms fail to succeed. We consider a theoretical investigation of such behavior as an interesting open problem to explore.

5.4

Computational complexity

e (n) , cf. (15), is an inverse of a (d 1d 2 × d 1d 2 )-matrix
While the harmonic mean weight matrix W
and therefore in general a dense (d 1d 2 × d 1d 2 )-matrix, it is important to note that it never has to
be computed explicitly
implementation of HM-IRLS;neither is it necessary to compute its
 in an(n)
1
(n)
−1
(n)
e ) = U (Σ )2−p U (n)∗ ⊕ V (n) (Σ(n) )2−p V (n)∗ explicitly.
inverse (W
2
Indeed, as it can be seen in (13) and by the definition of the Kronecker sum (55), the harmonic
mean weight matrix appears just as the linear operator (W (n) )−1 on the space of matrices Md1 ×d2 ,
whose action consists of a left- and right-sided matrix multiplication, cf. (12). Therefore, the application of (W (n) )−1 is O(d 1d 2 (d 1 + d 2 )) by the naive matrix multiplication algorithm, and can be
easily parallelized.
While this useful observation is helpful for the implementation of HM-IRLS, it is not true for
(n)
AM-IRLS, as the action of (W(arith)
)−1 , the inverse of the arithmetic mean weight matrix at iteration
n, is not representable as a sum of left- and right-sided matrix multiplication. This means that even
the execution of a fixed number of iterations of HM-IRLS is faster than computational advantage
over AM-IRLS.
f(n)−1 ◦ Φ∗ ∈ Mm×m depends on the linear measurement operator Φ.
The cost to compute Φ ◦ W
In the matrix completion setting (23), no additional arithmetic operations have to be performed, as
f(n)−1 ◦ Φ∗ is a
Φ is a just a selection operator in this case, and for HM-IRLS, this means that Φ ◦ W
sparse matrix.
Thus, the algorithm HM-IRLS consists of basically of two computational steps per iteration:
The computation of the SVD of the d 1 × d 2 -matrix X (n) and the solution of the linearly constrained
least squares problem in (13). The first is of time complexity O(d 1d 2 min(d 1 , d 2 )). The time complexity
of the second depends on Φ, but is dominated by the inversion of a symmetric, m × m sparse linear
system in the matrix completion setting, if m is the number of given entries. This has a worst case
time complexity of O(max(d 1 , d 2 )3r 3 ) if ρ is just a constant oversampling factor.
For the matrix completion case, this allows us to recover low-rank matrices up to, e.g., d 1 = d 2 =
3000 on a single machine given very few entries with HM-IRLS.
Acceleration possibilities and extensions
To tackle higher dimensionalities in reasonable runtimes, a key strategy could be to address the
computational bottleneck of HM-IRLS, the solution of the m × m linear system in (13), by using
iterative methods. For IRLS algorithms designed for the related sparse recovery problem, the usage
of conjugate gradient (CG) methods is discussed in [FPRW16]. By coupling the accuracy of the CG
solutions to the outer IRLS iteration and using appropriate preconditioning, the authors obtain a
competitive solver for the sparse recovery problem, also providing a convergence analysis. Similar
ideas could be used for an acceleration of HM-IRLS.

24

It is interesting to see if further computational improvements can be achieved by combining the
ideas of HM-IRLS with the usage of truncated and randomized SVDs [HMT11], replacing the full
SVDs of the X (n) that are needed to define the linear operator (W (n) )−1 in Algorithm 1.

6

Theoretical analysis

For the theoretical analysis of HM-IRLS, we introduce the following auxiliary functional Jp , leading to a variational interpretation of the algorithm. In the whole section, we denote d = min(d 1 , d 2 )
and D = max(d 1 , d 2 ).
Definition 13. Let 0 < p ≤ 1. Given a full rank matrix Z ∈ Md1 ×d2 , let
i −1 


h
1
e (Z ) := 2 Id2 ⊗ (ZZ ∗ ) 21 (ZZ ∗ ) 12 ⊕ (Z ∗Z ) 12
W
(Z ∗Z ) 2 ⊗ Id1 ∈ Hd1d2 ×d1d2
e associated to Z .
be the harmonic mean matrix W
We define the auxiliary functional Jp : Md1 ×d2 × R ≥0 ×Md1 ×d2 → R ≥0 as

Jp (X , ϵ, Z ) :=



 p kX vec k 2

2

f(Z ))
`2 (W

+

d
Í

ϵ 2p
2

i=1

σi (Z ) +

2−p
2

d
Í
i=1

p

σi (Z ) (p−2)


 +∞


if rank(Z ) = d,
if rank(Z ) < d.

e of Definition 13 is just the harmonic mean of the matrices W
e1 :=
We note that the matrix W
1
1
1
1
∗
∗
∗
∗
e
2
2
2
2
Id2 ⊗ (ZZ ) and W2 = (Z Z ) ⊗ Id1 , as introduced in Section 2.3, if (ZZ ) and (Z Z ) are positive
1
1
e1 + W
e2 is invertible and as (A−1 + B −1 )−1 =
definite. Indeed, in this case, (ZZ ∗ ) 2 ⊕ (Z ∗Z ) 2 = W
A(A + B)−1 B for any positive definite matrices A and B of the same dimensions,

e (Z ) = 2W
e1 W
e1 + W
e2 −1W
e2 = 2(W
e1−1 + W
e2−1 )−1 .
W
(27)
e (Z ) as it is well-defined for any full-rank Z ∈ Md1 ×d2 and as
We use the more general definition W
1
it allows to handle the case of non-square matrices, i.e., the case d 1 ,d 2 , as in this case (ZZ ∗ ) 2 or
1
e + and W
e + of the matrices W
e1
(Z ∗Z ) 2 has to be singular. Using the Moore-Penrose pseudo inverse W
1
2
e2 , we can rewrite W
e (Z ) from Definition 13 as
and W

e (Z ) = 2W
e1 W
e1 + W
e2 −1W
e2 = 2(W
e1+ + W
e2+ )−1 .
W
With the auxiliary functional Jp at hand, we can interpret Algorithm 1 as an alternating minimization of the functional Jp (X , ϵ, Z ) with respect to its arguments X , ϵ and Z .
e (n+1) as the evaluation
In the following,
 we derive the formula (15) for the weight matrix W
(n+1)
(n+1)
e
e
e
W
=W Z
of W from Definition 13 at the minimizer
Z (n+1) = arg min Jp (X (n+1) , ϵ (n+1) , Z ),
Z ∈Md1 ×d2

(28)

with the minimizer being unique. Similarly, the formula (13) can be interpreted as
(n)
(n)
X (n+1) = arg min kX vec k`2 (W
f(Z (n) )) = arg min Jp (X , ϵ , Z )
X ∈Md1 ×d2
Φ(X )=Y

2

X ∈Md1 ×d2
Φ(X )=Y

(29)

These observations constitute the starting point of the convergence analysis of Algorithm 1,
which is detailed subsequently after the verification of the optimization steps.
25

6.1

Optimization of Jp with respect to Z and X

Í
We fix X ∈ Md1 ×d2 with singular value decomposition X = di=1 σi ui vi∗ , where ui ∈ Cd1 vi ∈ Cd2 are
the left and right singular vectors respectively and σi = σi (X ) denote its singular values for i ∈ [d].
Our objective in the following is the justification of formula (15). To yield the building blocks
e (n+1) , we consider the minimization problem
of the weight matrix W
arg min Jp (X , ϵ, Z )

(30)

Z ∈Md1 ×d2

for ϵ > 0.
Lemma 14. The unique minimizer of (30) is given by
Z opt =

d
Õ

(σi (X )2 + ϵ 2 )

p−2
2

ui vi∗ .

i=1

Furthermore, the value of Jp at the minimizer Z opt is
Jp (X , ϵ, Z opt ) =

d
Õ
p
p
(σi (X )2 + ϵ 2 ) 2 =: дϵ (X )

(31)

i=1

for p > 0.
The proof of this results is detailed in the appendix.
Remark 15. We note that the value of Jp (X , ϵ, Z opt ) can be interpreted as a smooth ϵ-perturbation of
a p-th power of a Schatten-p quasi-norm of the matrix X . In fact, for ϵ = 0 we have
p

p

Jp (X , 0, Z opt ) = kX kSp = д0 (X ).
Now, we show that our definition rule (13) of X (n+1) in Algorithm 1 can be interpreted as a
minimization of the auxiliary functional Jp with respect to the variable X . Additionally, this minimization step can be formulated as the solution of a weighted least squares problem with weight
e (n) . This is summarized in the following lemma.
matrix W
e (Z ) := 2([(ZZ ∗ ) 21 ]+ ⊕
Lemma 16. Let 0 < p ≤ 1. Given a full-rank matrix Z ∈ Md1 ×d2 , let W
1
[(Z ∗Z ) 2 ]+ )−1 ∈ Hd1d2 ×d1d2 be the matrix from Definition 13 and W −1 : Md1 ×d2 → Md1 ×d2 the linear
operator of its inverse
i
1
1
1h
W −1 (X ) :=
[(ZZ ∗ ) 2 ]+X + X [(Z ∗Z ) 2 ]+ .
2
Then the matrix
 
X opt = W −1 ◦ Φ∗ ◦ (Φ ◦ W −1 ◦ Φ∗ )−1 Y ∈ Md1 ×d2
is the unique minimizer of the optimization problems
arg min Jp (X , ϵ, Z ) = arg min kX vec k`2 (W
f) .
Φ(X )=Y

Φ(X )=Y

2

(32)

Moreover, a matrix X opt ∈ Md1 ×d2 is a minimizer of the minimization problem (32) if and only if it
fulfills the property
e (Z )(X opt )vec , H vec i`2 = 0 for all H ∈ N (Φ) and Φ(X opt ) = Y .
hW
In Appendix B.3, the interested reader can find a sketch of the proof of this lemma.
26

(33)

6.2

Basic properties of the algorithm and convergence results

In the following subsection, we will have a closer look at Algorithm 1 and point out some of its
properties, in particular, the boundedness of the iterates (X (n) )n ∈N and the fact that two consecutive
iterates are getting arbitrarily close as n → ∞. These results will be useful to develop finally the
proof of convergence and to determine the rate of convergence of Algorithm 1 under conditions
determined along the way.
Lemma 17. Let (X (n) , ϵ (n) )n ∈N be the sequence of iterates and smoothing parameters of Algorithm 1.
Í
Let X (n) = di=1 σi(n)ui(n)vi(n)∗ be the SVD of the n-th iterate X (n) . Let (Z (n) )n ∈N be a corresponding
sequence such that
d
Õ
p−2
Z (n) =
(σi(n)2 + ϵ (n)2 ) 2 ui(n)vi(n)∗
i=1

for n ∈ N. Then the following properties hold:
(a) Jp (X (n) , ϵ (n) , Z (n) ) ≥ Jp (X (n+1) , ϵ (n+1) , Z (n+1) ) for all n ≥ 1,
p

(b) kX (n) kSp ≤ Jp (X (1) , ϵ (0) , Z (0) ) =: Jp,0 for all n ≥ 1,
(c) The iterates X (n) , X (n+1) come arbitrarily close as n → ∞, i.e.,
lim k(X (n) − X (n+1) )vec k`22 = 0.
n→∞

At this point we notice that, assuming X (n) → X and ϵ (n) → ϵ for n → ∞ with the limit point
(X , ϵ) ∈ Md1 ×d2 × R ≥0 , it would follow that
p

Jp (X (n) , ϵ (n) , Z (n) ) → дϵ (X )
for n → ∞ by equation (31).
Now, let ϵ > 0, a measurement vector Y ∈ Cm and the linear operator Φ be given and consider
the optimization problem
p
min дϵ (X )
(34)
X ∈Md1 ×d2
Φ(X )=Y

p
Í
p
p
with дϵ (X ) = di=1 (σi (X )2 + ϵ 2 ) 2 and σi (X ) being the i-th singular value of X , cf. (31). If дϵ (X ) is
non-convex, which is the case for p < 1, one might practically only be able to find critical points of
the problem.

Lemma 18. Let X ∈ Md1 ×d2 be a matrix with the SVD such that X =
Íd
∗
i=1 σi u i v i , let ϵ > 0. If we define
e (X , ϵ) = 2
W

 Õ
d
d
 Õ
  −1
2−p
2−p
2
2 2
∗
2
2 2
∗
(σi + ϵ ) ui ui ⊕
(σi + ϵ ) vi vi
∈ Hd1d2 ×d1d2 ,
i=1

i=1

e (X (n) , ϵ (n) ) = W
e (n) , with W
e (n) defined as in Algorithm 1, cf. (10).
then W
Furthermore, X is a critical point of the optimization problem (34) if and only if
e (X , ϵ)X vec , H vec i`2 = 0 for all H ∈ N (Φ) and Φ(X ) = Y .
hW
p

In the case that дϵ is convex, i.e., if p = 1, (35) implies that X is the unique minimizer of (34).
27

(35)

Now, we have some basic properties of the algorithm at hand that allow us, together with the
strong nullspace property in Definition 4 to carry out the proof of the convergence result in Theorem 9. The proof is sketched in Appendix C using the results above.

6.3

Locally superlinear convergence

In the proof of Theorem 11 we use the following bound on perturbations of the singular value
decomposition, which is originally due to Wedin [Wed72]. It bounds the alignment of the subspaces
spanned by the singular vectors of two matrices by their norm distance, given a gap between the first
singular values of the one matrix and the last singular values of the other matrix that is sufficiently
pronounced.
Lemma 19 (Wedin’s bound [Ste06]). Let X and X̄ be two matrices of the same size and their singular
value decompositions
 


 
 Σ1 0 V1∗
 Σ̄1 0 V̄1∗
Ū
Ū
and
D̄
=
,
X = U1 U2
1
2
0 Σ2 V2∗
0 Σ̄2 V̄2∗
where the submatrices have the sizes of corresponding dimensions. Suppose that δ, α satisfying 0 <
δ ≤ α are such that α ≤ σmin (Σ1 ) and σmax (Σ̄2 ) < α − δ . Then
kŪ2∗U1 kS∞ ≤

√ kX − X̄ kS∞
√ kX − X̄ kS∞
and kV̄2∗V1 kS∞ ≤ 2
.
2
δ
δ

(36)

As a first step towards the proof of Theorem 11, we show the following lemma.
Lemma 20. Let (X (n) )n be the output sequence of Algorithm 1 for parameters Φ, Y , r and 0 < p ≤ 1,
and X 0 ∈ Md1 ×d2 be a matrix such that Φ(X 0 ) = Y .
(n+1)
(i) Let η 2r
be the best rank-2r approximation of η (n+1) = X (n+1) − X 0 . Then

kη

(n+1)

(n+1) 2p
− η 2r
kSp

≤2

2−p

 Õ
d
i=r +1

σi2 (X (n) )

+ϵ

(n)2 

p
2

 2−p

(n+1)
kη vec
k

2p
,
f(n) )
`2 (W

e (n) denotes the harmonic mean weight matrix from (10).
where W
(ii) Assume that the linear map Φ : Md1 ×d2 → Cm fulfills the strong Schatten-p NSP of order 2r with
constant γ 2r < 1. Then
2p

kη (n+1) kS2 ≤ 2p

2−p

γ 2r
r 2−p

 Õ
d
i=r +1

σi2 (X (n) ) + ϵ (n)2

 p2

 2−p

(n+1)
kη vec
k

2p
f(n) )
`2 (W

(iii) Under the same assumption as for (ii), it holds that
2p

kη (n+1) kSp ≤ (1 + γ 2r )2 22−p

 Õ
d
i=r +1

σi2 (X (n) ) + ϵ (n)2

28

 p2

 2−p

(n+1)
kη vec
k

2p
.
f(n) )
`2 (W

(37)

Proof of Lemma 20. (i) Let the X (n) = Ue(n) Σ(n)Ve(n)∗ be the (full) singular value decomposition of
e(n) ∈ Ud1 and Ve(n) ∈ Ud2 are unitary matrices and Σ(n) = diag(σ1 (X (n) ), . . . , σr (X (n) )) ∈
X (n) , i.e., U
Md1 ×d2 . We define UT(n) ∈ Md1 ×r as the matrix of the first r columns of Ue(n) and UT(n)
∈ Md1 ×(d1 −r ) as
c


(n)
e(n) = U (n) U (n) , and similarly V and V (n) .
the matrix of its last d 1 − r columns, so that U
T
Tc
T
Tc
As Id1 = UT(n)UT(n)∗ + UT(n)
UT(n)∗
and Id2 = VT(n)VT(n)∗ + VT(n)
VT(n)∗
, we note that
c
c
c
c
UT(n)
UT(n)∗
η (n+1)VT(n)
VT(n)∗
= η (n+1) − UT(n)UT(n)∗η (n+1) + UT(n)
UT(n)∗
η (n+1)VT(n)VT(n)∗ ,
c
c
c
c
c
c
while UT(n)UT(n)∗η (n+1) + UT(n)
UT(n)∗
η (n+1)VT(n)VT(n)∗ has a rank of at most 2r . This implies that
c
c
(n+1)
kη (n+1) − η 2r
kSp ≤ kUT(n)
UT(n)∗
η (n+1)VT(n)
VT(n)∗
kSp = kUT(n)∗
η (n+1)VT(n)
kSp .
c
c
c
c
c
c

(38)

Using the definitions of Ue(n) and Ve(n) , we write the harmonic mean weight matrices of the n-th
iteration (10) as

e (n) = 2(Ve(n) ⊗ Ue(n) ) Σd(n)2−p ⊕ Σd(n)2−p −1 (Ve(n) ⊗ Ue(n) )∗ ,
(39)
W
1
2
(n)

(n)

where Σd1 ∈ Md1 ×d1 and Σd2 ∈ Md2 ×d2 are the diagonal matrices with the smoothed singular values
of X (n) from (11), but filled up with zeros if necessary. Using the abbreviation
1

(n+1)
e (n) 2 η vec
Ω := (Ve(n) ⊗ Ue(n) )∗W
∈ Cd1d2 ,

(40)

(n)2−p
(n)2−p  1/2
(n+1)
(n+1)
e (n)− 21 W
e (n) 12 η vec
η vec
=W
⊕ Σd2
Ω
= 2−1/2 (Ve(n) ⊗ Ue(n) ) Σd1


2−p
2−p
(n)
(n)
= 2−1/2 (Ve(n) ⊗ Ue(n) ) (Id2 ⊗ Σd1 2 )D L + (Σd2 2 ⊗ Id1 )D R Ω

(41)

we rewrite

with the diagonal matrices D L , D R ∈ Md1d2 ×d1d2 such that


(D L )i+(j−1)d1,i+(j−1)d1 = 1 +
and
(D R )i+(j−1)d1,i+(j−1)d1 =

 σ 2 (X (n) ) + ϵ (n)2  2−p
 −1/2
2
j
σi2 (X (n) )

+ ϵ (n)2

  σ 2 (X (n) ) + ϵ (n)2  2−p
2
i
σ j2 (X (n) )

+ ϵ (n)2

+1

 −1/2

for i ∈ [d 1 ] and j ∈ [d 2 ]. This can be seen from the definitions of the Kronecker product ⊗ and the
Kronecker sum ⊕ (cf. Appendix A), as
 (n)2−p

(n)2−p  1/2
Σd1
⊕ Σd2
= (si + s j )1/2
i+(j−1)d 1,i+(j−1)d 1

= si (si + s j )

−1/2

+ s j (si + s j )−1/2 = si1/2 (1 +
(n)2−p

if s ` denotes the `-th diagonal entry of Σd2

(n)2−p

and Σd1

29

s j −1/2
si
)
+ s j1/2 ( + 1)−1/2
si
sj
for ` ∈ [max(d 1 , d 2 )].

(n)

2−p

If we write Σd1,Tc2 ∈ M (d1 −r )×(d1 −r ) for the diagonal matrix containing the d 1 − r last diagonal
(n)2−p

elements of Σd1

(n)

2−p

and Σd2,Tc2 ∈ M (d1 −r )×(d1 −r ) for the diagonal matrix containing the d 2 − r last
(n)2−p

diagonal elements of Σd2

, it follows from (41) that



p
−2
(n+1) (n) p
UT(n)∗
=
2
η
V
Tc Sp
c
p

2−p
(n) 2

≤ 2− 2 Σd1,Tc



UT(n)∗
Ue(n)
c

(D L Ω)mat



(n)
Σd1

2−p
2

p

(n)
(D L Ω)mat +(D R Ω)mat Σd2

+

Tc ,Tc S
p



(D R Ω)mat



Tc ,Tc

2−p
2



Ve(n)∗VT(n)
c

p
Sp

p

2−p
(n) 2

Σd2,Tc

Sp

with the notation that MTc ,Tc denotes the submatrix of M which contains the intersection of the last
d 1 − r rows of M with its last d 2 − r columns.
Now, Hölder’s inequality for Schatten-p quasinorms (e.g., [GGK00, Theorem 11.2]) can be used
to see that
2−p 
2−p p
p
p



(n)
(n)
(D L Ω)mat Tc ,Tc .
Σd1,Tc2 (D L Ω)mat Tc ,Tc
≤ ΣTc 2
(42)
S

Sp

S2

2p
2−p

Inserting the definition
(n)
ΣTc

2−p
2

p
S

=

2p
2−p

 Õ
d
i=r +1

σi2 (X (n) )

+ϵ

(n)2 

2p(2−p)
(2−p)4

 Õ
 2−p
 2−p
d
2
2
p
2
(n)
(n)2  2
σi (X )+ ϵ
=
i=r +1

allows us to rewrite the first factor, while the second factor can be bounded by


(D L Ω)mat



p
Tc ,Tc S 2

p
p
≤ kΩmat kS2 = k(Ve(n)
S2
(n+1) p
(n+1) p
e (n) 21 η vec
kW
k`2 = kη vec
k f(n) ,
` (W )

≤ (D L Ω)mat
=

1

p

(n+1)
e (n) 2 η vec
k`2
⊗ Ue(n) )∗W

2

as the matrix D L ∈ Md1d2 ×d1d2 from (41) fulfills kD L kS∞ ≤ 1 since its entries are bounded by 1; we
also recall the definition (40) of Ω and that Ve(n) and Ue(n) are unitary.
2−p p


(n)
p
The term (D R Ω)mat Tc ,Tc Σd2,Tc2
in the bound of UT(n)∗
η (n+1)VT(n)
can be estimated analSp
c
c
Sp

ogously. Combining this with (38), we obtain
kη

(n+1)

(n+1) 2p
− η 2r
kSp

≤2

−p

2−p 
d
  Õ
p
2
2
(n+1) 2p
2
(n)
(n)2  2
2
σi (X ) + ϵ
kη vec
k f(n) ,

`2 (W

i=r +1

)

concluding the proof of statement (i).
(ii) Using the strong Schatten-p null space property (18) of order 2r and that η (n+1) ∈ N (Φ), we
estimate
2p
kη (n+1) kS2 =

(n+1) 2
(n+1) 2  p
kη 2r
kS2+ kη (n+1)− η 2r
kS2 ≤
2−p

 γ 2/p + γ 2/p−1
2r

2r
2/p−1
(2r )

2−p

(n+1) 2
kη (n+1)− η 2r
kSp

p

γ
γ (γ 2r + 1)p (n+1)
(n+1) 2p
(n+1) 2p
≤ 2r
kη
− η 2r
kSp ≤ 2p 2−p2r 2−p kη (n+1) − η 2r
kSp ,
2−p
(2r )
2 r
30

where we use in the second inequality a version of Stechkin’s lemma [KKRT16, Lemma 3.1], which
leads to the estimate
2−p

kη

(n+1)

(n+1) 2
− η 2r
kS2

≤

(n+1)
kη 2r
kS2

(2r )2−p

2/p−1

p

(n+1)
kη (n+1) − η 2r
kSp ≤

γ 2r

(2r )2/p−1

(n+1) 2
kη (n+1) − η 2r
kSp .

2p

Combining the estimate for kη (n+1) kS2 with statement (i), this results in
2p
kη (n+1) kS2

≤

2−p
γ
p 2r
2 2−p
r

 Õ
d
i=r +1

σi2 (X (n) )

+ϵ

(n)2 

p
2

 2−p

(n+1)
kη vec
k

2p
,
f(n) )
`2 (W

which shows statement (ii).
(iii) For the third statement, we use the strong Schatten-p NSP (18) to see that
p

p

p

p

(n+1)
(n+1)
(n+1)
kη (n+1) kSp = kη 2r
kSp + kη (n+1) − η 2r
kSp ≤ (1 + γ 2r )kη (n+1) − η 2r
kSp ,

and combine this with statement (i).



Lemma 21. Let (X (n) )n be the output sequence of Algorithm 1 with parameters Φ, Y , r and 0 < p ≤ 1,
e (n) be the harmonic mean weight matrix matrix (10) for n ∈ N. Let X 0 ∈ Md1 ×d2 be a rank-r
and W
matrix such that Φ(X 0 ) = Y with condition number κ :=

σ1 (X 0 )
σr (X 0 ) .

(i) If (24) is fulfilled for iteration n, then η (n+1) = X (n) − X 0 fulfills
(n+1)
η vec

2p−p 2

2p
f(n) )
`2 (W

(n)
4p r p/2σr (X 0 )p(p−1) p kη kS∞
p
≤
κ
kη (n+1) kS2 .
2
2p
(n)
2p−p
(1 − ζ )
(ϵ )

(ii) Under the same assumption as for (i), it holds that
(n+1)
η vec

2p−p 2

2p
f(n) )
`2 (W

(n)
7p r p/2 max(r , d − r )p/2σr (X 0 )p(p−1) p kη kS∞
p
kη (n+1) kS∞
κ
≤
2
2p
(n)
2p−p
(1 − ζ )
(ϵ )

Proof of Lemma 21. (i) Recall that X (n+1) = arg min kX vec k 2
Φ(X )=Y
e (n) . As
matrix W

f(n) )
`2 (W

is the minimizer of the weighted

least squares problem with weight
η (n+1) = X (n+1) − X 0 is in the null space of the
measurement map Φ, it follows from Lemma 16 that
(n+1) (n+1)
(n+1)
e (n)X vec
e (n) (η (n+1) + X 0 )vec , η vec
0 = hW
, η vec i = hW
i,

which is equivalent to
(n+1)
η vec

2
f(n) )
`2 (W

(n+1) (n+1)
(n+1)
e (n)η vec
e (n) (X 0 )vec , η vec
= hW
, η vec i = −hW
i.

Using Hölder’s inequality, we can therefore estimate
(n+1)
η vec

2
f(n) )
`2 (W

(n+1)
e (n) (X 0 )vec , η vec
e (n) (X 0 )vec ]mat , η (n+1) iF
= −hW
i`2 = −h[W
 (n)

e (X 0 )vec
≤ W
kη (n+1) kS2 .
mat S 2

31

(43)

e (n) on X 0 in the matrix space as
To bound the first factor, we first rewrite the action of W
h
i

e (n)(X 0 )vec = 2[(Ve(n) ⊗ Ue(n) ) Σd(n)2−p ⊕ Σd(n)2−p −1 (Ve(n) ⊗ U
e(n) )∗ (X 0 )vec ]mat =
W
1
2
mat

e(n)∗X 0Ve(n) ) Ve(n)∗ ,
= Ue(n) H (n) ◦ (U
using (39) and Lemma 20 about the action of inverses of Kronecker sums, with the notation that
H (n) ∈ Md1 ×d2 such that
h
i
2−p
2−p −1
2
(n)
(n)2 2
2
(n)
(n)2 2
+
1
(σ
(X
)
+
ϵ
)
Hi(n)
=
2
1
(σ
(X
)
+
ϵ
)
{j
≤d
}
{i
≤d
}
j
i
j
for i ∈ [d 1 ], j ∈ [d 2 ], where 1 {i ≤d } = 1 if i ≤ d and 1 {i ≤d } = 0 otherwise. This enables us to estimate

 (n)

e (X 0 )vec
W
mat

2


e(n)∗X 0Ve(n) ) Ve(n)∗
= Ue(n) H (n)◦ (U

S2

UT(n)∗X 0VT(n) UT(n)∗X 0VT(n)
(n)
c
= H ◦ (n)∗
(n)
UTc X 0VT(n) UT(n)∗
X
V
0
Tc
c
(n)∗
(n)
= HT(n)
,T ◦ (UT X 0VT )

+ HT(n)
◦ (UT(n)∗
X 0VT(n) )
c ,T
c

2
S2
2
S2

!

2

∗
= H (n)◦ (Ue(n) X 0Ve(n) )

S2

2
S2

2
S2

(n)∗
(n)
+ HT(n)
,Tc ◦ (UT X 0VTc )

+ HT(n)
◦ (UT(n)∗
X 0VT(n)
)
c ,Tc
c
c

2

(44)

S2
2
S2

,

using the notation from the proof of Lemma 20. To bound the first summand, we calculate
(n)∗
(n)
HT(n)
,T ◦(UT X 0VT )

≤
≤
≤
≤

S2

(n)
(n)∗ (n) (n)
(n)∗ (n) (n)
≤ HT(n)
,T ◦(UT X VT ) + HT ,T ◦(−UT η VT )
S2

(n)∗ (n) (n)
◦ ΣT(n) + HT(n)
,T ◦ (UT η VT ) S
S2
2
 1/2
Õ
r
2
(n)
σi (X )
r
+ max |Hi,(n)j |kUT(n)∗η (n)VT(n) kS2

2−p
2
i, j=1
(n) ) + ϵ (n)2
i=1 σi (X
p−2
√ p−1 (n)
rσr (X ) + (σr2 (X (n) ) + ϵ (n)2 ) 2 kUT(n)∗η (n)VT(n) kS2

√ p−1 (n)
√
√ p−2
p−2
rσr (X ) + σr (X (n) ) r kη (n) kS∞ = rσr (X (n) ) σr (X (n) )

S2

HT(n)
,T


+ kη (n) kS∞ ,

denoting ΣT(n) = diag(σi (X (n) ))ri=1 and that the matrices UT(n) and VT(n) contain the first r left resp.
√
right singular vectors of X (n) in the second inequality, together with the estimates kX kS1 ≤ r kX kS2 ≤
r kX kS∞ for (r × r )-matrices X .
With the notations sr0 := σr (X 0 ) and s 10 := σ1 (X 0 ), we note that
σr (X (n) ) ≥ sr0 (1 − ζ ),
as the assumption (24) implies that
sr0 = σr (X 0 ) = σr (X (n) − η (n) ) ≤ σr (X (n) ) + σ1 (η (n) ) ≤ σr (X (n) ) + ζ sr0 ,
using [Ber09, Proposition 9.6.8] in the first inequality.
32

Therefore, we can bound the first summand of (44) such that
(n)∗
(n)
HT(n)
,T ◦ (UT X 0VT )

S2

√
√
≤ r (sr0 (1 − ζ ))p−2 [sr0 (1 − ζ ) + ζ sr0 ] = r (sr0 )p−1 (1 − ζ )p−2 .

For the second summand in the estimate of

 (n)

e (X 0 )vec
W

sumption (24) are used to compute

2
mat S
2

(45)

, similar arguments and again as-

=0
(n)∗
(n)
HT(n)
,Tc ◦ (UT X 0VTc )
(n)∗ (n) (n)
HT(n)
,Tc ◦ (UT η VTc )

≤ 

S2

S2

2kUT(n)∗η (n)VT(n)
kF
c
(X (n) )2

+ ϵ (n)2 )

≤ HT(n)
,Tc

2−p
2

≤

z
}|
{
(n)∗ (n) (n)
◦ (UT X VTc )

max

i ∈[r ]
j ∈ {r +1, ...,d 2 }

S2

+

|Hi,(n)j |kUT(n)∗η (n)VT(n)
kS2
c

(46)

(n)∗ (n) (n)
(n) p−2
 ≤ 2σr (X ) kUT η VTc kS2

(σr
√
√ 0
≤ 2 r (sr (1 − ζ ))p−2 kη (n) kS∞ ≤ 2ζ r (sr0 )p−1 (1 − ζ )p−2 .
From exactly the same arguments it follows that also
HT(n)
◦ (UT(n)∗
X 0VT(n) )
c ,T
c

S2

√
≤ 2ζ r (sr0 )p−1 (1 − ζ )p−2 .

It remains to bound the last summand HT(n)
◦ (UT(n)∗
X 0VT(n)
)
c ,Tc
c
c
◦ (UT(n)∗
X 0VT(n)
)
HT(n)
c
c ,Tc
c

S2

≤

max

i ∈ {r +1, ...,d 1 }
j ∈ {r +1, ...,d 2 }

2
S2

(47)

. We see that

X 0VT(n)
kS2
Hi,(n)j kUT(n)∗
c
c

kS∞
UT0 kS∞ kS 0 kS2 kVT0∗VT(n)
kS2 ≤ (ϵ (n) )p−2 kUT(n)∗
X 0VT(n)
≤ (ϵ (n) )p−2 kUT(n)∗
c
c
c
c
√ (n)
√ (n)
0
√ (n) 2 (n) p−2
2kη kS∞ √ 0 2kη kS∞
−2 0 −1 s 1
≤ (ϵ (n) )p−2
rs
=
2
r
kη
k
(ϵ
)
(1
−
ζ
)
(s
)
r
1
S∞
(1 − ζ )sr0
(1 − ζ )sr0
sr0

(48)

where Hölder’s inequality for Schatten norms was used in the third inequality. In the fourth inequality, Wedin’s singular value perturbation bound of Lemma 19 is used with the choice Z = X 0 ,
Z = X (n) , α = sr0 and δ = (1 − ζ )sr0 , and finally ϵ (n) ≤ ζ sr0 in the last inequality, which is implied by
the rule (14) for ϵ (n) together with assumption (24).
Summarizing the estimates (45)–(48), we conclude that

  
kη (n) kS4∞ (n) 2p−4 0 −2p s 10 2
r (sr0 )2p−2
2
1 + 8ζ + 4
(ϵ )
(sr )
≤
S2
(1 − ζ )4−2p
(1 − ζ )2p
sr0
4−2p
2p   

kη (n) kS∞ kη (n) kS∞ s 10 2
r (sr0 )2p−2
2
2p
=
(1 + 8ζ )(1 − ζ ) + 4 (n) 4−2p
(1 − ζ )4
(sr0 )2p sr0
(ϵ )
4−2p


 (n) 4−2p 
kη (n) kS∞ 2p 2
r (sr0 )2p−2
13r (sr0 )2p−2 kη kS∞ 2
≤
9
+
4
ζ
κ
≤
κ ,
(1 − ζ )4
(1 − ζ )4
(ϵ (n) )4−2p
(ϵ (n) )4−2p
 (n)

e (X 0 )vec
W
mat

2

33

as 0 < ζ < 1, ϵ (n) ≤ σr +1 (X (n) ) = kXT(n)
kS∞ ≤ kη (n) kS∞ and using the assumption (24) in the second
c
inequality. This concludes the proof of Lemma 21(i) together with inequality (43) as 13p/2 ≤ 16p/2 =
4p .
(ii) For the second statement of Lemma 21, we proceed similarly as before, but note that by
Hölder’s inequality, also
(n+1)
η vec

2
`2

f(n) )
(W

 (n)

e (X 0 )vec
W

≤

mat S 1

kη (n+1) kS∞ ,

(49)

cf. (43). Furthermore
 (n)

e (X 0 )vec
W
mat

S1

(n)∗
(n)
≤ HT(n)
,T ◦ (UT X 0VT )

+

HT(n)
c ,T

◦

(n)∗
(n)
+ HT(n)
,Tc ◦ (UT X 0VTc )

S1

(UT(n)∗
X 0VT(n) )
c
S1

HT(n)
c ,Tc

+

◦

S1

(UT(n)∗
X 0VT(n)
)
c
c
S1

.

(50)

The four Schatten-1 norms can then be estimated by max(r, (d − r ))1/2 times the corresponding
Schatten-2 norms. Using then again inequalities (45)–(48), we conclude the proof of (ii).

Proof of Theorem 11. First we note that
d
Õ
i=r +1

σi2 (X (n) ) + ϵ

p
(n)2  2

! 2−p

p2

≤ 2p− 2 (d − r )2−p σr +1 (X (n) )p(2−p)

(51)

as ϵ (n) ≤ σr +1 (X (n+1) ) due to the choice of ϵ (n) in (14). We proceed by induction over n ≥ n.
Lemma 20(ii) and Lemma 21(ii) imply together with (51) that for n = n,
2p

p
kη (n+1) kS∞

≤

kη (n+1) kS2

p
kη (n+1) kS∞

2−p

≤ 2p γ 2r 2p−
2−p

≤ 25p γ 2r

p2
2

 d − r  2−p/2 7p r p (s 0 )p(p−1)
r

(1 − ζ )2p

r

 d − r  2−p/2 r p (s 0 )p(p−1)
r

(1 − ζ )2p

r

2p−p 2

κ p kη (n) kS∞

(52)

p(2−p)

κ p kη (n) kS∞

as σr +1 (X (n) ) = ϵ (n) by assumption for n = n.
Similarly, by Lemma 20(iii), Lemma 21(ii) and (51), the error in the Schatten-p quasinorm fulfills
2p

kη (n+1) kSp ≤ (1 + γ 2r )2 22+2p d − r

 2−p r p/2 (sr0 )p(p−1) p (n) p(2−p) (n+1) p
κ kη kS∞ kη
kS2
(1 − ζ )2p

(53)

for n = n. Using the strong Schatten-p null space property of order 2r for the operator Φ, we see
from the arguments of the proof of Lemma 20(ii) that
p

1−p/2

p

kη (n) kS∞ ≤ kη (n) kS2 ≤
p

and also kη (n+1) kS2 ≤
p

1−p/2

2p−1γ 2r
r 1−p/2

2p−1γ 2r

r 1−p/2

p

kη (n) kSp

p

p

kη (n+1) kSp . Inserting that in (53) and dividing by kη (n+1) kSp , we obtain
2−p

kη (n+1) kSp ≤ 24p (1 + γ 2r )2γ 2r

 d − r  2−p r p/2 (s 0 )p(p−1)
r

r

(1
34

− ζ )2p

p(1−p)

κ p kη (n) kS∞

p

kη (n) kSp .

Under the assumption that (25) holds, it follows from this and (52) that
p

p

p

p

kη (n+1) kS∞ ≤ kη (n) kS∞ and kη (n+1) kSp ≤ kη (n) kSp

(54)

for n = n, which also entails the statement of Theorem 11 for this iteration.
0
0
Let now n 0 > n such that (54) is true for all n with n 0 > n ≥ n. If σr +1 (X (n ) ) ≤ ϵ (n −1) , then
0
0
ϵ (n ) = σr +1 (X (n ) ) and the arguments from above show (54) also for n = n 0.
0
0
0
00
00
Otherwise σr +1 (X (n ) ) > ϵ (n −1) and there exists n 0 > n 00 ≥ n such that ϵ (n ) = ϵ (n ) = σr +1 (X (n ) ).
Then
0
2−p  Õ
d  2
 p2 2−p r p/2 max(r, d − r )p/2
γ
σi (X (n ) )
p (n 0 ) p(2−p)
(n 0 +1) p
p 2r
+
1
κ
kη kS∞
kη
kS∞≤ 14 2−p
00
0
r
ϵ (n )2
(sr )p(1−p) (1 − ζ )2p
i=r +1
and we compute
 Õ
0
d  2
σ (X (n ) )
i

ϵ (n

i=r +1

 kη (n0) k p

Sp

≤

≤

00 )2

ϵ (n

00 )p

+1

 p2  2−p

+ (d − r )

 2−p

00 )p

ϵ (n

i=r +1
 kη (n00) k p
Sp

00 )p

+ (d − r )

 2−p

 2−p

+ (d − r )
00
ϵ (n )p
 2−p 

3 + γ 2r 2−p
+ (d − r )
≤
(d − r )2−p ,
1 − γ 2r
≤

 2(1 + γ )kX (n00) k p
2r
Tc
Sp
(1 − γ 2r )ϵ (n

≤

 Õ
0
p
d
σi (X (n ) )

using that X 0 is a matrix of rank at most r in the second inequality, the inductive hypothesis in
the third and an analogue of (61) for a Schatten-p quasinorm on the left hand side (cf. [KKRT16,
Lemma 3.2] for the corresponding result for p = 1) in the last inequality. The latter argument uses
the assumption on the null space property. This shows that
p

p(2−p)

kη (n +1) kS∞ ≤ µ kη (n ) kS∞
0

for
2−p

e
µ := 24p γ 2r
p

0

 (3 + γ )(d − r )  2−p r p/2 (s 0 )p(p−1)


p
2r
r
p
p
2 , (1 + γ )2 ,
κ
max
2
(d
−
r
)
2r
(1 − γ 2r )r
(1 − ζ )2p
p

and kη (n +1) kS∞ ≤ kη (n ) kS∞ under the assumption (25) of Theorem 11, as e
µ ≤ µ with µ as in (26).
Indeed e
µ ≤ µ since
0

0

  d − r  2−p
 d − r  2−p/2

p
r p/2 ≤ 2p (1 + γ 2r )2
rp.
max 2p (d − r ) 2 , (1 + γ 2r )2
r
r
p

p

The same argument shows that kη (n +1) kSp ≤ kη (n ) kSp , which finishes the proof.
0

0



Remark 22. We note that the weight matrices of the previous IRLS approaches IRLS-col and
IRLS-row [FRW11, MF12] at iteration n could be expressed in our notation as
(n)

Id2 ⊗ WL(n) := Id2 ⊗ U (n) (Σd1 )p−2U (n)∗
and

(n)

WR(n) ⊗ Id1 := V (n) (Σd2 )p−2V (n)∗ ⊗ Id1 ,
35

respectively, cf. Section 2.2, if X (n) = U (n) Σ(n)V (n)∗ = UT(n) ΣT(n)VT(n)∗ + UT(n)
ΣT(n)
V (n)∗ is the SVD of the
c
c Tc

iterate X (n) with UT(n) and VT(n) containing the r first left- and right singular vectors.
Now let

T (n) := UT(n)Z 1∗ + Z 2VT(n)∗ : Z 1 ∈ Md1 ×r , Z 2 ∈ Md2 ×r

be the tangent space of the smooth manifold of rank-r matrices at the best rank-r approximation
UT(n) ΣT(n)VT(n)∗ of X (n) , or, put differently, the direct sum of the row and column spaces of UT(n) ΣT(n)VT(n)∗ .
The fact that left- or right-sided weight matrices do not lead to algorithms with superlinear convergence rates for p < 1 can be explained by noting that there are always parts of the space T (n) that
are equipped with too large weights if X (n) = U (n) Σ(n)V (n)∗ is already approximately low-rank. In
particular, proceeding as in (44), we obtain for Id2 ⊗ WL(n)


Id2 ⊗ WL(n) (X 0 )vec



2
mat S
2

(n)  p−2

= ΣT

UT(n)∗X 0VT(n)

 p−2 (n)∗
UTc X 0VT(n)
+ ΣTc (n)

2

+

S2

+

S2

(n)

2

(n)  p−2

ΣT

UT(n)∗X 0VT(n)
c

(n)  p−2 (n)∗
ΣTc
UTc X 0VT(n)
c
(n)

2
S2

2
S2

(n)

if ΣT denotes the diagonal matrix with the first r non-zero entries of Σd1 and ΣTc the one of the
remaining entries.
Here, the third of the four summands would become too large for p < 1 to allow for a superlinear
convergence when the last d − r singular values of X (n) approach zero. An analogous argument can be
used for the right-sided weight matrix WR(n) ⊗ Id1 and, notably, also for arithmetic mean weight matrices
(n)
W(arith)
= Id2 ⊗ WL(n) + WR(n) ⊗ Id1 , cf. Section 2.3.

7

Acknowledgments

The two authors acknowledge the support and hospitality of the Hausdorff Research Institute for
Mathematics (HIM) during the early stage of this work within the HIM Trimester Program ”Mathematics of Signal Processing”. C.K. is supported by the German Research Foundation (DFG) in the
context of the Emmy Noether Junior Research Group “Randomized Sensing and Quantization of
Signals and Images” (KR 4512/1-1) and the ERC Starting Grant “High-Dimensional Sparse Optimal
Control” (HDSPCONTR - 306274). J.S. is supported by the DFG through the D-A-CH project no.
I1669-N26 and through the international research training group IGDK 1754 “Optimization and Numerical Analysis for Partial Differential Equations with Nonsmooth Structures”. The authors thank
Ke Wei for providing code of his implementations. They also thank Massimo Fornasier for helpful
discussions.

A

Kronecker and Hadamard products

For two matrices A = (ai j )i ∈[d1 ], j ∈[d3 ] ∈ Cd1 ×d3 and B ∈ Cd2 ×d4 , we call the matrix representation of
their tensor product with respect to the standard bases the Kronecker product A ⊗ B ∈ Cd1 ·d2 ×d3 ·d4 .
By its definition, A ⊗ B is a block matrix of d 2 × d 4 blocks whose block of index (i, j) ∈ [d 1 ] × [d 3 ] is
the matrix ai j B ∈ Rd2 ×d4 . This implies, e.g., for A ∈ Cd1 ×d3 with d 1 = 2 and d 3 = 3 that




a 11 B a 12 B a 13 B
a 11 a 12 a 13
A⊗B =
⊗B =
.
a 21 a 22 a 23
a 21 B a 22 B a 23 B
36

The Kronecker product is useful for the elegant formulation of matrix equations involving left
and right matrix multiplications with the variable X , as
AX B ∗ = Y

if and only if

(B ⊗ A)X vec = Yvec .

We define the Hadamard product A ◦ B ∈ Cd1 ×d2 of two matrices A ∈ Cd1 ×d2 and B ∈ Cd1 ×d2 as
their entry-wise product
(A ◦ B)i, j = Ai, j Bi, j
with i ∈ [d 1 ] and j ∈ [d 2 ]. The Hadamard product is also known as Schur product in the literature.
Furthermore, if d 1 = d 3 and d 2 = d 4 , we define the Kronecker sum A ⊕ B ∈ Cd1d2 ×d1d2 of two
matrices A ∈ Cd1 ×d1 and B ∈ Cd2 ×d2 as the matrix
A ⊕ B = (Id2 ⊗ A) + (B ⊗ Id1 ).

(55)

Note that equations of the form AX + X B ∗ = Y can be rewritten as
(A ⊕ B)X vec = Yvec ,
using again the vectorizations of X and Y . An explicit formula that expresses the inverse (A ⊕ B)−1
of the Kronecker sum A ⊕ B is provided by the following lemma.
Lemma 23 ([Jam68]). Let A ∈ Hd1 ×d1 and B ∈ Hd2 ×d2 , where one of the matrices is positive definite
and the other positive semidefinite. If we denote the singular vectors of A by ui ∈ Cd1 , i ∈ [d 1 ], its
singular values by σi , i ∈ [d 1 ] and the singular vectors resp. values of B by v j ∈ Cd2 resp. µ j , j ∈ [d 2 ],
then
d1 Õ
d2 v v ∗ ⊗ u u ∗
Õ
j j
i i
−1
.
(56)
(A ⊕ B) =
σi + µ j
i=1 j=1
Furthermore, the action of (A ⊕ B)−1 on the matrix space Md1 ×d2 can be written as



(A ⊕ B)−1Z vec mat = U H ◦ (U ∗ZV ) V ∗ .

(57)

for Z ∈ Md1 ×d2 , U = [u 1 , . . . , ud1 ], and V = [v 1 , . . . , vd2 ] and the matrix H ∈ Md1 ×d2 with the entries
Hi, j = (σi + µ j )−1 , i ∈ [d 1 ], j ∈ [d 2 ].

B
B.1

Proofs of preliminary statements in Section 6
Proof of Lemma 14: Main part

First, we define the function
p
f X,ϵ (Z )

= Jp (X , ϵ, Z ) =



 p kX vec k 2

2

f(Z ))
`2 (W

+

ϵ 2p
2

d
Í
i=1


 +∞


σi (Z ) +

2−p
2

d
Í
i=1

p

σi (Z ) (p−2)

if rank(Z ) = d,
if rank(Z ) < d,

for X ∈ Md1 ×d2 , ϵ > 0 fixed and with Z ∈ Md1 ×d2 as its only argument. We note that the set
p
of minimizers of f X,ϵ (Z ) does not contain an instance Z with rank smaller than d as the value of
p
f X,ϵ (Z ) is infinite at such points and, therefore, it is sufficient to search for minimizers on the set

Ω = Z ∈ Md1 ×d2 | rank(Z ) = d of matrices with rank d. We observe that the set Ω is an open set
and that we have that
37

k→∞

p

(a) f X,ϵ (Z ) is lower semicontinuous, which means that any sequence (Z k )k ∈N with Z k −→ Z
p
p
fulfills lim inf f X,ϵ (Z k ) ≥ f X,ϵ (Z ),
k→∞

p

(b) f X,ϵ (Z ) ≥ α for all Z ∈ Md1 ×d2 for some constant α,
k→∞

p

p

k→∞

(c) f X,ϵ (Z ) is coercive, i.e., for any sequence (Z k )k ∈N with kZ k kF −→ ∞, we have f X,ϵ (Z k ) −→
∞.
p

Property (a) is true as f X,ϵ (Z ))|Ω is a concatenation of an indicator function of an open set,
which is lower-semicontinuous and a sum of continuous functions on Ω. Property (b) is obviously
true for the choice α = 0.
d
p
ϵ 2p
ϵ 2p Í
ϵ 2p
To justify point (c), we note that f X,ϵ (Z ) > 2
σi (Z ) = 2 kZ kS1 ≥ 2 kZ kF and therefore,
i=1

coercivity nis clear from its definition.
o As a consequence from (a) and (c), it is also true that the level
p
sets LC = Z ∈ Md1 ×d2 | f X,ϵ (Z ) ≤ C are closed and bounded and therefore, compact.
Via the direct method of calculus of variations, we conclude from the properties (a) - (c) that
p
p
f X,ϵ (Z ) has at least one global minimizer belonging to the set of critical points of f X,ϵ (Z ) [Dac89,
Theorem 1].
p
To characterize the set of critical points of f X,ϵ (Z ), its derivative with respect to Z is calculated
explicitly and equated with zero in Appendix B.2. The solution of the resulting equation reveals that
p−2
Í
Í
Z opt = di=1 (σi2 (X )+ϵ 2 ) 2 ui vi∗ =: di=1 e
σi ui vi∗ is the only critical point and consequently the unique
Íd
p
L := Íd e
∗
R
global minimizer of f X,ϵ (Z ). We define the matrices Wopt
σi vi vi∗ ,
i=1 σi u i u i and Wopt :=
i=1 e

−1
R )−1 ⊕ (W L )−1
e (Z opt ) = 2 (Wopt
and note that W
with Definition 13. To verify the second part of
opt
the theorem, we simply plug the optimal solution Z opt into the functional Jp and compute using
(56) that
d
d
p
p
ϵ 2p Õ
2 − p Õ p−2
kX vec k`2 (W
+
e
σ
+
e
σ
i
2 f(Z opt ))
2
2 i=1
2 i=1 i
!
#
"
d2 Õ
d1 u u ∗ ⊗ v v ∗
d
d
d
p
Õ
j
k
pÕ 2
ϵ 2p Õ
2 − p Õ p−2
j
k
=
(u
⊗
v
)
+
e
σ
+
e
σ
σi (X )(ui∗ ⊗ vi∗ )2
i
i
i
2 i=1
2 i=1
2 i=1 i
e
σk−1 + e
σ j−1
k =1 j=1

Jp (X , ϵ, Z opt ) =

ii

d
d
p
pÕ 2
2 − p Õ p−2
=
(σi (X ) + ϵ 2 )e
σi +
e
σi
2 i=1
2 i=1

=
=

d
d
p−2
p
pÕ 2
2−p Õ 2
(σi (X ) + ϵ 2 )(σi2 (X ) + ϵ 2 ) 2 +
(σi (X ) + ϵ 2 ) 2
2 i=1
2 i=1
d
Õ
i=1

B.2

p

(σi2 (X ) + ϵ 2 ) 2 .
p

Proof of Lemma 14: Critical points of fX ,ϵ

Let us without loss of generality consider the case d = d 1 = d 2 and define
Ω = {Z ∈ Md ×d s.t. rank(Z ) = d } .
38

e (Z ) can be then rewritten as
As already mentioned in (27), the harmonic mean matrix W

e (Z ) = 2W
e1 W
e1 + W
e2 −1W
e2 = 2(W
e1−1 + W
e2−1 )−1
W
e2 = (Z ∗Z ) 12 ⊗ Id . For Z ∈ Ω, we reformulate
e1 := Id ⊗ (ZZ ∗ ) 21 and W
for Z ∈ Ω with the definitions W
the auxiliary functional such that
p

f X,ϵ (Z ) = J p (X , ϵ, Z ) =
=

d
d
p
ϵ 2p Õ
p
2−p Õ
(p−2)
kX vec k`2 (W
+
σ
(Z
)
+
σ
(Z
)
i
i
f
(Z ))
2
2
2 i=1
2 i=1

p
p
ϵ 2p
2−p
∗ 1/2 2
∗
2(p−2) k 2 .
kX vec k`2 (W
+
k(Z
Z
)
k
+
k(Z
Z
)
F
F
2 f(Z ))
2
2
2

p

To identify the set of critical points of f X,ϵ (Z ) located in Ω, we compute its derivative with respect
to Z using the derivative rules (7), (12), (13), (15), (16), (18), (20) in Chapter 8.2 and Theorem 3 in
Chapter 8.4 of [MN99] in the following. Using the notation of [MN99], we calculate
p
∂ f X,ϵ (Z )



 pϵ 2  
p  ∗ e e −1 e
∗ − 12
∗
∗ − 12 ∗
tr Z (Z Z ) ∂Z + tr((Z Z ) Z ∂Z )
= − tr X vecW ∂W W X vec +
2
4




4−p
4−p
p
∗
∗
∗
∗
2(p−2)
2(p−2)
− tr Z (Z Z )
∂Z + tr((Z Z )
Z ∂Z )
4

where

i

i
1 h
1 h ∗ −3 ∗
−1
∗ − 12
∗
∗ − 32
∗ − 21
e
2
∂W = ∂ (ZZ ) ⊕ (Z Z )
=−
(Z Z ) Z ∂Z + ∂Z Z (Z Z )
⊗ Id1
2
4

i
1h
∗ − 23 ∗
∗ − 32
∗
Id ⊗ ∂Z (ZZ ) Z + (ZZ ) Z ∂Z
−
.
4 2
We can reformulate the first term as follows using the cyclicity of the trace,

(58)


 p h 
p  ∗ e e −1 e
∗
∗ − 32 ∗
e
e
− tr X vecW ∂W W X vec =
tr (W X vec )mat (W X vec )mat (Z Z ) Z ∂Z
2
8

3
e X vec )∗mat (W
e X vec )mat ∂Z ∗
+ tr Z (Z ∗Z )− 2 (W


3
e X vec )mat (W
e X vec )∗mat ∂Z
+ tr Z ∗ (ZZ ∗ )− 2 (W

i
e X vec )mat (W
e X vec )∗mat (ZZ ∗ )− 32 Z ∂Z ∗ .
+ tr (W
p

To determine the critical points of f X,ϵ (Z ), we summarize the calculations above, rearrange the

39

terms and equate the derivative with zero, such that
p
∂ f X,ϵ (Z )

p h e
e X vec )mat (Z ∗Z )− 23 Z ∗ + Z ∗ (ZZ ∗ )− 32 (W
e X vec )∗mat
e X vec )mat (W
= tr (W X vec )∗mat (W
8
i

4−p
1
+2ϵ 2 (Z ∗Z )− 2 Z ∗ − 2(Z ∗Z ) 2(p−2) Z ∗ ∂Z
3
p h
e X vec )∗mat (W
e X vec )mat + (W
e X vec )mat (W
e X vec )∗mat (ZZ ∗ )− 23 Z
+ tr Z (Z ∗Z )− 2 (W
8

4−p i
2
∗ − 12
∗
∗
2(p−2)
+2ϵ Z (Z Z ) − 2Z (Z Z )
∂Z
:=

p
p
p
tr (A∂Z ) + tr (A∗ ∂Z ∗ ) = tr ((A ⊕ A)∂Z ) = 0,
8
8
8

where

h
e X vec )∗mat (W
e X vec )mat (Z ∗Z )− 23 Z ∗ + Z ∗ (ZZ ∗ )− 32 (W
e X vec )mat (W
e X vec )∗mat
A = (W
i
4−p
2 ∗ − 21 ∗
∗
∗
2(p−2)
+2ϵ (Z Z ) Z − 2(Z Z )
Z .

(59)

and hence an easy calculation as in [Duc] gives
p

∂ f X,ϵ (Z )
∂Z

=

p
8

tr ((A ⊕ A)∂Z ) p
= (A ⊕ A) = 0.
∂Z
8

Now we have to find Z such that A ⊕ A = 0. This implies that all eigenvalues of A ⊕ A =
A ⊗ Id + Id ⊗ A are equal to zero. The eigenvalues of the Kronecker sum of two matrices A1 and A2
with eigenvalues λs and µ t with s, t ∈ [d] are the sum of the eigenvalues λs + µ t . As in our case
A = A1 = A2 this means that all eigenvalues of A itself have to be zero. This is only possible if A is
the zero matrix.
Let Z = U ΣV ∗ ∈ Md ×d with U , V ∈ Ud and Σ ∈ Md ×d , where Σ = diag(σ ) is a diagonal matrix
2
with ascending entries. We define the matrix H = Hi, j = σ −1 +σ
−1 for i = 1, . . . , d, j = 1, . . . , d
i

j

corresponding to the result of reshaping the
 diagonal of 2(Σ ⊕ Σ) into a d × d-matrix. Using (57), we
e X vec )mat = U H ◦ (U ∗XV ) V ∗ and denote B := H ◦ (U ∗XV ).
can express (W
Plugging the decomposition Z = U ΣV ∗ into (59), we can therefore calculate

A = 0 ⇔ (U BV ∗ )∗ (U BV ∗ )(V Σ2V ∗ )−3/2 (U ΣV ∗ )∗ + (U ΣV ∗ )∗ (U Σ2U ∗ )∗ )−3/2 (U BV ∗ )(U BV ∗ )∗
4−p

+ 2ϵ 2 (V Σ2V ∗ )−1/2 (U ΣV ∗ )∗ − 2(V Σ2V ∗ ) 2(p−2) (U ΣV ∗ )∗ = 0
2

⇔ V B ∗ BΣ−2U ∗ + V Σ−2 BB ∗U ∗ + 2ϵ 2V Id U ∗ − 2V Σ p−2 U ∗ = 0
2

⇔ B ∗ BΣ−2 + Σ−2 BB ∗ + 2ϵ 2 Id − 2Σ p−2 = 0.

(60)
We now note that
− 2Σ
is diagonal and therefore,
+
is diagonal as
∗
−2
∗
2
well. Moreover, observe that B B + Σ BB Σ is again a diagonal matrix and has a symmetric first
summand B ∗ B. As the sum or difference of symmetric matrices is again symmetric also the second
summand Σ−2 BB ∗ Σ2 has to be symmetric, i.e., Σ−2 BB ∗ Σ2 = (Σ−2 BB ∗ Σ2 )∗ = Σ2 BB ∗ Σ−2 . We conclude
that it has to hold that BB ∗ Σ4 = Σ4 BB ∗ and hence Σ4 and BB ∗ commute.
2ϵ 2 Id

2
p−2

B ∗ BΣ−2

40

Σ−2 BB ∗

This is only possible if either Σ is a multiple of the identity or if BB ∗ is diagonal. Assuming the
first case, (60) would imply that also BB ∗ and B ∗ B have to be a multiple of the identity. Therefore, this
first case, where Σ is a multiple of the identity is a special case of the second possible scenario, where
BB ∗ is diagonal. Hence, it suffices to further consider the more general second case. (Considerations
for B ∗ B can be carried out analogously.)
Diagonality of BB ∗ only occurs if B is either orthonormal or diagonal. Assuming orthonormality
would lead to contradictions with the equations in (60). Hence B = H ◦(U ∗XV ) can only be diagonal.
Let now be X = Ū S¯V̄ ∗ the singular value decomposition of X . As H has no zero entries due to
the full rank of W , this implies the diagonality of U ∗Ū S¯V̄ ∗V . Consequently, U and V can only be
chosen such that P = [U ∗Ū ]d ×d and P ∗ = [V̄ ∗V ]d ×d for a permutation matrix P ∈ Ud . The reshuffled
indexing corresponding to P is denoted by p(i) ∈ [d] for i ∈ [d]. Having in mind that Hii = σi for
i ∈ [d], we obtain
2

¯ ∗ ))∗ (H ◦ (P SP
¯ ∗ ))Σ−2 + Σ−2 (H ◦ (P SP
¯ ∗ ))(H ◦ (P SP
¯ ∗ ))∗ + 2ϵ 2 Id − 2Σ p−2 = 0
(H ◦ (P SP
2

2
⇔ 2s̄p(i)
+ 2ϵ 2 = 2σip−2 for all i ∈ [d]
2
⇔ σi = (s̄p(i)
+ ϵ 2)

p−2
2

for all i ∈ [d].

As the diagonal of Σ was assumed to have ascending entries and the diagonal of S¯ has descending
entries, the permutation matrix P has to be equal to the identity matrix. From P = Id , it follows that
p−2
U = Ū and V = V̄ and hence Σ = (S¯2 + ϵ 2 Id ) 2 .
We summarize our calculations by stating that
Z opt = Ū ΣV̄ ∗ = Ū (S¯2 + ϵ 2 Id )

p−2
2

V̄ ∗

p

is the only critical point of f X,ϵ on the domain Ω.
e (Z ) is adapted by introducing
The results extend for the case d 1 , d 2 , where the definition of W
∗
1/2
the Moore-Penrose pseudo inverse of (ZZ )

e (Z ) = 2W
e1 W
e1 + W
e2 −1W
e2 = 2(W
e1+ + W
e2−1 )−1 .
W
The corresponding derivative rule as pointed out in [MN99, Chapter 8.4, Theorem 5] can be used
for the calculation in (58).

B.3

Proof of Lemma 16

The equality of the optimization problems (32) can easily be seen by the fact that only the first
e (Z ) = 2([(Z ∗Z ) 21 ]+ ⊕
summand of Jp (X , ϵ, Z ) depends on X . Now, it is important to show first that W
1
[(ZZ ∗ ) 2 ]+ )−1 is positive definite as minimizing Jp (X , ϵ, Z ) then reduces to minimizing a quadratic
Í
form. Let Z = di=1 σi ui vi∗ , where ui , vi for i ∈ [d] are the left and right singular vectors, respectively,
Í
and σi for i ∈ [d] are the singular values of Z . Since Z ∗Z = di=1 σi2vi vi∗  0, also the generalized
Í
1
1
inverse root fulfills [(ZZ ∗ ) 2 ]+  0 and for ZZ ∗ = di=1 σi2ui ui∗  0, it follows that [(ZZ ∗ ) 2 ]+  0.
1
1
We stress that at least one of the matrices (ZZ ∗ ) 2 and (Z ∗Z ) 2 is positive definite and hence also
e (Z )  0. With the fact that W
e (Z )  0, the statement can be proven analogously to the results in
W
[FRW11, Lemma 5.1].
41

B.4

Proof of Lemma 17

(a) With the minimization property that defines X (n+1) in (29), the inequality ϵ (n+1) ≤ ϵ (n) , and the
minimization property that defines Z (n+1) in (28) and Lemma 14 the monotonicity follows from
Jp (X (n) , ϵ (n) , Z (n) ) ≥ Jp (X (n+1) , ϵ (n) , Z (n) ) ≥ Jp (X (n+1) , ϵ (n+1) , Z (n) )
≥ Jp (X (n+1) , ϵ (n+1) , Z (n+1) )
(b) Using Lemma 14 and the monotonicity property of (a) for all n ∈ N, we see that
p

p

kX (n) kSp ≤ дϵ (n) (X (n) ) = Jp (X (n) , ϵ (n) , Z (n) ) ≤ Jp (X (1) , ϵ (0) , Z (0) ),
(c) The proof follows analogously to [FRW11, Proposition 6.1] where only the technical calculation
p e (n) −1
to bound σ1 ((W
) ) requires to take into account that the spectrum of a Kronecker sum A ⊕ B
consists of the pairwise sum of the spectra of A and B [Ber09, Proposition 7.2.3].

B.5

Proof of Lemma 18

e (X (n) , ϵ (n) ) = W
e (n) is clear from the definition of W
e (X , ϵ) and (10). To show
The first statement W
the necessity of (35), let X ∈ Md1 ×d2 be a critical point of (34). Without loss of generality, let us
p
assume that d 1 ≤ d 2 . In this case, a short calculation shows that дϵ (X ) = tr (XX ∗ + ϵ 2 Id1 )p/2 . It
follows from the matrix derivative rules of [MN99, Chapter 8.2, (7),(15),(18) and (20)] that
p
∇дϵ (X )

2

= p(XX + ϵ Id1 )
∗

p−2
2

d
Õ
p−2
X =p
(σi2 + ϵ 2 ) 2 σi ui vi∗ ,
i=1

Í
using the singular value decomposition X = di=1 σi ui vi∗ in the last equality. Using the Kronecker


p
e (X , ϵ)X vec
sum inversion formula (56), we see that ∇дϵ (X ) = p W
. The proof can be continued
mat
analogously to [DDFG10, Lemma 5.2].

C

Proof of Theorem 9

For statement (i) of the convergence result of Algorithm 1, we use the following reverse triangle
inequalities implied by the strong Schatten-p NSP: Let X , X 0 ∈ Md1 ×d2 such that Φ(X −X 0) = 0. Then
p

kX 0 − X kF ≤

1−p/2

2p γr
1  0 p
p
+
2β
(X
)
kX
k
−
kX
k
r
Sp ,
Sp
Sp
r 1−p/2 1 − γr

(61)

where βr (X )Sp is defined in (22). This inequality can be proven using an adaptation of the proof of
the corresponding result for `p -minimization in [GPYZ15, Theorem 13] and the generalization of
Mirksy’s singular value inequality to concave functions [Aud14, Fou18]. Furthermore, the proof of
the similar statement in [KKRT16, Theorem 12] can be adapted to show (61).
The further part of the proof of (i) as well as (ii) follow analogously to [FRW11, Theorem 6.11]
and [DDFG10, Theorem 5.3] using the preliminary results deduced in Section 6. Statement (iii) is a
direct consequence of Theorem 11, which is proven in Section 6.3.

42

References
[AR15] A. Ahmed and J. Romberg. Compressive multiplexing of correlated signals. IEEE Trans.
Inf. Theory, 61(1):479–498, 2015.
[Aud14] K. M. R. Audenaert. A generalisation of Mirsky’s singular value inequalities. preprint,
arXiv:1410.4941 [math.FA], 2014.
[Ber09] D. S. Bernstein. Matrix Mathematics: Theory, Facts, and Formulas (Second Edition).
Princeton University Press, 2009.
[BNS16] S. Bhojanapalli, B. Neyshabur, and N. Srebro. Global optimality of local search for low
rank matrix recovery. In Advances in Neural Information Processing Systems (NIPS),
pages 3873–3881, 2016.
[BTW15] J. D. Blanchard, J. Tanner, and K. Wei. CGIHT: conjugate gradient iterative hard thresholding for compressed sensing and matrix completion. Inf. Inference, 4(4):289–327, 2015.
[CDK15] J. A. Chavez-Dominguez and D. Kutzarova. Stability of low-rank matrix recovery and
its connections to Banach space geometry. J. Math. Anal. Appl., 427(1):320–335, 2015.
[CESV13] E. J. Candès, Y. Eldar, T. Strohmer, and V. Voroninski. Phase retrieval via matrix completion. SIAM J. Imag. Sci., 6(1):199–225, 2013.
[Cha07] R. Chartrand. Exact reconstructions of sparse signals via nonconvex minimization. IEEE
Signal Process. Lett., 14:707–710, 2007.
[CLS15] E. J. Candès, X. Li, and M. Soltanolkotabi. Phase retrieval via Wirtinger flow: Theory
and algorithms. IEEE Trans. Inf. Theory, 61(4):1985–2007, 2015.
[CP11] E. J. Candès and Y. Plan. Tight oracle inequalities for low-rank matrix recovery from a
minimal number of noisy random measurements. IEEE Trans. Inf. Theory, 57(4):2342–
2359, April 2011.
[CR09] E. J. Candès and B. Recht. Exact matrix completion via convex optimization. Found.
Comput. Math., 9(6):717–772, 2009.
[CSV13] E. J. Candès, T. Strohmer, and V. Voroninski. PhaseLift: Exact and Stable Signal Recovery
from Magnitude Measurements via Convex Programming. Commun. Pure Appl. Math.,
66(8):1241–1274, 2013.
[Dac89] B. Dacorogna. Direct Methods in the Calculus of Variations. Springer, New York, 1989.
[DDFG10] I. Daubechies, R. DeVore, M. Fornasier, and C.S. Güntürk. Iteratively reweighted least
squares minimization for sparse recovery. Commun. Pure Appl. Math., 63:1–38, 2010.
[DGM13] D. L. Donoho, M. Gavish, and A. Montanari. The phase transition of matrix recovery
from Gaussian measurements matches the minimax MSE of matrix denoising. Proc. Nat.
Acad. Sci. U.S.A., 110(21):8405–8410, 2013.
[DR16] M. A. Davenport and J. Romberg. An overview of low-rank matrix recovery from incomplete observations. IEEE J. Sel. Topics Signal Process., 10:608–622, 06 2016.
43

[Duc] J. Duchi. Properties of the Trace and Matrix Derivatives. Available electronically
at https://web.stanford.edu/∼jduchi/projects/matrix prop.
pdf.
[ENP12] Y.C. Eldar, D. Needell, and Y. Plan. Uniqueness conditions for low-rank matrix recovery.
Appl. Comput. Harmon. Anal., 33(2):309–314, 2012.
[Faz02] M. Fazel. Matrix rank minimization with applications. Ph.D. Thesis, Electrical Engineering Department, Stanford University, 2002.
[Fou18] S. Foucart. Concave Mirsky Inequality and Low-Rank Recovery. SIAM J. Matrix Anal.
Appl., 39(1):99–103, 2018.
[FPRW16] M. Fornasier, S. Peter, H. Rauhut, and S. Worm. Conjugate gradient acceleration of
iteratively re-weighted least squares methods. Comput. Optim. Appl., 65(1):205–259,
2016.
[FR13] S. Foucart and H. Rauhut. A Mathematical Introduction to Compressive Sensing. Applied
and Numerical Harmonic Analysis. Birkhäuser/Springer, New York, 2013.
[FRW11] M. Fornasier, H. Rauhut, and R. Ward. Low-rank matrix recovery via iteratively
reweighted least squares minimization. SIAM J. Optim., 21(4):1614–1640, 2011. code
from https://github.com/rward314/IRLSM].
[GB14] M. Grant and S. Boyd. CVX: Matlab software for disciplined convex programming,
version 2.1. http://cvxr.com/cvx, March 2014.
[GGK00] I. Gohberg, S. Goldberg, and N. Krupnik. Traces and determinants of linear operators,
volume 116 of Operator Theory: Advances and Applications. Birkhäuser, Basel, 2000.
[GKK15] D. Gross, F. Krahmer, and R. Kueng. A partial derandomization of phaselift using spherical designs. J. Fourier Anal. Appl., 21(2):229–266, 2015.
[GLF+ 10] D. Gross, Y.-K. Liu, S. T. Flammia, S. Becker, and J. Eisert. Quantum state tomography
via compressed sensing. Phys. Rev. Lett., 105:150401, 2010.
[GLM16] R. Ge, J. D. Lee, and T. Ma. Matrix completion has no spurious local minimum. In
Advances in Neural Information Processing Systems (NIPS), pages 2973–2981, 2016.
[GNOT92] D. Goldberg, D. Nichols, B. M. Oki, and D. Terry. Using collaborative filtering to weave
an information tapestry. Commun. ACM, 35(12):61–70, 1992.
[GPYZ15] Y. Gao, J. Peng, S. Yue, and Y. Zhao. On the null space property of `q -minimization for
0 < q ≤ 1 in compressed sensing. J. Funct. Spaces, 2015:4203–4215, 2015.
[Gro11] D. Gross. Recovering low-rank matrices from few coefficients in any basis. IEEE Trans.
Inf. Theory, 57(3):1548–1566, 2011.
[HH09] J. P. Haldar and D. Hernando. Rank-constrained solutions to linear matrix equations
using powerfactorization. IEEE Signal Process. Lett., 16(7):584–587, July 2009. [using
AltMin (Alternating Minimization) algorithm].
44

[HMT11] N. Halko, P. G. Martinsson, and J. A. Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Rev.,
53(2):217–288, 2011.
[Jam68] A. Jameson. Solution of the equation ax + xb = c by inversion of an m x m or n x n
matrix. SIAM J. Appl. Math., 16(5):1020–1023, 1968.
[JMD10] P. Jain, Raghu M., and I. S. Dhillon. Guaranteed rank minimization via singular value
projection. In Advances in Neural Information Processing Systems (NIPS), pages 937–945,
2010.
[JNS13] P. Jain, P. Netrapalli, and S. Sanghavi. Low-rank matrix completion using alternating
minimization. In Proc. ACM Symp. Theory Comput. (STOC), pages 665–674, Palo Alto,
CA, USA, June 2013.
[KC14] A. Kyrillidis and V. Cevher. Matrix recipes for hard thresholding methods. J. Math.
Imaging Vision, 48(2):235–265, 2014. [using Matrix ALPS II (’Matrix ALgrebraic PursuitS II’) algorithm, code from http://akyrillidis.github.io/
projects/].
[KKRT16] M. Kabanava, R. Kueng, H. Rauhut, and U. Terstiege. Stable low-rank matrix recovery
via null space properties. Inf. Inference, 5(4):405–441, 2016.
[KS17] C. Kümmerle and J. Sigl. Harmonic Mean Iteratively Reweighted Least Squares for
low-rank matrix recovery. In 12th International Conference on Sampling Theory and
Applications (SampTA), pages 489–493, 2017.
[KTT15] F. J. Király, L. Theran, and R. Tomioka. The Algebraic Combinatorial Approach for LowRank Matrix Completion. J. Mach. Learn. Res., 16:1391–1436, 2015.
[LHV13] Z. Liu, A. Hansson, and L. Vandenberghe. Nuclear norm system identification with
missing inputs and outputs. Systems Control Lett., 62(8):605–612, 2013.
[LV10] Z. Liu and L. Vandenberghe. Interior-point method for nuclear norm approximation
with application to system identification. SIAM J. Matrix Anal. Appl., 31(3):1235–1256,
2010.
[MF12] K. Mohan and M. Fazel. Iterative reweighted algorithms for matrix rank minimization.
J. Mach. Learn. Res., 13(1):3441–3473, 2012. [using IRLS-MF (’IRLS-p’) algorithm, code
from https://faculty.washington.edu/mfazel/].
[MMBS13] B. Mishra, G. Meyer, F. Bach, and R. Sepulchre. Low-rank optimization with trace norm
penalty. SIAM J. Optim., 23(4):2124–2149, 2013.
[MN99] J.R. Magnus and H. Neudecker. Matrix Differential Calculus with Applications in Statistics
and Econometrics. Wiley Series in Probability and Statistics. Wiley, 1999.
[OMFH11] S. Oymak, K. Mohan, M. Fazel, and B. Hassibi. A simplified approach to recovery conditions for low rank matrices. In Proceedings of the IEEE International Symposium on
Information Theory (ISIT), pages 2318–2322, 2011.
45

[PABN16] D.L. Pimentel-Alarcón, N. Boston, and R. D. Nowak. A Characterization of Deterministic Sampling Patterns for Low-Rank Matrix Completion. preprint, arXiv:1503.02596v3
[stat.ML], October 2016.
[PKCS16] D. Park, A. Kyrillidis, C. Caramanis, and S. Sanghavi.
Finding Low-rank Solutions to Matrix Problems, Efficiently and Provably. preprint, arXiv:1606.03168
[math.OC], [using BFGD (’Bi-Factored Gradient Descent’) algorithm, code from http:
//akyrillidis.github.io/projects/], 2016.
[RFP10] B. Recht, M. Fazel, and P. A. Parrilo. Guaranteed minimum-rank solutions of linear
matrix equations via nuclear norm minimization. SIAM Rev., 52(3):471–501, 2010.
[RXH11] B. Recht, W. Xu, and B. Hassibi. Null space conditions and thresholds for rank minimization. Math. Program., 127(1):175–202, 2011.
[SL16] R. Sun and Z. Q. Luo. Guaranteed matrix completion via non-convex factorization. IEEE
Trans. Inf. Theory, 62(11):6535–6579, 2016.
[SRJ05] N. Srebro, J. Rennie, and T. S. Jaakkola. Maximum-margin matrix factorization. In
Advances in Neural Information Processing Systems (NIPS), pages 1329–1336, 2005.
[SS16] É. Schost and P.-J. Spaenlehauer. A quadratically convergent algorithm for structured
low-rank approximation. Found. Comput. Math., 16(2):457–492, 2016.
[Ste06] M. Stewart. Perturbation of the SVD in the presence of small singular values. Linear
Algebra Appl., 419(1):53–77, 2006.
[TBS+ 15] S. Tu, R. Boczar, M. Simchowitz, M. Soltanolkotabi, and B. Recht. Low-rank Solutions
of Linear Matrix Equations via Procrustes Flow. preprint, arXiv:1507.03566 [math.OC],
2015.
[TW13] J. Tanner and K. Wei. Normalized Iterative Hard Thresholding for Matrix Completion.
SIAM J. Sci. Comput., 35(5):S104–S125, 2013.
[TW16] J. Tanner and K. Wei. Low rank matrix completion by alternating steepest descent
methods. Appl. Comput. Harmon. Anal., 40(2):417–429, 2016. [using ASD (’Alternating
Steepest Descent’) algorithm, code from https://www.math.ucdavis.edu/
∼kewei/publications.html].
[Van13] B. Vandereycken.
Low-rank matrix completion by Riemannian optimization.
SIAM J. Optim., 23(2):1214–1236, 2013.
[using Riemann Opt (’Riemannian Optimization’) algorithm, code from http://www.unige.ch/math/
vandereycken/matrix completion.html].
[WCCL16] K. Wei, J.-F. Cai, T. F. Chan, and S. Leung. Guarantees of Riemannian Optimization for
Low Rank Matrix Recovery. SIAM J. Matrix Anal. Appl., 37(3):1198–1222, 2016.
[Wed72] P.-Å. Wedin. Perturbation bounds in connection with singular value decomposition.
BIT, 12(1):99–111, 1972.

46

[WYZ12] Z. Wen, W. Yin, and Y. Zhang. Solving a low-rank factorization model for matrix completion by a nonlinear successive over-relaxation algorithm. Math. Program. Comput.,
4(4):333–361, 2012.
[ZL15] Q. Zheng and J. Lafferty. A convergent gradient descent algorithm for rank minimization and semidefinite programming from random linear measurements. In Advances in
Neural Information Processing Systems (NIPS), pages 109–117, 2015.

47

