remarkRemark hypothesisHypothesis claimClaim

# Recovering Simultaneously Structured Data via Non-Convex Iteratively
Reweighted Least Squares

Christian Kümmerle Department of Computer Science, University of North
Carolina at Charlotte, Charlotte, NC 28223, USA (, http://ckuemmerle.com).
<EMAIL_ADDRESS>Johannes Maly Department of Mathematics, Ludwig-
Maximilians-Universität München, 80799 Munich, Germany (, https://johannes-
maly.github.io<EMAIL_ADDRESS>
# Recovering Simultaneously Structured Data via Non-Convex Iteratively
Reweighted Least Squares

Christian Kümmerle Department of Computer Science, University of North
Carolina at Charlotte, Charlotte, NC 28223, USA (, http://ckuemmerle.com).
<EMAIL_ADDRESS>Johannes Maly Department of Mathematics, Ludwig-
Maximilians-Universität München, 80799 Munich, Germany (, https://johannes-
maly.github.io<EMAIL_ADDRESS>
###### Abstract

We propose a new algorithm for the problem of recovering data that adheres to
multiple, heterogeneous low-dimensional structures from linear observations.
Focusing on data matrices that are simultaneously row-sparse and low-rank, we
propose and analyze an iteratively reweighted least squares (IRLS) algorithm
that is able to leverage both structures. In particular, it optimizes a
combination of non-convex surrogates for row-sparsity and rank, a balancing of
which is built into the algorithm. We prove locally quadratic convergence of
the iterates to a simultaneously structured data matrix in a regime of minimal
sample complexity (up to constants and a logarithmic factor), which is known
to be impossible for a combination of convex surrogates. In experiments, we
show that the IRLS method exhibits favorable empirical convergence,
identifying simultaneously row-sparse and low-rank matrices from fewer
measurements than state-of-the-art methods.

##  1 Introduction

Reconstructing an image from (noisy) linear observations is maybe the most
relevant inverse problem for modern image processing and appears in various
applications like medical imaging and astronomy [BBDM21]. If the latent image
is $n$-dimensional, for $n\in\mathbb{N}$, it is well-known that $\Omega(n)$
observations are required for robust identification in general. In practice,
imaging problems are however often ill-posed, i.e., the number of observations
is smaller than $n$ or the operator creating the observations is defective
[Tik63, Vog02]. In such situations, the fundamental lower bound of $\Omega(n)$
can be relaxed by leveraging structural priors of the latent image in the
reconstruction process.

Of the various priors that are used for solving ill-posed inverse problems in
the literature, sparsity111A vector $\mathbf{x}\in\mathbb{R}^{n}$ is called
$s$-sparse if $\mathbf{x}$ has at most $s$ non-zero entries. For a matrix
$\mathbf{X}\in\mathbb{R}^{n_{1}\times n_{2}}$ there are various ways to count
the level of sparsity. In this work, we use the most common definition and
call $\mathbf{X}$ $s$-row-sparse (resp. -column-sparse) if at most $s$ rows
(resp. columns) are non-zero. and low-rankness are most prevalent. This
prominent role can be explained with their competitive performance in imaging
tasks and the rigorous mathematical analysis they allow [BDE09, MBP+14]. For
instance, consider the recovery of an $n_{1}\times n_{2}$-dimensional image
$\mathbf{X}_{\star}\in\mathbb{R}^{n_{1}\times n_{2}}$ from linear observations

(1)
$\displaystyle\mathbf{y}=\mathcal{A}(\mathbf{X}_{\star})+\boldsymbol{\eta}\in\mathbb{R}^{m},$

where $\mathcal{A}\colon\mathbb{R}^{n_{1}\times n_{2}}\to\mathbb{R}^{m}$ is a
linear operator modeling the impulse response of the sensing device and
$\boldsymbol{\eta}\in\mathbb{R}^{m}$ models additive noise. Whereas this
problem is ill-posed for $m<n_{1}n_{2}$, it has been established [FR13, CP11,
RFP10] that it becomes well-posed if $\mathbf{X}_{\star}$ is sparse or of low
rank. The aforementioned works prove that $m=\Omega(s_{1}s_{2})$ observations
suffice for robust reconstruction if $\mathbf{X}_{\star}$ is $s_{1}$-row-
sparse and $s_{2}$-column-sparse, and that $m=\Omega(r(n_{1}+n_{2}))$
observations suffice if $\mathbf{X}_{\star}$ is a rank-$r$ matrix. These
bounds, which relax the general lower bound of $m=\Omega(n_{1}n_{2})$, agree
with the degrees of freedom of sparse and low-rank matrices, respectively.

A number of computationally challenging problems in signal processing and
machine learning can be formulated as instances of Eq. 1 with
$\mathbf{X}_{\star}$ being _simultaneously structured_ , i.e.,
$\mathbf{X}_{\star}$ is both of rank $r$ and $s_{1}$-row-
sparse/$s_{2}$-column-sparse. Examples encompass sparse phase retrieval
[KST95, JOH17, CLM16, JH17, Sol19], sparse blind deconvolution [LLJB16, SC21],
hyperspectral imaging [HV19, GV12, TRB17, ZHL+22], sparse reduced-rank
regression [CH12, YGK20], and graph denoising and refinement [RSV12, ZZL22].
In these settings, the hope is that due to leveraging the simultaneous
structure, $\Omega(r(s_{1}+s_{2}))$ observations suffice to identify the data
matrix. For $r\ll s_{1},s_{2}\ll n_{1},n_{2}$, these bounds are significantly
smaller than the bounds for single-structured data matrices.

From an algorithmic point of view, however, the simultaneously structured
recovery problem poses obstacles that are _not_ present for problems where
$\mathbf{X}_{\star}$ is only of low-rank _or_ (group) sparse: In the latter
case, variational methods [BB18] that formulate the reconstruction method in
terms of optimizing a suitable objective with a structural regularization term
involving $\ell_{p}$/$\ell_{2,p}$-(quasi-)norms and and $S_{p}$-Schatten
(quasi-)norms have been well-understood, leading to tractable algorithms in
the information theoretically optimal regime [Don06, Cha07, FR13].
For simultaneously structured problems, on the other hand, Oymak et al. showed
in [OJF+15] that a mere linear combination of _convex_ regularizers for
different sparsity structures — in our case, nuclear and $\ell_{2,1}$-norms —
cannot outperform recovery guarantees of the “best” one of them alone. While
this indicates that leveraging two priors at once is a way more intricate
problem than leveraging a single prior, it was also shown in [OJF+15] that
minimizing a linear combination of rank and row sparsity _can_ indeed lead to
guaranteed recovery from $\Omega(r(s_{1}+s_{2}))$ measurements. The downside
is that the combination of these _non-convex and discontinuous_ quantities
does not lend itself directly to practical optimization algorithms, and to the
best of our knowledge, so far, there have been no works directly tackling the
optimization of a combination of _non-convex surrogates_ that come with any
sort of convergence guarantees.

###  1.1 Contribution

In this work, we approach the reconstruction of simultaneously sparse and low-
rank matrices by leveraging the positive results of [OJF+15] for non-convex
regularizers. To this end, we introduce a family of non-convex, but
continuously differentiable regularizers that are tailored to the recovery
problem for simultaneously structured data. The resulting objectives lend
themselves to efficient optimization by a novel algorithm from the class of
_iteratively reweighted least squares ( IRLS)_ [DDFG10, FRW11, MF12, BS15,
APS19, KMV21, KMVS21], the convergence of which we analyze in the information
theoretically (near-)optimal regime. Specifically, our main contributions are
threefold:

  1. (i)

In Algorithm 1, we propose a novel IRLS method that is tailored to leveraging
both structures of the latent solution, sparsity and low-rankness, at once.
The core components of the algorithm are the weight operator defined in
Definition 2.1 and the update of the smoothing parameters in Eq. 12. Notably,
the algorithm automatically balances between its low-rank and its sparsity
promoting terms, leading to a reliable identification of $s_{1}$-row-sparse
and rank-$r$ ground truths222To enhance readability of the presented proofs
and results, we restrict ourselves to row-sparsity of $\mathbf{X}_{\star}$
here. It is straight-forward to generalize the arguments to column-sparsity as
well..

  2. (ii)

Under the assumption that $\mathcal{A}$ behaves almost isometrically on the
set of row-sparse and low-rank matrices, we show in Theorem 2.4 that locally
Algorithm 1 exhibits quadratic convergence towards $\mathbf{X}_{\star}$. Note
that if $\mathcal{A}$ is, e.g., a Gaussian operator, the isometry assumption
is (up to log factors) fulfilled in the information theoretic (near-)optimal
regime $m=\Omega(r(s_{1}+n_{2}))$ [LWB18].

  3. (iii)

Finally, in Section 2.3 we identify the underlying family of objectives that
are minimized by Algorithm 1. To make this precise, we define for $\tau>0$ and
$e$ denoting Euler’s number the real-valued function
$f_{\tau}:\mathbb{R}\to\mathbb{R}$ such that

(2)
$\begin{split}f_{\tau}(t)=\begin{cases}\frac{1}{2}\tau^{2}\log(et^{2}/\tau^{2}),&\text{
if }|t|>\tau,\\\ \frac{1}{2}t^{2},&\text{ if
}|t|\leq\tau,\end{cases},\end{split}$

which is quadratic around the origin and otherwise logarithmic in its
argument. Using this definition, we define for $\varepsilon>0$ the
_$(\varepsilon-)$ smoothed $\log$-determinant_ objective
$\mathcal{F}_{\textit{lr},\varepsilon}:\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}$ and for $\delta>0$ the _$(\delta-$)smoothed sum of
logarithmic row-wise $\ell_{2}$-norms_ objective
$\mathcal{F}_{\textit{sp},\delta}:\mathbb{R}^{n_{1}\times n_{2}}\to\mathbb{R}$
such that

(3)
$\displaystyle\mathcal{F}_{\textit{lr},\varepsilon}(\mathbf{X})=\sum_{r=1}^{\min\\{n_{1},n_{2}\\}}f_{\varepsilon}(\sigma_{r}(\mathbf{X})),\qquad\qquad\mathcal{F}_{\textit{sp},\delta}(\mathbf{X})=\sum_{i=1}^{n_{1}}f_{\delta}(\|\mathbf{X}_{i,:}\|_{2}).$

Combining the above, we further define the _$(\varepsilon,\delta-)$ smoothed
logarithmic surrogate_ objective
$\mathcal{F}_{\varepsilon,\delta}:\mathbb{R}^{n_{1}\times n_{2}}\to\mathbb{R}$
as

(4)
$\displaystyle\mathcal{F}_{\varepsilon,\delta}(\mathbf{X}):=\mathcal{F}_{\textit{lr},\varepsilon}(\mathbf{X})+\mathcal{F}_{\textit{sp},\delta}(\mathbf{X}).$

In Theorem 2.5, we prove for _any_ $\mathcal{A}$ that the iterates of
Algorithm 1 minimize quadratic majorizations of
$\mathcal{F}_{\varepsilon,\delta}$ and form a non-increasing sequence on
$\mathcal{F}_{\varepsilon,\delta}$. To the best of our knowledge, the proposed
method is the so far only approach for recovering simultaneously sparse and
low-rank matrices which combines local (quadratic) convergence with a rigorous
variational interpretation.

The numerical simulations in Section 3 support our theoretical findings and
provide empirical evidence for the efficacy of the proposed method.

###  1.2 Related Work

In this section, we review literature relevant for our work.

#### Sparse and Low-Rank Recovery

Whereas leveraging a single matrix structure like sparsity _or_ low-rankness
in the reconstruction process can easily be obtained by convex regularizers
[RFP10, CRPW12], Oymak et al. [OJF+15] showed that, if one is interested in
near-optimal sampling rates, one cannot expect comparably simple solutions for
identifying simultaneously structured objects; a minimization of Eq. 4 with
_convex_ terms $\mathcal{F}_{\textit{lr},\varepsilon}(\cdot)$ and
$\mathcal{F}_{\textit{sp},\delta}(\cdot)$ would be only as good as using the
one structure that is information theoretically more favorable. A closely
related problem333Since observations in SPCA are provided from noisy samples
of the underlying distribution, whereas in our case the matrix itself is
observed indirectly, it is however hard to directly compare results from
sparse and low-rank matrix reconstruction with corresponding results for SPCA.
that appears in statistical literature under the name _Sparse Principal
Component Analysis_ (SPCA) [ZHT06, dGJL05] is known to be NP-hard in general
[MI17]. Due to this intrinsic hardness of reconstructing simultaneously sparse
and low-rank matrices there exist only a handful approaches that come with
rigorous theoretical analysis.
The first line of works [BR16, FGJR20] aims to overcome the aforementioned
limitations of purely convex methods in a neat way. They assume that the
operator $\mathcal{A}$ has a nested structure such that basic solvers for low-
rank resp. row/column-sparse recovery can be applied in two consecutive steps.
Despite being an elegant idea, this approach clearly restricts possible
choices for $\mathcal{A}$ and is of hardly any practical use.
In a second line of work, Lee et al. [LWB18] consider general impulse response
operators that satisfy a suitable restricted isometry property for
$s_{1}$-row- and $s_{2}$-column-sparse rank-$r$ matrices. They propose and
analyze a highly efficient, greedy method, the so-called Sparse Power
Factorization (SPF) which is a modified version of power factorization [JNS13]
and uses hard thresholding pursuit [Fou11] to enforce sparsity in addition. In
particular, they show that if $\mathbf{X}_{\star}$ is rank-$R$, has
$s_{1}$-sparse columns and $s_{2}$-sparse rows, then $m\gtrsim
R(s_{1}+s_{2})\log(\max\\{en_{1}/s_{1},en_{2}/s_{2}\\})$ Gaussian observations
suffice for robust recovery, which is up to the log-factor at the information
theoretical limit we discussed above. The result however assumes a low noise
level and requires that SPF is initialized by a, in general, intractable
method. Only in the special case that $\mathbf{X}_{\star}$ is spiky, which
means that the norms of non-zero rows/columns exhibit a fast decay, a
tractable substitute for the initialization method works provably. The
analysis of SPF has been extended to the blind deconvolution setup in
[LLJB16]. In [YGK20], a related approach that combines gradient descent of a
smooth objective with hard thresholding is considered, for which the authors
show linear convergence from a suitable intialization if the measurement
operator satisfies a restricted strong convexity and smoothness assumption.
A third line of work, approaches the problem from a variational point of view.
In [FMN21, Mal21] the authors aim at enhancing robustness of recovery by
alternating minimization of an $\ell_{1}$-norm based multi-penalty functional.
In essence, the theoretical results bound the reconstruction error of global
minimizers of the proposed functional depending on the number of observations.
Although the authors only provide local convergence guarantees for the
proposed alternating methods, the theoretical error bounds for global
minimizers hold for arbitrarily large noise magnitudes and a wider class of
ground-truth matrices than the one considered in [LWB18].
The works [FGJR20, EKPU22], which build upon generalized projection operators
to modify iterative hard thresholding to the simultaneous setting, share the
lack of global convergence guarantees.
Finally, the alternative approach of using optimally weighted sums or maxima
of convex regularizers [KSJ19] requires optimal tuning of the parameters under
knowledge of the ground-truth.

#### Iteratively Reweighted Least Squares

The herein proposed iteratively reweighted least squares algorithm builds on a
long line of research on IRLS going back to Weiszfeld’s algorithm proposed in
the 1930s for a facility location problem [Wei37, BS15]. IRLS is a practical
framework for the optimization of non-smooth, possibly non-convex, high-
dimensional objectives that minimizes quadratic models which majorize these
objectives. Due to its ease of implementation and favorable data-efficiency,
it has been widely used in compressed sensing [GR97, CY08, DDFG10, LXY13,
FPRW16, KMVS21], robust statistics [HW77, AH15, MGJK19], computer vision
[CG17, LC22, SWL22], low-rank matrix recovery and completion [FRW11, MF12,
KS18, KMV21], and in inverse problems involving group sparsity [CHHL14, ZB15,
CHLH18]. Recently, it has been shown [LK23] that dictionary learning
techniques can be incorporated into IRLS schemes for sparse and low-rank
recovery to allow the learning of a sparsifying dictionary while recovering
the solution. Whereas IRLS can be considered as a type of majorize-minimize
algorithm [Lan16], optimal performance is achieved if intertwined with a
smoothing strategy for the original objective, in which case globally linear
(for convex objectives) [DDFG10, APS19, MGJK19, KMVS21, PKV22] and locally
superlinear (for non-convex objectives) [DDFG10, KS18, KMV21, PKV22]
convergence rates have been shown under suitable conditions on the linear
operator $\mathcal{A}$.

However, there has only been little work on IRLS optimizing a sum of
heterogenous objectives [SM17] — including the combination of low-rank
promoting and sparsity-promoting objectives — nor on the convergence analysis
of any such methods. The sole algorithmically related approach for our setting
has been studied in [Che18], where a method has been derived in a sparse
Bayesian learning framework, the main step of which amounts to the
minimization of weighted least squares problems. Whereas the algorithm of
[Che18] showcases that such a method can empiricially identify simultaneously
structured matrices from a small number of measurements, no convergence
guarantees or rates have been provided in the information-theoretically
optimal regime. Furthermore, [Che18] only focuses on general sparsity rather
than row or column sparsity.

###  1.3 Notation

We denote matrices and vectors by bold upper- and lower-case letters to
distinguish them from scalars and functions. We furthermore denote the $i$-th
row of a matrix $\mathbf{Z}\in\mathbb{R}^{n_{1}\times n_{2}}$ by
$\mathbf{Z}_{i,:}$ and the $j$-th column of $\mathbf{Z}$ by
$\mathbf{Z}_{:,j}$. We abbreviate $n=\min\\{n_{1},n_{2}\\}$. We denote the
$r$-th singular value of a matrix $\mathbf{Z}\in\mathbb{R}^{n_{1}\times
n_{2}}$ by $\sigma_{r}(\mathbf{Z})$. Likewise, we denote the in
$\ell_{2}$-norm $s$-largest row of $\mathbf{Z}$ by $\rho_{s}(\mathbf{Z})$. (To
determine $\rho_{s}(\mathbf{Z})$, we form the in $\ell_{2}$-norm non-
increasing rearrangement of the rows of $\mathbf{Z}$ and, by convention, sort
rows with equal norm according to the row-index.) We use $\circ$ to denote the
Hadamard (or Schur) product, i.e., the entry-wise product of two
vectors/matrices.
We denote the Euclidean $\ell_{2}$-norm of a vector
$\mathbf{z}\in\mathbb{R}^{n}$ by $\|\mathbf{z}\|_{2}$. For
$\mathbf{Z}\in\mathbb{R}^{n_{1}\times n_{2}}$, the matrix norms we use
encompass the operator norm $\|\mathbf{Z}\|\leavevmode\nobreak\
:=\leavevmode\nobreak\
\sup_{\|\mathbf{w}\|_{2}=1}\|\mathbf{Z}\mathbf{w}\|_{2}$, the row-sum
$p$-quasinorm $\|\mathbf{Z}\|_{p,2}\leavevmode\nobreak\ :=\leavevmode\nobreak\
\big{(}\sum_{i=1}^{n_{1}}\|\mathbf{Z}_{i,:}\|_{2}^{p}\big{)}^{1/p}$, the row-
max norm $\|\mathbf{Z}\|_{\infty,2}\leavevmode\nobreak\ :=\leavevmode\nobreak\
\max_{i\in[n_{1}]}\|\mathbf{Z}_{i,:}\|_{2}$, and the Schatten-p quasinorm
$\|\mathbf{Z}\|_{S_{p}}\leavevmode\nobreak\ :=\leavevmode\nobreak\
\left(\sum_{r=1}^{n}\sigma_{r}(\mathbf{Z})^{p}\right)^{\nicefrac{{1}}{{p}}}$.
Note that two special cases of Schatten quasinorms are the nuclear norm
$\|\mathbf{Z}\|_{*}\leavevmode\nobreak\ :=\leavevmode\nobreak\
\|\mathbf{Z}\|_{S_{1}}$ and the Frobenius norm
$\|\mathbf{Z}\|_{F}\leavevmode\nobreak\ :=\leavevmode\nobreak\
\|\mathbf{Z}\|_{S_{2}}$.

##  2 IRLS for Sparse and Low-Rank Reconstruction

Recall that we are interested in recovering a rank-$r$ and $s$-row-sparse
matrix $\mathbf{X}_{\star}\in\mathbb{R}^{n_{1}\times n_{2}}$ from $m$ linear
observations

(5)
$\displaystyle\mathbf{y}=\mathcal{A}(\mathbf{X}_{\star})\in\mathbb{R}^{m},$

i.e., $\mathcal{A}\colon\mathbb{R}^{n_{1}\times n_{2}}\to\mathbb{R}^{m}$ is
linear. We write
$\mathbf{X}_{\star}\in\mathcal{M}_{r,s}^{n_{1},n_{2}}:=\mathcal{M}_{r}^{n_{1},n_{2}}\cap\mathcal{N}_{s}^{n_{1},n_{2}}$,
where $\mathcal{M}_{r}^{n_{1},n_{2}}\subset\mathbb{R}^{n_{1}\times n_{2}}$
denotes the set of matrices with rank at most $r$ and
$\mathcal{N}_{s}^{n_{1},n_{2}}\subset\mathbb{R}^{n_{1}\times n_{2}}$ denotes
the set of matrices with at most $s$ non-zero rows. For convenience, we
suppress the indices $n_{1}$ and $n_{2}$ whenever the ambient dimension is
clear from the context. In particular, we know that
$\mathbf{X}_{\star}=\mathbf{U}_{\star}\boldsymbol{\Sigma}_{\star}\mathbf{V}_{\star}^{*}$,
where $\mathbf{U}_{\star}\in\mathcal{N}_{s}^{n_{1},r}$,
$\boldsymbol{\Sigma}_{\star}\in\mathbb{R}^{r\times r}$, and
$\mathbf{V}_{\star}\in\mathbb{R}^{n_{2}\times r}$ denote the reduced SVD of
$\mathbf{X}_{\star}$. Furthermore, the row supports of $\mathbf{X}_{\star}$
and $\mathbf{U}_{\star}$ (the index sets of non-zero rows of
$\mathbf{X}_{\star}$ resp. $\mathbf{U}_{\star}$) are identical, i.e.,
$\mathrm{supp}(\mathbf{X}_{\star})=\mathrm{supp}(\mathbf{U}_{\star})=S_{\star}\subset[n_{1}]:=\\{1,\dots,n_{1}\\}$.

###  2.1 How to Combine Sparse and Low-Rank Weighting

As discussed in Section 1, the challenge in designing a reconstruction method
for Eq. 5 lies in simultaneously leveraging both structures of
$\mathbf{X}_{\star}$ in order to achieve the optimal sample complexity of
$m\approx r(s+n_{2})$. To this end, we propose a novel IRLS-based approach in
Algorithm 1. The key ingredient for computing the $(k+1)$-st iterate
$\mathbf{X}^{(k+1)}\in\mathbb{R}^{n_{1}\times n_{2}}$ in Algorithm 1 is the
multi-structural weight operator
$W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}:\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}^{n_{1}\times n_{2}}$ of Definition 2.1, which depends on
the current iterate $\mathbf{X}^{(k)}$.

######  Definition 2.1.

For $\delta_{k},\varepsilon_{k}>0$,
$\sigma_{i}^{k}:=\sigma_{i}(\mathbf{X}^{(k)})$, and
$\mathbf{X}^{(k)}\in\mathbb{R}^{n_{1}\times n_{2}}$, let (6)
$r_{k}:=|\\{i\in[n]:\sigma_{i}^{(k)}>\varepsilon_{k}\\}|$ denote the number of
singular values $\sigma^{(k)}=\big{(}\sigma_{i}^{(k)}\big{)}_{i=1}^{r_{k}}$ of
$\mathbf{X}^{(k)}$ larger than $\varepsilon_{k}$ and furthermore be (7)
$s_{k}:=\big{|}\big{\\{}i\in[n_{1}]:\|\mathbf{X}_{i,:}^{(k)}\|_{2}>\delta_{k}\big{\\}}\big{|}$
the number of rows of $\mathbf{X}^{(k)}$ with $\ell_{2}$-norm larger than
$\delta_{k}$. Define the $r_{k}$ left and right singular vectors of
$\mathbf{X}^{(k)}$ as columns of $\mathbf{U}\in\ R^{n_{1}\times r_{k}}$ and
$\mathbf{V}\in\mathbb{R}^{n_{2}\times r_{k}}$, respectively, corresponding to
the leading singular values
$\boldsymbol{\sigma}^{(k)}=\big{(}\sigma_{i}^{(k)}\big{)}_{i=1}^{r_{k}}:=\left(\sigma_{i}(\mathbf{X}^{(k)})\right)_{i=1}^{r_{k}}$.
We define the _weight operator_
$W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}:\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}^{n_{1}\times n_{2}}$ at iteration $k$ of Algorithm 1 as
(8) $\displaystyle
W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}(\mathbf{Z})=W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{Z})+\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbf{Z}$
where
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}:\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}^{n_{1}\times n_{2}}$ is its low-rank promoting part (9)
$\begin{split}W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{Z})=\begin{bmatrix}\mathbf{U}&\mathbf{U}_{\perp}\end{bmatrix}\boldsymbol{\Sigma}_{\varepsilon_{k}}^{-1}\begin{bmatrix}\mathbf{U}^{*}\\\
\mathbf{U}_{\perp}^{*}\end{bmatrix}\mathbf{Z}\begin{bmatrix}\mathbf{V}&\mathbf{V}_{\perp}\end{bmatrix}\boldsymbol{\Sigma}_{\varepsilon_{k}}^{-1}\begin{bmatrix}\mathbf{V}^{*}\\\
\mathbf{V}_{\perp}^{*}\end{bmatrix}\end{split}$ with
$\boldsymbol{\Sigma}_{\varepsilon_{k}}=\max(\sigma_{i}^{(k)}/\varepsilon_{k},1)$
and
$\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\in\mathbb{R}^{n_{1}\times
n_{1}}$ is its sparsity-promoting part, which is diagonal with (10)
$\displaystyle\left(\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\right)_{ii}=\max\left(\big{\|}(\mathbf{X}^{(k)})_{i,:}\big{\|}_{2}^{2}/\delta_{k}^{2},1\right)^{-1},\qquad\text{for
all $i\in[n_{1}]$}.$ The matrices $\mathbf{U}_{\perp}$ and
$\mathbf{V}_{\perp}$ are arbitrary complementary orthogonal bases for
$\mathbf{U}$ and $\mathbf{V}$ that do do not need to be computed in Algorithm
1.

Algorithm 1 IRLS for simultaneously low-rank rand row-sparse matrices

1: Input: Linear operator $\mathcal{A}\colon\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}^{m}$, data $\mathbf{y}\in\mathbb{R}^{m}$, rank and
sparsity estimates $\widetilde{r}$ and $\widetilde{s}$.

2: Initialize $k=0$,
$W_{\mathbf{X}^{(0)},\varepsilon_{0},\delta_{0}}=\mathbf{Id}$ and set
$\delta_{k},\varepsilon_{k}=\infty$.

3: for $k=1$ to $K$ do

4: Weighted Least Squares: Update iterate $\mathbf{X}^{(k)}$ by  (11)
$\mathbf{X}^{(k)}=\operatorname*{arg\,min}\limits_{\mathbf{X}:\mathcal{A}(\mathbf{X})=\mathbf{y}}\langle\mathbf{X},W_{\mathbf{X}^{(k-1)},\varepsilon_{k-1},\delta_{k-1}}(\mathbf{X})\rangle.$

5: Update Smoothing : Compute $\widetilde{r}+1$-st singular value
$\sigma_{\widetilde{r}+1}^{(k)}:=\sigma_{\widetilde{r}+1}(\mathbf{X}^{(k)})$
and $(\widetilde{s}+1)$-st largest row $\ell_{2}$-norm
$\rho_{\widetilde{s}+1}(\mathbf{X}^{(k)})$ of $\mathbf{X}^{(k)}$, update  (12)
$\varepsilon_{k}=\min\left(\varepsilon_{k-1},\sigma_{\widetilde{s}+1}^{(k)}\right),\quad\delta_{k}=\min\left(\delta_{k-1},\rho_{\widetilde{s}+1}(\mathbf{X}^{(k)})\right).$

6: Update Weight Operator: For $r_{k}$ and $s_{k}$ as in Eq. 6 and Eq. 7,

  * •

compute first $r_{k}$ singular triplets $\sigma^{(k)}\in\mathbb{R}^{r_{k}}$,
$\mathbf{U}\in\mathbb{R}^{n_{1}\times r_{k}}$ and
$\mathbf{V}\in\mathbb{R}^{n_{2}\times r_{k}}$,

  * •

compute $W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}$ in Eq. 8 via
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}$ in Eq. 9 and
$\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}$ in Eq. 10.

7: end for

8: Output: $\mathbf{X}^{(K)}$.

The high-level idea of Algorithm 1, as for other IRLS methods, is to minimize
quadratic functionals, which we call
$\mathcal{Q}_{\textit{lr},\varepsilon_{k}}(\,\cdot\,|\mathbf{X}^{(k)}):\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}$ and
$\mathcal{Q}_{\textit{sp},\delta_{k}}(\,\cdot\,|\mathbf{X}^{(k)}):\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}$ and define them by

(13)
$\displaystyle\begin{split}\mathcal{Q}_{\textit{lr},\varepsilon_{k}}(\mathbf{Z}|\mathbf{X}^{(k)})\\!&:=\\!\mathcal{F}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)})+\langle\nabla\mathcal{F}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)}),\mathbf{Z}-\mathbf{X}^{(k)}\rangle+\frac{1}{2}\langle\mathbf{Z}-\mathbf{X}^{(k)},W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{Z}-\mathbf{X}^{(k)})\rangle,\\\
\mathcal{Q}_{\textit{sp},\delta_{k}}(\mathbf{Z}|\mathbf{X}^{(k)})&:=\mathcal{F}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k)})+\langle\nabla\mathcal{F}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k)}),\mathbf{Z}-\mathbf{X}^{(k)}\rangle+\frac{1}{2}\langle\mathbf{Z}-\mathbf{X}^{(k)},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}(\mathbf{Z}-\mathbf{X}^{(k)})\rangle,\end{split}$

that majorize $\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\cdot)$ (see Theorem
2.5 below) for any iteration $k$. This minimization leads to the weighted
least squares problem Eq. 11 in Algorithm 1. This step can be implemented by
standard numerical linear algebra (see the supplementary material for a
discussion of its computational complexity). As a second ingredient of the
method, the smoothing parameters $\varepsilon_{k}$ and $\delta_{k}$ of
$\mathcal{F}_{\varepsilon_{k},\delta_{k}}$ are updated (i.e., decreased) in
step Eq. 12 before the weight operator is updated according to the current
iterate information. For the weight operator update, it is only necessary to
compute row norms and leading singular triplets of $\mathbf{X}^{(k)}$.

######  Remark 2.2.

The particular form of the low-rank promoting part of the weight operator
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}$ in Eq. 9 is due to
[Küm19, KMV21] and captures optimally spectral information both in the column
and row space, unlike prior work on low-rank IRLS [FRW11, MF12], while
retaining the property that the induced quadratic model
$\mathcal{Q}_{\textit{lr},\varepsilon_{k}}(\,\cdot\,|\mathbf{X}^{(k)})$
majorizes $\mathcal{F}_{\textit{lr},\varepsilon_{\varepsilon_{k}}}(\cdot)$
(see proof of Theorem 2.5). This choice is critical to enable a fast local
rate as established in Theorem 2.4.

###  2.2 Local Quadratic Convergence of IRLS

Our first main result states that Algorithm 1 exhibits quadratic convergence
in a local neighborhood of $\mathbf{X}_{\star}$, a property Algorithm 1 shares
with several methods from the IRLS family. We only need to assume that
$\mathcal{A}$ acts almost isometrically on the set $\mathcal{M}_{r,s}$.

######  Definition 2.3.

We say that a linear operator $\mathcal{A}\colon\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}^{m}$ satisfies the rank-$r$ and row-$s$-sparse restricted
isometry property (or $(r,s)$-RIP) with RIP-constant $\delta\in(0,1)$ if

$\displaystyle(1-\delta)\left\|\mathbf{Z}\right\|_{F}^{2}\leq\left\|\mathcal{A}(\mathbf{Z})\right\|_{2}^{2}\leq(1+\delta)\left\|\mathbf{Z}\right\|_{F}^{2},$

for all $\mathbf{Z}\in\mathcal{M}_{r,s}$.

It is worth highlighting that Gaussian operators satisfy the above RIP with
high-probability if $m\geq cr(s+n_{2})\log(en_{1}/s)$, for some absolute
constant $c>0$, see for instance [LWB18]. Up to log-factors, this is at the
information theoretic limit which we discussed in the beginning. The
convergence result for Algorithm 1 now reads as follows.

######  Theorem 2.4 (Local Quadratic Convergence).

Let $\mathbf{X}_{\star}\in\mathcal{M}_{r,s}$ be a fixed ground-truth matrix
that is $s$-row-sparse and of rank $r$. Let linear observations
$\mathbf{y}=\mathcal{A}(\mathbf{X}_{\star})$ be given and assume that
$\mathcal{A}$ has the $(r,s)$-RIP with $\delta\in(0,1)$. Assume that the
$k$-th iterate $\mathbf{X}^{(k)}$ of Algorithm 1 with $\widetilde{r}=r$ and
$\widetilde{s}=s$ updates the smoothing parameters in Eq. 12 such that one of
the statements $\varepsilon_{k}=\sigma_{r+1}(\mathbf{X}^{(k)})$ or
$\delta_{k}=\rho_{s+1}(\mathbf{X}^{(k)})$ is true, and that $r_{k}\geq r$ and
$s_{k}\geq s$.

If

(14)
$\displaystyle\varepsilon_{k}\leq\frac{1}{48}\sigma_{r}(\mathbf{X}_{\star})$

and if $\mathbf{X}^{(k)}$ satisfies

(15) $\displaystyle\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|$
$\displaystyle\leq\frac{1}{48\sqrt{n}c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{3}}\min\left\\{\frac{\sigma_{r}(\mathbf{X}_{\star})}{r},\frac{\rho_{s}(\mathbf{X}_{\star})}{s}\right\\}$

where $c_{\left\|\mathcal{A}\right\|_{2\to
2}}=\sqrt{1+\tfrac{\left\|\mathcal{A}\right\|_{2\to 2}^{2}}{(1-\delta)}}$ and
$n=\min\\{n_{1},n_{2}\\}$, then the local convergence rate is quadratic in the
sense that

$\displaystyle\|\mathbf{X}^{(k+1)}-\mathbf{X}_{\star}\|\leq\min\\{\mu\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|^{2},0.9\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|\\},$

for

(16) $\displaystyle\mu=4.179c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{2}\Big{(}\frac{5r}{\sigma_{r}(\mathbf{X}_{\star})}+\frac{2s}{\rho_{s}(\mathbf{X}_{\star})}\Big{)},$

and $\mathbf{X}^{(k+\ell)}\overset{\ell\to\infty}{\to}\mathbf{X}_{\star}$.

The proof of Theorem 2.4 is presented in the supplementary material. To the
best of our knowledge, so far no other method exists for recovering
simultaneously sparse and low-rank matrices that exhibits local quadratic
convergence. In particular, the state-of-the-art competitor methods [LWB18,
FGJR20, Mal21, EKPU22] reach a local linear error decay at best.
On the other hand, Eq. 15 is rather pessimistic since for Gaussian
$\mathcal{A}$ the constant $c_{\left\|\mathcal{A}\right\|_{2\to 2}}$ scales
like $\sqrt{(n_{1}n_{2})/m}$, which means that the right-hand side of Eq. 15
behaves like $m^{3/2}/(n(n_{1}n_{2})^{3/2})$, whereas we observe quadratic
convergence in experiments within an empirically much larger convergence
radius. Closing this gap between theory and practical performance is future
work.
It is noteworthy that the theory in [LWB18] — to our knowledge the only other
related work explicitly characterizing the convergence radius — holds on a
neighborhood of $\mathbf{X}_{\star}$ that is independent of the ambient
dimension. The authors of [LWB18] however assume that the RIP-constant decays
with the conditioning number $\kappa$ of $\mathcal{A}$, a quantity that might
be large in applications. Hence, ignoring log-factors the sufficient number of
measurements in [LWB18] scales like $m=\Omega(\kappa^{2}r(s_{1}+n_{2}))$. In
contrast, Theorem 2.4 works for any RIP-constant less than one which means for
$m=\Omega(r(s_{1}+n_{2}))$.

###  2.3 IRLS as Quadratic Majorize-Minimize Algorithm

With Theorem 2.4, we have provided a local convergence theorem that quantifies
the behavior of Algorithm 1 in a small neighbourhood of the simultaneously
row-sparse and low-rank ground-truth
$\mathbf{X}_{\star}\in\mathbb{R}^{n_{1}\times n_{2}}$. The result is based on
sufficient regularity of the measurement operator $\mathcal{A}$, which in turn
is satisfied with high probability if $\mathcal{A}$ consists of sufficiently
generic random linear observations that concentrate around their mean.
In this section we establish that, for _any_ measurement operator
$\mathcal{A}$, Algorithm 1 can be interpreted within the framework of
iteratively reweighted least squares (IRLS) algorithms [DDFG10, MF12, PKV22],
which implies a strong connection to the minimization of a suitable smoothened
objective function. In our case, the objective
$\mathcal{F}_{\varepsilon,\delta}$ in Eq. 4 is a linear combination of sum-of-
logarithms terms penalizing both non-zero singular values [FHB03, MF12,
CESV13] as well as non-zero rows of a matrix $\mathbf{X}$ [KASL21].
We show in Theorem 2.5 below that the IRLS algorithm Algorithm 1 studied in
this paper is based on minimizing at each iteration quadratic models that
majorize $\mathcal{F}_{\varepsilon,\delta}$, and furthermore, that the
iterates $(\mathbf{X}^{(k)})_{k\geq 1}$ of Algorithm 1 define a non-increasing
sequence
$\left(\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k)})\right)_{k\geq
1}$ with respect to the objective $\mathcal{F}_{\varepsilon,\delta}$ of Eq. 4.
The proof combines the fact that for fixed smoothing parameters
$\varepsilon_{k}$ and $\delta_{k}$, the weighted least squares and weight
update steps Algorithm 1 can be interpreted as a step of a Majorize-Minimize
algorithm [SBP16, Lan16], with a decrease in the underlying objective Eq. 4
for updated smoothing parameters.

######  Theorem 2.5.

Let $\mathbf{y}\in\mathbb{R}^{m}$, let the linear operator
$\mathcal{A}\colon\mathbb{R}^{n_{1}\times n_{2}}\to\mathbb{R}^{m}$ be
arbitrary. If $(\mathbf{X}^{(k)})_{k\geq 1}$ is a sequence of iterates of
Algorithm 1 and $(\delta_{k})_{k\geq 1}$ and $(\varepsilon_{k})_{k\geq 1}$ are
the sequences of smoothing parameters as defined therein, then the following
statements hold.

  1. 1.

The quadratic model functions
$\mathcal{Q}_{\textit{lr},\varepsilon_{k}}(\,\cdot\,|\mathbf{X}^{(k)})$ and
$\mathcal{Q}_{\textit{sp},\delta_{k}}(\,\cdot\,|\mathbf{X}^{(k)})$ defined in
Eq. 13 globally majorize the $(\varepsilon_{k},\delta_{k})-$smoothed
logarithmic surrogate objective $\mathcal{F}_{\varepsilon_{k},\delta_{k}}$,
i.e., for _any_ $\mathbf{Z}\in\mathbb{R}^{n_{1}\times n_{2}}$, it holds that

(17)
$\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{Z})\leq\mathcal{Q}_{\textit{lr},\varepsilon_{k}}(\mathbf{Z}|\mathbf{X}^{(k)})+\mathcal{Q}_{\textit{sp},\delta_{k}}(\mathbf{Z}|\mathbf{X}^{(k)}).$

  2. 2.

The sequence
$\left(\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k)})\right)_{k\geq
1}$ is non-increasing.

  3. 3.

If $\overline{\varepsilon}:=\lim_{k\to\infty}\varepsilon_{k}>0$ and
$\overline{\delta}:=\lim_{k\to\infty}\delta_{k}>0$, then
$\lim_{k\to\infty}\|\mathbf{X}^{(k)}-\mathbf{X}^{(k+1)}\|_{F}=0$. Furthermore,
in this case, every accumulation point of $(\mathbf{X}^{(k)})_{k\geq 1}$ is a
stationary point of

$\min\limits_{\mathbf{X}:\mathcal{A}(\mathbf{X})=\mathbf{y}}\mathcal{F}_{\overline{\varepsilon},\overline{\delta}}(\mathbf{X}).$

The proof of Theorem 2.5 is presented in the supplementary material.

##  3 Numerical Evaluation

In this section, we explore the empirical performance of IRLS in view of the
theoretical results of Theorems 2.4 and 2.5, and compare its ability to
recover simultaneous low-rank and row-sparse data matrices with the state-of-
the-art methods Sparse Power Factorization (SPF) [LWB18] and Riemannian
adaptive iterative hard thresholding (RiemAdaIHT)
[EKPU22], which are among the methods with the best empirical performance
reported in the literature. The method ATLAS [FMN21] and its successor [Mal21]
are not used in our empirical studies since they are tailored to robust
recovery and yield suboptimal performance when seeking high-precision
reconstruction in low noise scenarios. We use spectral initialization for SPF
and RiemAdaIHT. The weight operator of IRLS is initialized by the identity as
described in Algorithm 1, solving an unweighted least squares problem in the
first iteration.

#### Performance in Low-Measurement Regime

Figures 1 and 2 show the empirical probability of successful recovery when
recovering $s$-row sparse ground-truths
$\mathbf{X}_{\star}\in\mathbb{R}^{256\times 40}$ of rank $r=1$ (resp. $r=5$)
from Gaussian measurements under oracle knowledge on $r$ and $s$. The results
are averaged over 64 random trials. As both figures illustrate, the region of
success of IRLS comes closest to the information theoretic limit of
$r(s+n_{2}-r)$ which is highlighted by a red line, requiring a significantly
lower oversampling factor than the baseline methods.

$10$$20$$30$$40$$50$$100$$150$$200$Row sparsity $s$Nr. of measurements $m$

$10$$20$$30$$40$$50$$100$$150$$200$

$10$$20$$30$$40$$50$$100$$150$$200$$0$$0.2$$0.4$$0.6$$0.8$$1$

Figure 1: Left column: RiemAdaIHT, center: SPF, right: IRLS. Phase transition
experiments with $n_{1}=256$, $n_{2}=40$, $r=1$, Gaussian measurements.
Algorithmic hyperparameters informed by model order knowledge (i.e.,
$\widetilde{r}=r$ and $\widetilde{s}=s$ for IRLS). White corresponds to
empirical success rate of $1$, black to $0$.

$20$$40$$60$$400$$600$$800$$1{,}000$Number of measurements $m$

$20$$40$$60$$400$$600$$800$$1{,}000$

$20$$40$$60$$400$$600$$800$$1{,}000$

$20$$40$$60$$400$$600$$800$$1{,}000$Row sparsity $s$

$20$$40$$60$$400$$600$$800$$1{,}000$

$20$$40$$60$$400$$600$$800$$1{,}000$

Figure 2: Left column: RiemAdaIHT, center: SPF, right: IRLS. First row: As in
Figure 1, but for data matrix $\mathbf{X}_{\star}$ of rank $r=5$. Second row:
As first row, but hyper-parameters $r$ and $s$ are overestimated as
$\widetilde{r}=2r=10$, $\widetilde{s}=\lfloor 1.5s\rfloor$.

#### Sensitivity to Parameter Choice

In applications of our setting, the quantities $r$ and $s$ might be unknown or
difficult to estimate. In the second row of Figure 2, we repeat the experiment
of the first row (rank-$5$ ground truth), but run the algorithms with rank and
sparsity estimates of $\widetilde{r}=2r$ and $\widetilde{s}=\lfloor
1.5s\rfloor$. Whereas all considered methods suffer a deterioration of
performance, we observe that IRLS deteriorates relatively the least by a large
margin. Furthermore, we observe that even if IRLS does not recovery
$\mathbf{X}_{\star}$, it converges typically to a matrix that is still low-
rank and row-sparse (with larger $r$ and $s$) satisfying the data constraint,
while the other methods fail to convergence to such a matrix.

#### Convergence Behavior

Finally, we examine the convergence rate of the iterates to validate the
theoretical prediction of Theorem 2.4 in the setting of Figure 2. Figure 3
depicts in log-scale the approximation error over the iterates of SPF,
RiemAdaIHT, and IRLS. We observe that the IRLS indeed exhibits empirical
quadratic convergence within a few iterations (around $10$), whereas the other
methods clearly only exhibit linear convergence. The experiment further
suggests that the rather pessimistic size of the convergence radius
established by Theorem 2.4 could possibly be improved by future
investigations.

$0$$7$$14$$21$$28$$35$$42$$49$$56$$63$$70$$77$$84$$91$$98$$105$$112$$10^{-11}$$10^{-8}$$10^{-5}$$10^{-2}$Iteration
$k$Rel. Frobenius errorIRLS-LowRankRowSparseSparse Power Factorization
[LWB18]Riemannian Adaptive IHT [EKPU22] Figure 3: Comparison of convergence
rate. Setting as in Figure 2 with $s=40$ and $m=1125$.

##  4 Conclusion, Limitations and Future Work

### Conclusion

In this paper, we adapted the IRLS framework to the problem of recovering
simultaneously structured matrices from linear observations focusing on the
special case of row sparsity and low-rankness. Our convergence guarantee
Theorem 2.4 is hereby the first one for any method minimizing combinations of
structural surrogate objectives that holds in the information-theoretic near-
optimal regime and exhibits local quadratic convergence. The numerical
experiments we conducted for synthetic data suggest that, due to its weak
dependence on the choice of hyperparameters, IRLS in the form of Algorithm 1
can be a practical method for identifying simultaneously structured data even
in difficult problem instances.

### Limitations and Future Work

As in the case of established IRLS methods that optimize non-convex surrogate
objectives representing a single structure [DDFG10, KS18, KMV21, PKV22], the
radius of guaranteed quadratic convergence in Theorem 2.4 is the most
restrictive assumption. Beyond the interpretation in terms of surrogate
minimization as presented in Theorem 2.5, which holds without any assumptions
on the initialization, our method shares a lack of a global convergence
guarantees with other non-convex IRLS algorithms [MF12, KMV21, PKV22].

The generalization and application of our framework to combinations of
structures beyond rank and row- (or column-)sparsity lies outside the scope of
the present paper, but could involve subspace-structured low-rankness [FHB03,
CC14, YHC18, SOF20] or analysis sparsity [EA06, RYF19]. A generalization of
the presented IRLS framework to higher-order objects such as low-rank tensors
is of future interest as convexifications of structure-promoting objectives
face similar challenges [MHWG14, OJF+15, YKH19] in this case.

In parameter-efficient deep learning, both sparse [FC19, HABN+21, FZS22,
VTCZ+22] and low-rank [WXW+17, WLK+20, SZK+22] weight parameter models have
gained considerable attention due to the challenges of training and storing,
e.g., in large transformer-based models [VSP+17, BMR+20]. It will be of
interest to study whether in this non-linear case IRLS-like preconditioning of
the parameter space can find network weights that are simultaneously sparse
and low-rank, and could potentially lead to further increases in efficiency.

## References

  * [AH15] K. Aftab and R. Hartley.  Convergence of iteratively re-weighted least squares to robust m-estimators.  In IEEE Winter Conference on Applications of Computer Vision (WACV 2015), pages 480–487, 2015.
  * [APS19] D. Adil, R. Peng, and S. Sachdeva.  Fast, Provably convergent IRLS Algorithm for p-norm Linear Regression.  Advances in Neural Information Processing Systems (NeurIPS 2019), 2019.
  * [ARR14] A. Ahmed, B. Recht, and J. Romberg.  Blind deconvolution using convex programming.  IEEE Transactions on Information Theory, 60(3):1711–1732, 2014\.
  * [BB18] M. Benning and M. Burger.  Modern regularization methods for inverse problems.  Acta Numerica, 27:1–111, 2018.
  * [BBDM21] M. Bertero, P. Boccacci, and C. De Mol.  Introduction to inverse problems in imaging.  CRC press, 2021.
  * [BDE09] A. M. Bruckstein, D. L. Donoho, and M. Elad.  From sparse solutions of systems of equations to sparse modeling of signals and images.  SIAM review, 51(1):34–81, 2009.
  * [BMR+20] T. Brown, B. Mann, N. Ryder, M. Subbiah, J. D. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell, et al.  Language models are few-shot learners.  Advances in Neural Information Processing Systems (NeurIPS 2020), 33:1877–1901, 2020.
  * [Bou23] N. Boumal.  An introduction to optimization on smooth manifolds.  Cambridge University Press, 2023.
  * [BR16] S. Bahmani and J. Romberg.  Near-optimal estimation of simultaneously sparse and low-rank matrices from nested linear measurements.  Information and Inference: A Journal of the IMA, 5(3):331–351, 2016\.
  * [BS15] A. Beck and S. Sabach.  Weiszfeld’s method: Old and new results.  Journal of Optimization Theory and Applications, 164(1):1–40, 2015\.
  * [CC14] Y. Chen and Y. Chi.  Robust Spectral Compressed Sensing via Structured Matrix Completion.  IEEE Transactions on Information Theory, 60(10):6576–6601, 2014\.
  * [CE17] P. Campisi and K. Egiazarian.  Blind Image Deconvolution: Theory and Applications.  CRC Press, 2017.
  * [CESV13] E. J. Candès, Y. Eldar, T. Strohmer, and V. Voroninski.  Phase Retrieval via Matrix Completion.  SIAM J. Imag. Sci., 6(1):199–225, 2013.
  * [CG17] A. Chatterjee and V. M. Govindu.  Robust relative rotation averaging.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 40(4):958–972, 2017.
  * [CH12] L. Chen and J. Z. Huang.  Sparse reduced-rank regression for simultaneous dimension reduction and variable selection.  Journal of the American Statistical Association, 107(500):1533–1545, 2012.
  * [Cha07] R. Chartrand.  Exact reconstruction of sparse signals via nonconvex minimization.  IEEE Signal Processing Letters, 14(10):707–710, 2007.
  * [Che18] W. Chen.  Simultaneously sparse and low-rank matrix reconstruction via nonconvex and nonseparable regularization.  IEEE Transactions on Signal Processing, 66(20):5313–5323, 2018\.
  * [CHHL14] C. Chen, J. Huang, L. He, and H. Li.  Preconditioning for accelerated iteratively reweighted least squares in structured sparsity reconstruction.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR 2014), pages 2713–2720, 2014.
  * [CHLH18] C. Chen, L. He, H. Li, and J. Huang.  Fast iteratively reweighted least squares algorithms for analysis-based sparse reconstruction.  Medical Image Analysis, 49:141–152, 2018.
  * [CLM16] T. T. Cai, X. Li, and Z. Ma.  Optimal rates of convergence for noisy sparse phase retrieval via thresholded Wirtinger flow.  The Annals of Statistics, pages 2221–2251, 2016.
  * [CP11] E. J. Candè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, 2011.
  * [CRPW12] V. Chandrasekaran, B. Recht, P. A. Parrilo, and A. S. Willsky.  The Convex Geometry of Linear Inverse Problems.  Foundations of Computational Mathematics, 12(6):805–849, 2012.
  * [CY08] R. Chartrand and W. Yin.  Iteratively reweighted algorithms for compressive sensing.  In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3869–3872, 2008.
  * [CZ18] T. T. Cai and A. Zhang.  Rate-optimal perturbation bounds for singular subspaces with applications to high-dimensional statistics.  The Annals of Statistics, 46(1):60–89, 2018.
  * [Dax10] A. Dax.  On extremum properties of orthogonal quotients matrices.  Linear Algebra and its Applications, 432(5):1234–1257, 2010.
  * [DDFG10] I. Daubechies, R. DeVore, M. Fornasier, and C. Güntürk.  Iteratively Reweighted Least Squares Minimization for Sparse Recovery.  Commun. Pure Appl. Math., 63:1–38, 2010.
  * [dGJL05] A. d’Aspremont, L. E. Ghaoui, M. I. Jordan, and G. R. Lanckriet.  A direct formulation for sparse PCA using semidefinite programming.  In Advances in Neural Information Processing Systems (NIPS 2005), pages 41–48, 2005.
  * [Don06] D. L. Donoho.  Compressed sensing.  IEEE Transactions on Information Theory, 52(4):1289–1306, 2006\.
  * [EA06] M. Elad and M. Aharon.  Image Denoising Via Sparse and Redundant Representations Over Learned Dictionaries.  IEEE Transactions on Image Processing, 15(12):3736–3745, 2006.
  * [EKPU22] H. Eisenmann, F. Krahmer, M. Pfeffer, and A. Uschmajew.  Riemannian thresholding methods for row-sparse and low-rank matrix recovery.  Numerical Algorithms, pages 1–25, 2022.
  * [FC19] J. Frankle and M. Carbin.  The Lottery Ticket Hypothesis: Finding Sparse, Trainable Neural Networks.  In International Conference on Learning Representations, 2019.
  * [FGJR20] S. Foucart, R. Gribonval, L. Jacques, and H. Rauhut.  Jointly low-rank and bisparse recovery: Questions and partial answers.  Analysis and Applications, 18(01):25–48, 2020.
  * [FHB03] M. Fazel, H. Hindi, and S. P. Boyd.  Log-det heuristic for matrix rank minimization with applications to Hankel and Euclidean distance matrices.  In Proceedings of the 2003 American Control Conference, volume 3, pages 2156–2162. IEEE, 2003.
  * [FMN21] M. Fornasier, J. Maly, and V. Naumova.  Robust recovery of low-rank matrices with non-orthogonal sparse decomposition from incomplete measurements.  Applied Mathematics and Computation, 392, 2021.
  * [Fou11] S. Foucart.  Hard thresholding pursuit: an algorithm for compressive sensing.  SIAM Journal on Numerical Analysis, 49(6):2543–2563, 2011.
  * [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.  Birkhäuser Basel, 2013.
  * [FRW11] M. Fornasier, H. Rauhut, and R. Ward.  Low-rank matrix recovery via iteratively reweighted least squares minimization.  SIAM Journal on Optimization, 21(4):1614–1640, 2011.
  * [FZS22] W. Fedus, B. Zoph, and N. Shazeer.  Switch transformers: Scaling to trillion parameter models with simple and efficient sparsity.  The Journal of Machine Learning Research, 23(1):5232–5270, 2022\.
  * [GR97] I. Gorodnitsky and B. Rao.  Sparse signal reconstruction from limited data using FOCUSS: A re-weighted minimum norm algorithm.  IEEE Trans. Signal Process., 45(3):600–616, 1997.
  * [GV12] M. Golbabaee and P. Vandergheynst.  Hyperspectral Image Compressed Sensing via Low-Rank and Joint-Sparse Matrix Recovery.  Proceedings of ICASSP 2012, 2012.
  * [HABN+21] T. Hoefler, D. Alistarh, T. Ben-Nun, N. Dryden, and A. Peste.  Sparsity in deep learning: Pruning and growth for efficient inference and training in neural networks.  The Journal of Machine Learning Research, 22(1):10882–11005, 2021\.
  * [HV19] B. D. Haeffele and R. Vidal.  Structured low-rank matrix factorization: Global optimality, algorithms, and applications.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 42(6):1468–1482, 2019.
  * [HW77] P. W. Holland and R. E. Welsch.  Robust regression using iteratively reweighted least-squares.  Communications in Statistics - Theory and Methods, 6(9):813–827, 1977.
  * [JC93] S. M. Jefferies and J. C. Christou.  Restoration of astronomical images by iterative blind deconvolution.  The Astrophysical Journal, 415:862, 1993.
  * [JH17] G. Jagatap and C. Hegde.  Fast, sample-efficient algorithms for structured phase retrieval.  Advances in Neural Information Processing Systems (NeurIPS 2017), 30, 2017.
  * [JNS13] P. Jain, P. Netrapalli, and S. Sanghavi.  Low-rank matrix completion using alternating minimization.  Proceedings of the forty-fifth Annual ACM Symposium on Theory of Computing, pages 665–674, 2013.
  * [JOH17] K. Jaganathan, S. Oymak, and B. Hassibi.  Sparse phase retrieval: Uniqueness guarantees and recovery algorithms.  IEEE Transactions on Signal Processing, 65(9):2402–2410, 2017.
  * [KASL21] C. Ke, M. Ahn, S. Shin, and Y. Lou.  Iteratively Reweighted Group Lasso Based on Log-Composite Regularization.  SIAM Journal on Scientific Computing, 43(5):S655–S678, 2021.
  * [KMV21] C. Kümmerle and C. Mayrink Verdun.  A Scalable Second Order Method for Ill-Conditioned Matrix Completion from Few Samples.  In International Conference on Machine Learning (ICML 2021), 2021\.
  * [KMVS21] C. Kümmerle, C. Mayrink Verdun, and D. Stöger.  Iteratively Reweighted Least Squares for Basis Pursuit with Global Linear Convergence Rate.  Advances in Neural Information Processing Systems (NeurIPS 2021), 2021.
  * [KS18] C. Kümmerle and J. Sigl.  Harmonic mean iteratively reweighted least squares for low-rank matrix recovery.  The Journal of Machine Learning Research, 19(1):1815–1863, 2018\.
  * [KSJ19] M. Kliesch, S. J. Szarek, and P. Jung.  Simultaneous Structures in Convex Signal Recovery—Revisiting the Convex Combination of Norms.  Frontiers in Applied Mathematics and Statistics, 5:23, 2019.
  * [KST95] M. V. Klibanov, P. E. Sacks, and A. V. Tikhonravov.  The phase retrieval problem.  Inverse problems, 11(1):1, 1995.
  * [Küm19] C. Kümmerle.  Understanding and Enhancing Data Recovery Algorithms: From Noise-Blind Sparse Recovery to Reweighted Methods for Low-Rank Matrix Optimization.  Dissertation, Technical University of Munich, 2019.  Available at https://mediatum.ub.tum.de/doc/1521436/1521436.pdf.
  * [Lan16] K. Lange.  MM optimization algorithms.  SIAM, 2016.
  * [LC22] S. H. Lee and J. Civera.  HARA: A hierarchical approach for robust rotation averaging.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR 2022), pages 15777–15786, 2022.
  * [LK23] S. Lefkimmiatis and I. S. Koshelev.  Learning Sparse and Low-Rank Priors for Image Recovery via Iterative Reweighted Least Squares Minimization.  In The Eleventh International Conference on Learning Representations, 2023.
  * [LLJB16] K. Lee, Y. Li, M. Junge, and Y. Bresler.  Blind recovery of sparse signals from subsampled convolution.  IEEE Transactions on Information Theory, 63(2):802–821, 2016.
  * [LLSW19] X. Li, S. Ling, T. Strohmer, and K. Wei.  Rapid, robust, and reliable blind deconvolution via nonconvex optimization.  Applied and Computational Harmonic Analysis, 47(3):893–934, 2019\.
  * [LS05] A. S. Lewis and H. S. Sendov.  Nonsmooth Analysis of Singular Values. Part I: Theory.  Set-Valued Analysis, 13(3):213–241, 2005.
  * [LW20] H. Lyu and R. Wang.  An exact sin $\Theta$ formula for matrix perturbation analysis and its applications.  arXiv preprint arXiv:2011.07669, 2020.
  * [LWB18] K. Lee, Y. Wu, and Y. Bresler.  Near-Optimal Compressed Sensing of a Class of Sparse Low-Rank Matrices Via Sparse Power Factorization.  IEEE Transactions on Information Theory, 64(3):1666–1698, 2018\.
  * [LXY13] M.-J. Lai, Y. Xu, and W. Yin.  Improved iteratively reweighted least squares for unconstrained smoothed $\ell_{q}$-minimization.  SIAM Journal on Numerical Analysis, 51(2):927–957, 2013.
  * [Mal21] J. Maly.  Robust Sensing of Low-Rank Matrices with Non-Orthogonal Sparse Decomposition.  arXiv preprint arXiv:2103.05523, 2021.
  * [MBP+14] J. Mairal, F. Bach, J. Ponce, et al.  Sparse modeling for image and vision processing.  Foundations and Trends® in Computer Graphics and Vision, 8(2-3):85–283, 2014.
  * [MF12] K. Mohan and M. Fazel.  Iterative Reweighted Algorithms for Matrix Rank Minimization.  J. Mach. Learn. Res., 13(1):3441–3473, 2012.
  * [MGJK19] B. Mukhoty, G. Gopakumar, P. Jain, and P. Kar.  Globally-convergent Iteratively Reweighted Least Squares for Robust Regression Problems.  In International Conference on Artificial Intelligence and Statistics, pages 313–322, 2019.
  * [MHWG14] C. Mu, B. Huang, J. Wright, and D. Goldfarb.  Square Deal: Lower Bounds and Improved Relaxations for Tensor Recovery.  In Proceedings of the 31st International Conference on Machine Learning (ICML 2014), volume 32 of Proceedings of Machine Learning Research, pages 73–81, 2014.
  * [MI17] M. Magdon-Ismail.  NP-hardness and inapproximability of sparse PCA.  Information Processing Letters, 126:35–38, 2017.
  * [MWCC20] C. Ma, K. Wang, Y. Chi, and Y. Chen.  Implicit Regularization in Nonconvex Statistical Estimation: Gradient Descent Converges Linearly for Phase Retrieval, Matrix Completion, and Blind Deconvolution.  Foundations of Computational Mathematics, 2020.
  * [ODBP15] P. Ochs, A. Dosovitskiy, T. Brox, and T. Pock.  On iteratively reweighted algorithms for nonsmooth nonconvex optimization in computer vision.  SIAM Journal on Imaging Sciences, 8(1):331–372, 2015.
  * [OJF+15] S. Oymak, A. Jalali, M. Fazel, Y. C. Eldar, and B. Hassibi.  Simultaneously structured models with application to sparse and low-rank matrices.  IEEE Transactions on Information Theory, 61(5):2886–2908, 2015\.
  * [PKV22] L. Peng, C. Kümmerle, and R. Vidal.  Global Linear and Local Superlinear Convergence of IRLS for Non-Smooth Robust Regression.  In Advances in Neural Information Processing Systems (NeurIPS 2022), volume 35, pages 28972–28987, 2022.
  * [PKV23] L. Peng, C. Kümmerle, and R. Vidal.  On the Convergence of IRLS and Its Variants in Outlier-Robust Estimation.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR 2023), 2023.
  * [RFP10] B. Recht, M. Fazel, and P. A. Parrilo.  Guaranteed Minimum-Rank Solutions of Linear Matrix Equations via Nuclear Norm Minimization.  SIAM Review, 52(3):471–501, 2010.
  * [RSV12] E. Richard, P.-A. Savalle, and N. Vayatis.  Estimation of simultaneously sparse and low rank matrices.  In Proceedings of the 29th International Conference on Machine Learning (ICML 2012), pages 51–58, 2012.
  * [RYF19] S. Ravishankar, J. C. Ye, and J. A. Fessler.  Image reconstruction: From sparsity to data-adaptive methods and machine learning.  Proceedings of the IEEE, 108(1):86–109, 2019.
  * [SBP16] Y. Sun, P. Babu, and D. P. Palomar.  Majorization-minimization algorithms in signal processing, communications, and machine learning.  IEEE Transactions on Signal Processing, 65(3):794–816, 2016.
  * [SC21] L. Shi and Y. Chi.  Manifold gradient descent solves multi-channel sparse blind deconvolution provably and efficiently.  IEEE Transactions on Information Theory, 67(7):4784–4811, 2021\.
  * [SM17] M. Samejima and Y. Matsushita.  Fast general norm approximation via iteratively reweighted least squares.  In Computer Vision–ACCV 2016 Workshops: ACCV 2016 International Workshops, Taipei, Taiwan, November 20-24, 2016, pages 207–221. Springer, 2017\.
  * [SOF20] Y. Sun, S. Oymak, and M. Fazel.  Finite Sample System Identification: Optimal Rates and the Role of Regularization.  In Proceedings of the 2nd Conference on Learning for Dynamics and Control, volume 120, pages 16–25. PMLR, 10–11 Jun 2020.
  * [Sol19] M. Soltanolkotabi.  Structured signal recovery from quadratic measurements: Breaking sample complexity barriers via nonconvex optimization.  IEEE Transactions on Information Theory, 65(4):2374–2400, 2019\.
  * [SWL22] Y. Shi, C. M. Wyeth, and G. Lerman.  Robust Group Synchronization via Quadratic Programming.  In International Conference on Machine Learning (ICML 2022), pages 20095–20105. PMLR, 2022.
  * [SZK+22] S. Schotthöfer, E. Zangrando, J. Kusch, G. Ceruti, and F. Tudisco.  Low-rank lottery tickets: finding efficient low-rank neural networks via matrix differential equations.  In Advances in Neural Information Processing Systems (NeurIPS 2022), volume 35, pages 20051–20063, 2022.
  * [Tik63] A. N. Tikhonov.  Solution of incorrectly formulated problems and the regularization method.  Soviet Math., 4:1035–1038, 1963.
  * [TRB17] C. G. Tsinos, A. A. Rontogiannis, and K. Berberidis.  Distributed blind hyperspectral unmixing via joint sparsity and low-rank constrained non-negative matrix factorization.  IEEE Transactions on Computational Imaging, 3(2):160–174, 2017\.
  * [Van13] B. Vandereycken.  Low-Rank Matrix Completion by Riemannian Optimization.  SIAM J. Optim., 23(2), 2013.
  * [Vog02] C. R. Vogel.  Computational methods for inverse problems.  SIAM, 2002.
  * [VSP+17] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin.  Attention is all you need.  Advances in Neural Information Processing Systems (NeurIPS 2017), 30, 2017.
  * [VTCZ+22] M. Varma T, X. Chen, Z. Zhang, T. Chen, S. Venugopalan, and Z. Wang.  Sparse Winning Tickets are Data-Efficient Image Recognizers.  Advances in Neural Information Processing Systems (NeurIPS 2022), 35:4652–4666, 2022.
  * [Wei37] E. Weiszfeld.  Sur le point pour lequel la somme des distances de n points donnés est minimum.  Tohoku Mathematical Journal, First Series, 43:355–386, 1937.
  * [WLK+20] S. Wang, B. Z. Li, M. Khabsa, H. Fang, and H. Ma.  Linformer: Self-attention with linear complexity.  arXiv preprint arXiv:2006.04768, 2020.
  * [Woo50] M. A. Woodbury.  Inverting modified matrices.  Memorandum report, 42(106):336, 1950.
  * [WXW+17] W. Wen, C. Xu, C. Wu, Y. Wang, Y. Chen, and H. Li.  Coordinating filters for faster deep neural networks.  In Proceedings of the IEEE International Conference on Computer Vision (CVPR 2017), pages 658–666, 2017.
  * [YGK20] M. Yu, V. Gupta, and M. Kolar.  Recovery of simultaneous low rank and two-way sparse coefficient matrices, a nonconvex approach.  Electronic Journal of Statistics, 14:413–457, 2020.
  * [YHC18] J. C. Ye, Y. Han, and E. Cha.  Deep convolutional framelets: A general deep learning framework for inverse problems.  SIAM Journal on Imaging Sciences, 11(2):991–1048, 2018.
  * [YKH19] Q. Yao, J. T.-Y. Kwok, and B. Han.  Efficient nonconvex regularized tensor completion with structure-aware proximal iterations.  In International Conference on Machine Learning (ICML 2019), pages 7035–7044. PMLR, 2019.
  * [ZB15] Z. Zeinalkhani and A. H. Banihashemi.  Iterative Reweighted $\ell_{2}$/$\ell_{1}$ Recovery Algorithms for Compressed Sensing of Block Sparse Signals.  IEEE Transactions on Signal Processing, 63(17):4516–4531, 2015\.
  * [ZHL+22] Y. Zhang, L.-T. Huang, Y. Li, K. Zhang, and C. Yin.  Low-rank and sparse matrix recovery for hyperspectral image reconstruction using Bayesian learning.  Sensors, 22(1):343, 2022.
  * [ZHT06] H. Zou, T. Hastie, and R. Tibshirani.  Sparse Principal Component Analysis.  Journal of Computational and Graphical Statistics, 15(2):265–286, 2006.
  * [ZZL22] Z. Zhang, Z. Zhai, and L. Li.  Graph Refinement via Simultaneously Low-Rank and Sparse Approximation.  SIAM Journal on Scientific Computing, 44(3):A1525–A1553, 2022.

## Supplementary material for _Recovering Simultaneously Structured Data via
Non-Convex Iteratively Reweighted Least Squares_

This supplement is structured as follows.

  * •

Appendix A presents some details about the experimental setup as well as
additional numerical experiments.

  * •

Appendix B discusses the computational complexity of Algorithm 1.

  * •

Section C.1 presents the proof of Theorem 2.4.

  * •

Section C.2 presents the proof of Theorem 2.5.

  * •

Appendix D details some technical results that are used in Sections C.1 and
C.2.

##  Appendix A Experimental Setup and Supplementary Experiments

In this section, we elaborate on the detailed experimental setup that was used
in Section 3 of the main paper. Furthermore, we provide additional experiments
comparing the behavior of the three methods studied in Section 3 for linear
measurement operators $\mathcal{A}$ that are closer to operators that can be
encountered in applications of simultaneous low-rank and group-sparse
recovery. Finally, we shed light on the evolution of the objective function
Eq. 4 of IRLS (Algorithm 1), including in situations where the algorithm does
not manage to recover the ground truth.

###  A.1 Experimental Setup

The experiments of Section 3 were conducted using MATLAB implementations of
the three algorithms on different Linux machines using MATLAB versions R2019b
or R2022b. In total, the preparation and execution of the experiments used
approximately 1200 CPU hours. The CPU models used in the simulations are Dual
18-Core Intel Xeon Gold 6154, Dual 24-Core Intel Xeon Gold 6248R, Dual 8-Core
Intel Xeon E5-2667, 28-Core Intel Xeon E5-2690 v3, 64-Core Intel Xeon Phi KNL
7210-F. For Sparse Power Factorization (SPF) [LWB18], we used our custom
implementation of [LWB18, Algorithm 4 "rSPF_HTP"] and for Riemannian adaptive
iterative hard thresholding (RiemAdaIHT) [EKPU22], we used an implementation
provided to us by Max Pfeffer in private communications. We refer to Appendix
B for implementation details for the IRLS method Algorithm 1.

In all phase transition experiments, we define _successful recovery_ such that
the relative Frobenius error
$\frac{\left\|\mathbf{X}^{(K)}-\mathbf{X}_{\star}\right\|_{F}}{\left\|\mathbf{X}_{\star}\right\|_{F}}$
of the iterate $\mathbf{X}^{(K)}$ returned by the algorithm relative to the
simultaneously low-rank and row-sparse ground truth matrix
$\mathbf{X}_{\star}$ is smaller than the threshold $10^{-4}$. As stopping
criteria, we used the criterion that the relative change of Frobenius norm
satisfies
$\frac{\left\|\mathbf{X}^{(k)}-\mathbf{X}^{(k-1)}\right\|_{F}}{\left\|\mathbf{X}^{(k)}\right\|_{F}}<\operatorname{tol}$
for IRLS, the change in the matrix factors norms satisfy
$\left\|\mathbf{U}_{k}-\mathbf{U}_{k-1}\right\|<\operatorname{tol}$ and
$\left\|\mathbf{V}_{k}-\mathbf{V}_{k-1}\right\|<\operatorname{tol}$ for SPF,
and the norm of the Riemannian gradient in RiemAdaIHT being smaller than
$\operatorname{tol}$ for $\operatorname{tol}=10^{-10}$, or if a maximal number
of iterations is reached. This iteration threshold was chosen as
$\operatorname{max\\_iter}=250$ for IRLS and SPF and as
$\operatorname{max\\_iter}=2000$ for RiemAdaIHT, reflecting the fact that
RiemAdaIHT is a gradient-type method which might need many iterations to reach
a high-accuracy solution. The parameters were chosen so that the stopping
criteria do not prevent a method’s iterates reaching the recovery threshold if
they were to reach $\mathbf{X}_{\star}$ eventually.

In the experiments, we chose random ground truths
$\mathbf{X}_{\star}\in\mathbb{R}^{n_{1}\times n_{2}}$ of rank $r$ and row-
sparsity $s$ such that
$\mathbf{X}_{\star}=\tilde{\mathbf{X}}_{\star}/\left\|\tilde{\mathbf{X}}_{\star}\right\|_{F}$,
where
$\tilde{\mathbf{X}}_{\star}=\mathbf{U}_{\star}\mathrm{diag}(\mathbf{d}_{\star})\mathbf{V}_{\star}^{*}$,
and where $\mathbf{U}_{\star}\in\mathbb{R}^{n_{1}\times r}$ is a matrix with
$s$ non-zero rows whose location is chosen uniformly at random and whose
entries are drawn from i.i.d. standard Gaussian random variables, $\mathbf{d}$
has i.i.d. standard Gaussian entries and
$\mathbf{V}_{\star}\in\mathbb{R}^{n_{2}\times r}$ has likewise i.i.d. standard
Gaussian entries.

###  A.2 Random Rank-One Measurements

In Section 3, we considered only measurement operator
$\mathcal{A}:\mathbb{R}^{n_{1}\times n_{2}}\to\mathbb{R}^{m}$ whose matrix
representation consists of i.i.d. Gaussian entries, i.e., operators such that
there are independent matrices $\mathbf{A}_{1},\ldots\mathbf{A}_{m}$ with
i.i.d. standard Gaussian entries such that

$\mathcal{A}(\mathbf{X})_{j}=\langle\mathbf{A}_{j},\mathbf{X}\rangle_{F}$

for any $\mathbf{X}\in\mathbb{R}^{n_{1}\times n_{2}}$. While it is known that
such Gaussian measurement operators satisfy the $(r,s)$-RIP of Section 3,
which is the basis of our convergence theorem Theorem 2.4, in a regime of a
near-optimal number of measurements with high probability, practically
relevant measurement operators are often more structured; another downside of
dense Gaussian measurements is that it is computationally expensive to
implement their action on matrices.

In relevant applications of our setup, however, e.g., in sparse phase
retrieval [JOH17, CLM16, JH17] or blind deconvolution [LLJB16, SC21], the
measurement operator consists of rank-one measurements. For this reason, we
now conduct experiments in settings related to the ones depicted in Figure 1
and Figure 2 Section 3, but for _random rank-one measurements_ where the
action of $\mathcal{A}:\mathbb{R}^{n_{1}\times n_{2}}\to\mathbb{R}^{m}$ on
$\mathbf{X}$ can be written as

$\mathcal{A}(\mathbf{X})_{j}=\langle\mathbf{a}_{j}\mathbf{b}_{j}^{*},\mathbf{X}\rangle_{F}$

for each $j=1,\ldots,m$, where $\mathbf{a}_{j},\mathbf{b}_{j}$ are independent
random standard Gaussian vectors. In Figure 4, we report the phase transition
performance of RiemAdaIHT, SPF and IRLS for $(256\times 40)$-dimensional
ground truths of different row-sparsities and different ranks if we are given
such random rank-one measurements.

$10$$20$$30$$40$$50$$100$$150$$200$Number of measurements $m$

$10$$20$$30$$40$$50$$100$$150$$200$

$10$$20$$30$$40$$50$$100$$150$$200$

$20$$40$$60$$400$$600$$800$$1{,}000$Row sparsity $s$$0$$0.2$$0.4$$0.6$$0.8$$1$

$20$$40$$60$$400$$600$$800$$1{,}000$

$20$$40$$60$$400$$600$$800$$1{,}000$

Figure 4: Left column: RiemAdaIHT, center: SPF, right: IRLS. Success rates for
the recovery of low-rank and row-sparse matrices from random rank-one
measurements. First row: Rank-$1$ ground truth $\mathbf{X}_{\star}$ (cf.
Figure 1. Second row: Rank-$5$ ground truth $\mathbf{X}_{\star}$ (cf. Figure
2).

We observe in Figure 4 that compared to the setting of dense Gaussian
measurements, the phase transitions of all three algorithms deteriorate
slightly; especially for $r=1$, one can observe that the transition between no
success and high empirical success rate is extends across a larger area. IRLS
performs clearly best for both $r=1$ and $r=5$, whereas SPF has the second
best performance for $r=5$. For $r=1$, it is somewhat unclear whether
RiemAdaIHT or SPF performs better.

###  A.3 Random Rank-One Measurements

We now revisit the experiments of Section A.2 for a third measurement setup
motivated from blind deconvolution problems [ARR14, LLJB16, LLSW19, MWCC20,
SC21, EKPU22], which are prevalent in astronomy, medical imaging and
communications engineering [JC93, CE17]. In particular, in these settings, if
$\mathbf{z}\in\mathbb{R}^{m}$ is an (unknown) signal and
$\mathbf{w}\in\mathbb{R}^{m}$ is an (unknown) convolution kernel, assume we
are given the entries of their convolution
$\widetilde{y}=\mathbf{z}\ast\mathbf{w}$. If we know that
$\mathbf{z}=\mathbf{A}\mathbf{u}$ for some known matrix
$\mathbf{A}\in\mathbb{R}^{m\times n_{1}}$ and an $s$-sparse vector
$\mathbf{u}\in\mathbb{R}^{n_{1}}$ and $\mathbf{w}=\mathbf{B}\mathbf{v}$ for
some known matrix $\mathbf{B}\in\mathbb{R}^{m\times n_{2}}$ and arbitrary
vector $\mathbf{v}\in\mathbb{R}^{n_{2}}$, applying the discrete Fourier
transform (represented via the DFT matrix $\mathbf{F}\in\mathbb{C}^{m\times
m}$), we can write the coordinates of

$\mathbf{y}=\mathbf{F}\widetilde{y}=\mathrm{diag}(\mathbf{F}\mathbf{z})\mathbf{F}\mathbf{w}=\mathrm{diag}(\mathbf{F}\mathbf{A}\mathbf{u})\mathbf{F}\mathbf{B}\mathbf{v}$

as

$\mathbf{y}_{j}=\mathcal{A}(\mathbf{u}\mathbf{v}^{*})_{j}=\langle(\mathbf{F}\mathbf{A})_{j,:}^{*}\overline{\mathbf{F}\mathbf{B}}_{j,:},\mathbf{u}\mathbf{v}^{*}\rangle_{F}$

for each $j=1,\ldots,m$, which allows us to write the problem as a
simultaneously rank-$1$ and $s$-row sparse recovery problem from Fourier-type
measurements.

$10$$20$$30$$40$$50$$100$$150$$200$Number of measurements $m$

$10$$20$$30$$40$$50$$100$$150$$200$

$10$$20$$30$$40$$50$$100$$150$$200$

$20$$40$$60$$400$$600$$800$$1{,}000$Row sparsity $s$$0$$0.2$$0.4$$0.6$$0.8$$1$

$20$$40$$60$$400$$600$$800$$1{,}000$

$20$$40$$60$$400$$600$$800$$1{,}000$

Figure 5: Left column: RiemAdaIHT, center: SPF, right: IRLS. Success rates for
the recovery of low-rank and row-sparse matrices from Fourier rank-one
measurements. First row: Rank-$1$ ground truth $\mathbf{X}_{\star}$. Second
row: Rank-$5$ ground truth $\mathbf{X}_{\star}$ (cf. Figure 2).

In Figure 5, we report the results of simulations with $\mathbf{A}$ and
$\mathbf{B}$ chosen generically as standard real Gaussians for these Fourier-
based rank-$1$ measurements (including for rank-$5$ ground truths, which goes
beyond a blind deconvolution setting). We observe that the transition from no
recovery to exact recovery for an increasing number of measurement (with fixed
dimension parameters $s$, $n_{1}$ and $n_{2}$) happens earlier than for the
random Gaussian rank-one measurements of Section A.2, but slightly later than
for dense Gaussian measurements. Again, IRLS exhibits the best empirical data-
efficiency with sharpest phase transition curves.

As a summary, we observe that IRLS is able to recovery simultaneously low-rank
and row-sparse matrices empirically from fewer measurements than state-of-the-
art methods for a variety of linear measurement operators, including in cases
where the RIP assumption of Definition 2.3 is not satisfied and in cases that
are relevant for applications.

###  A.4 Evolution of Objective Values

While Theorem 2.4 guarantees local convergence if the measurement operator
$\mathcal{A}$ is generic enough and contains enough measurements (RIP-
assumption), it is instructive to study the behavior of Algorithm 1 in
situations where there are _not_ enough measurements available to identify a
specific low-rank and row-sparse ground truth $\mathbf{X}_{\star}$ which
respect to which the measurements have been taken.

In this setting, Theorem 2.5 guarantees that the behavior of the IRLS methods
is still benign as the sequence of $\varepsilon$\- and $\delta$-smoothed log-
objectives
$\left(\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k)})\right)_{k\geq
1}$ from Eq. 4 is non-increasing. In Figure 6, we illustrate the evolution of
the relative Frobenius error of an iterate to the ground truth
$\mathbf{X}_{\star}$, the $(\varepsilon_{k},\delta_{k})$-smoothed logarithmic
surrogates $\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k)})$ as
well as of the rank and row-sparsity parts
$\mathcal{F}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)})$ and
$\mathcal{F}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k)})$ of the objective,
respectively, in two typical situations.

In particular, we can see the evolution of these four quantities in the
setting of data of dimensionality $n_{1}=128$, $n_{2}\in\\{20,40\\}$, $s=20$
and $r=5$ created as in the other experiments, where a number of $m=875$ and
$m=175$ (corresponding to an oversampling factor of $3.0$ and $1.0$,
respectively) dense Gaussian measurements are provided to Algorithm 1.

In the left plot of Figure 6, which corresponds to setting of abundant
measurements, we observe that the four quantities all track each other
relatively well on a semilogarithmic scale (note that we plot the square roots
of the objective values to match the order of the (unsquared) relative
Frobenius error), converging to values between $10^{-13}$ and $10^{-11}$ (at
which point the stopping criterion of the method applies) within $12$
iterations.

In the second plot of Figure 6, the number of measurements exactly matches the
number of degrees of freedom of the ground truth, in which case the
$\mathbf{X}^{(k)}$ does _not_ converge to $\mathbf{X}_{\star}$. However, it
can be seen that Algorithm 1 still finds very meaningful solutions: It can be
seen that within $86$ iterations,
$\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k)})$ converges to
$\approx 10^{-12}$ (since
$\sqrt{\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k)})}\approx
10^{-6}$) in a manner that is partially “staircase-like“: After $20$ initial
iterations where $\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k)})$
decreases significantly at each iteration, its decrease is dominated by
relatively sudden, alternating drops of the (blue) sparsity objective
$\mathcal{F}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k)})$ and the (red) rank
objective $\mathcal{F}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)})$, which
typically do not occur simultaneously.

This illustrates the _self-balancing_ property of the two objective terms in
the IRLS objective
$\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k)})$: while the final
iterate at iteration $k=86$ is not of the target row-sparsity $s=20$ and
$r=5$, it is still $20$-row sparse and has essentially rank $6$. This means
that Algorithm 1 has found an alternative parsimonious solution to the
simultaneous low-rank and row-sparse recovery problem that is just slightly
less parsimonious.

Arguably, this robust performance in the low-data regime of IRLS is rather
unique, and to the best of our knowledge, not shared by methods such as SPF or
RiemAdaIHT, which typically breakdown in such a regime.

$0$$2$$4$$6$$8$$10$$12$$10^{-14}$$10^{-10}$$10^{-6}$$10^{-2}$$10^{2}$iteration
$k$Rank Objective
$\sqrt{\mathcal{F}_{lr,\varepsilon_{k}}(\mathbf{X}^{(k)})}$Sparsity Objective
$\sqrt{\mathcal{F}_{sp,\delta_{k}}(\mathbf{X}^{(k)})}$IRLS Objective
$\sqrt{\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k)})}$Rel. Frob.
error $\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|_{F}/\|\mathbf{X}_{\star}\|_{F}$

$0$$20$$40$$60$$80$$10^{-7}$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$iteration
$k$

Figure 6: Objective/ error quantities of iterates $\mathbf{X}^{(k)}$ for
iterations $k$. Left: Typical result for $n_{1}=128$, $n=40$, $m=875$. Right:
Typical result for $n_{1}=128$, $n=20$, $m=175$.

##  Appendix B Discussion of Computational Complexity

It is well-known that the solution of the linearly constrained weighted least
squares problem Eq. 11 can be written as

(18)
$\mathbf{X}^{(k)}=W_{k-1}^{-1}\mathcal{A}^{*}\left(\mathcal{A}W_{k-1}^{-1}\mathcal{A}^{*}\right)^{-1}\mathbf{y}$

where $W_{k-1}:=W_{\mathbf{X}^{(k-1)},\varepsilon_{k-1},\delta_{k-1}}$ is the
weight operator Eq. 8 of iteration $k-1$ [DDFG10, KMV21]. In [KMV21, Theorem
3.1 and Supplementary Material], it was shown that in the case of low-rank
matrix completion without the presence of a row-sparsity inducing term, this
weighted least squares problem can be solved by solving an equivalent, well-
conditioned linear system via an iterative solver that uses the application of
a system matrix whose matrix-vector products have time complexity of
$O(mr+r^{2}\max(n_{1},n_{2}))$.

In the case of Algorithm 1, the situations is slightly more involved as we
cannot provide an explicit formula for the inverse of the weight operator
$W_{k-1}$ as it amounts to the sum of the weight operators
$W_{\mathbf{X}^{(k-1)},\varepsilon_{k-1}}^{\textit{lr}}$ and
$\mathbf{W}_{\mathbf{X}^{(k-1)},\delta_{k-1}}^{\textit{sp}}$ that are
diagonalized by different, mutually incompatible bases. However, computing
this inverse is facilitated by the _Sherman-Morrison-Woodbury_ formula [Woo50]

$(\mathbf{E}\mathbf{C}\mathbf{F}^{*}+\mathbf{B})^{-1}=\mathbf{B}^{-1}-\mathbf{B}^{-1}\mathbf{E}(\mathbf{C}^{-1}+\mathbf{F}^{*}\mathbf{B}^{-1}\mathbf{E})^{-1}\mathbf{F}^{*}\mathbf{B}^{-1}$

for suitable matrices of compatible size $\mathbf{E},\mathbf{F}$ and
invertible $\mathbf{C},\mathbf{B}$ and the fact that both
$W_{\mathbf{X}^{(k-1)},\varepsilon_{k-1}}^{\textit{lr}}$ and
$\mathbf{W}_{\mathbf{X}^{(k-1)},\delta_{k-1}}^{\textit{sp}}$ exhibit a “low-
rank plus (scaled) identity“ or a “sparse diagonal plus (scaled) identity“
structure. After a simple application of the SMW formula, Eq. 18 can be
rewritten such that the computational bottleneck becomes the assembly and
inversion of a $O(r_{k}\max(n_{1},n_{2}))$ linear system. We refer to our
MATLAB implementation for further details.

With this implementation, which is not optimized for large-scale problems, the
computational cost becomes comparable to the implementations of SPF or
RiemAdaIHT we used. We leave further improvements and adaptations to large-
scale settings to future work.

##  Appendix C Proofs

The following two sections contain the proofs of our main results. Let us
begin with some helpful observations.

First note that the low-rank promoting part
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}:\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}^{n_{1}\times n_{2}}$ of our weight operator can be re-
written as

(19)
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{Z})=\begin{bmatrix}\mathbf{U}&\mathbf{U}_{\perp}\end{bmatrix}\left(\mathbf{H}(\boldsymbol{\sigma}^{(k)},\varepsilon_{k})\circ\left(\begin{bmatrix}\mathbf{U}^{*}\\\
\mathbf{U}_{\perp}^{*}\end{bmatrix}\mathbf{Z}\begin{bmatrix}\mathbf{V}&\mathbf{V}_{\perp}\end{bmatrix}\right)\right)\begin{bmatrix}\mathbf{V}^{*}\\\
\mathbf{V}_{\perp}^{*}\end{bmatrix},$

where

$\displaystyle\mathbf{H}(\boldsymbol{\sigma}^{(k)},\varepsilon_{k})$
$\displaystyle:=\left[\min\left(\varepsilon_{k}/\sigma_{i}^{(k)},1\right)\min\left(\varepsilon_{k}/\sigma_{j}^{(k)},1\right)\right]_{i,j=1}^{n_{1},n_{2}}$
$\displaystyle=\left[\begin{array}[]{c|c}\left(\frac{\varepsilon_{k}^{2}}{\sigma_{i}^{(k)}\sigma_{j}^{(k)}}\right)_{i,j=1}^{r_{k}}&\left(\frac{\varepsilon_{k}}{\sigma_{i}^{(k)}}\right)_{i,j=1}^{r_{k},d_{2}}\\\
\hline\cr\left(\frac{\varepsilon_{k}}{\sigma_{j}^{(k)}}\right)_{i,j=1}^{d_{1},r_{k}}&\mathbf{1}\end{array}\right]\in\mathbb{R}^{n_{1}\times
n_{2}}.$

Consequently, all weight operators in Definition 2.1 are self-adjoint and
positive. Whereas for $\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}$
this is obvious, for $W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}$ it
follows from the matrix representation

$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}=\Big{(}\begin{bmatrix}\mathbf{U}&\mathbf{U}_{\perp}\end{bmatrix}\otimes\begin{bmatrix}\mathbf{V}&\mathbf{V}_{\perp}\end{bmatrix}\Big{)}\mathbf{D}_{\mathbf{H}(\boldsymbol{\sigma}^{(k)},\varepsilon_{k})}\Big{(}\begin{bmatrix}\mathbf{U}&\mathbf{U}_{\perp}\end{bmatrix}^{*}\otimes\begin{bmatrix}\mathbf{V}&\mathbf{V}_{\perp}\end{bmatrix}^{*}\Big{)}$

where
$\mathbf{D}_{\mathbf{H}(\boldsymbol{\sigma}^{(k)},\varepsilon_{k})}\in\mathbb{R}^{n_{1}n_{2}\times
n_{1}n_{2}}$ is a diagonal matrix with the entries of
$\mathbf{H}(\boldsymbol{\sigma}^{(k)},\varepsilon_{k})$, which are all
positive, on its diagonal.

###  C.1 Proof of Theorem 2.4

Before approaching the proof of Theorem 2.4, let us collect various important
observations. In order to keep the presentation concise, we defer part of the
proofs to Appendix D.

For a rank-$r$ matrix
$\mathbf{Z}=\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{*}$, we define the
tangent space of the manifold of rank-$r$ matrices at $\mathbf{Z}$ as

(20) $\displaystyle
T_{\mathbf{U},\mathbf{V}}:=\\{\mathbf{U}\mathbf{Z}_{1}^{*}+\mathbf{Z}_{2}\mathbf{V}^{*}\colon\mathbf{Z}_{1}\in\mathbb{R}^{n_{2}\times
r},\mathbf{Z}_{2}\in\mathbb{R}^{n_{1}\times r}\\}.$

In a similar manner, we can define for
$\mathbf{Z}=\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{*}\in\mathcal{M}_{r,s}$
and $S=\mathrm{supp}(\mathbf{Z})=\mathrm{supp}(\mathbf{U})\subset[n_{1}]$ the
tangent space of $\mathcal{M}_{r}$ restricted to $S$ as

(21) $\displaystyle
T_{\mathbf{U},\mathbf{V},S}:=\\{\mathbf{U}\mathbf{Z}_{1}^{*}+\mathbf{Z}_{2}\mathbf{V}^{*}\colon\mathbf{Z}_{1}\in\mathbb{R}^{n_{2}\times
r},\mathbf{Z}_{2}\in\mathbb{R}^{n_{1}\times r}\text{ with
}\mathrm{supp}(\mathbf{Z}_{2})=S\\}.$

As the following lemma shows, orthogonal projections onto the sets
$\mathcal{M}_{r}^{n_{1},n_{2}}$, $\mathcal{N}_{s}^{n_{1},n_{2}}$,
$T_{\mathbf{U},\mathbf{V}}$, and $T_{\mathbf{U},\mathbf{V},S}$ can be
efficiently computed.

######  Lemma C.1.

We denote the projection operators onto $\mathcal{M}_{r}^{n_{1},n_{2}}$ and
$\mathcal{N}_{s}^{n_{1},n_{2}}$ by $\mathsf{T}_{r}$ and $\mathsf{H}_{s}$.
$\mathsf{T}_{r}$ truncates a matrix to the $r$ dominant singular values;
$\mathsf{H}_{s}$ sets all but the $s$ in $\ell_{2}$-norm largest rows to zero.
In case of ambiguities (multiple singular values/rows of same magnitude), by
convention we choose the $r$ (respectively $s$) with smallest index.

For $\mathbf{U}$ and $\mathbf{V}$ fixed, the orthogonal projection onto
$T_{\mathbf{U},\mathbf{V}}$ is given by

$\displaystyle\mathbb{P}_{\mathbf{U},\mathbf{V}}:=\mathbb{P}_{T_{\mathbf{U},\mathbf{V}}}\mathbf{Z}=\mathbf{U}\mathbf{U}^{*}\mathbf{Z}+\mathbf{Z}\mathbf{V}\mathbf{V}^{*}-\mathbf{U}\mathbf{U}^{*}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}.$

For $S\subset[n_{1}]$ and $\mathbf{U},\mathbf{V}$ fixed with
$\mathrm{supp}(\mathbf{U})=S$, the orthogonal projection onto
$T_{\mathbf{U},\mathbf{V},S}$ is given by

$\displaystyle\mathbb{P}_{\mathbf{U},\mathbf{V},S}:=\mathbb{P}_{T_{\mathbf{U},\mathbf{V},S}}\mathbf{Z}$
$\displaystyle=\mathbb{P}_{S}(\mathbf{U}\mathbf{U}^{*}\mathbf{Z}+\mathbf{Z}\mathbf{V}\mathbf{V}^{*}-\mathbf{U}\mathbf{U}^{*}\mathbf{Z}\mathbf{V}\mathbf{V}^{*})$
$\displaystyle=\mathbf{U}\mathbf{U}^{*}\mathbf{Z}+\mathbb{P}_{S}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}-\mathbf{U}\mathbf{U}^{*}\mathbf{Z}\mathbf{V}\mathbf{V}^{*},$

where $\mathbb{P}_{S}$ projects to the row support $S$, i.e., it sets all rows
to zero which are not indexed by $S$.

The proof of Lemma C.1 is provided in Section D.2. In contrast to the above
named projections, the projection onto $\mathcal{M}_{r,s}^{n_{1},n_{2}}$ is
not tractable. However, [EKPU22, Lemma 2.4] shows that locally
$\mathbb{P}_{\mathcal{M}_{r,s}}$ can be replaced by the concatenation of
$\mathsf{T}_{r}$ and $\mathsf{H}_{s}$, i.e., for
$\mathbf{Z}_{\star}\in\mathcal{M}_{r,s}$ and
$\mathbf{Z}\approx\mathbf{Z}_{\star}$, one has that

$\displaystyle\mathbb{P}_{\mathcal{M}_{r,s}}(\mathbf{Z})=\mathsf{T}_{r}(\mathsf{H}_{s}(\mathbf{Z})).$

For a matrix $\mathbf{X}\in\mathbb{R}^{n_{1}\times n_{2}}$ and $i\in[n_{1}]$,
we set $\rho_{i}(\mathbf{X})=\|(\mathbf{X})_{i^{\prime},:}\|_{2}$ where
$i^{\prime}$ is a row index corresponding to the $i$-th largest row of
$\mathbf{X}$ in $\ell_{2}$-norm. More precisely, if $\overline{\mathbf{X}}$ is
a decreasing rearrangement of $\mathbf{X}$ with rows ordered by magnitude in
$\ell_{2}$-norm, then
$\rho_{i}(\mathbf{X})=\|(\overline{\mathbf{X}})_{i,:}\|_{2}$. As the following
lemma shows, the quantity $\rho_{s}(\mathbf{X})$ determines a local
neighborhood of $\mathbf{X}$ on which $\mathsf{H}_{s}$ preserves the row-
support.

######  Lemma C.2.

Let $\mathbf{X}\in\mathbb{R}^{n_{1}\times n_{2}}$ be a matrix with row-support
$S\subset[n_{1}]$ and $|S|=s$. Then, for any
$\mathbf{Z}\in\mathbb{R}^{n_{1}\times n_{2}}$ with
$\|\mathbf{X}-\mathbf{Z}\|_{\infty,2}:=\max_{i\in[n_{1}]}\|\mathbf{X}_{i,:}-\mathbf{Z}_{i,:}\|_{2}\leq\frac{1}{2}\rho_{s}(\mathbf{X})$
the matrix $\mathsf{H}_{s}(\mathbf{Z})$ has row-support $S$.

######  Proof C.3.

Note that

$\max_{i\in[n_{1}]}\|(\mathbf{Z})_{i,:}-(\mathbf{X})_{i,:}\|_{2}\leq\frac{1}{2}\rho_{s}(\mathbf{X})$

implies that any non-zero row of $\mathbf{X}$ corresponds to a non-zero row of
$\mathsf{H}_{s}(\mathbf{Z})$ and hence yields the claim.

A first important observation is that if $\mathcal{A}$ has the $(r,s)$-RIP,
then the norm of kernel elements of $\mathcal{A}$ is bounded in the following
way.

######  Lemma C.4.

If $\mathcal{A}$ has the $(r,s)$-RIP with $\delta\in(0,1)$ and
$\mathbf{U}\in\mathbb{R}^{n_{1}\times
r}$,$\mathbf{V}\in\mathbb{R}^{n_{2}\times r}$ with
$\mathrm{supp}(\mathbf{U})=S$, $|S|\leq s$, then

$\displaystyle\left\|\boldsymbol{\Xi}\right\|_{F}\leq\sqrt{1+\frac{\left\|\mathcal{A}\right\|_{2\to
2}^{2}}{(1-\delta)}}\left\|\mathbb{P}_{T_{\mathbf{U},\mathbf{V},S}}^{\perp}(\boldsymbol{\Xi})\right\|_{F},$

for all $\boldsymbol{\Xi}\in\ker(\mathcal{A})$.

The proof of Lemma C.4 is presented in Section D.3.

######  Remark C.5.

If $\mathcal{A}$ is a Gaussian operator with standard deviation
$\sqrt{\frac{1}{m}}$, one has with high probability that
$\left\|\mathcal{A}\right\|_{2\to 2}^{2}\approx\frac{n_{1}n_{2}}{m}$.

We use of the following lemma to characterize the solution of the weighted
least squares problem Eq. 11. Its proof is analogous to [KMV21, Lemma B.7] and
[DDFG10, Lemma 5.2].

######  Lemma C.6.

Let $\mathcal{A}:\mathbb{R}^{n_{1}\times n_{2}}\to\mathbb{R}^{m}$ and
$\mathbf{y}\in\mathbb{R}^{m}$. Let
$W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}:\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}^{n_{1}\times n_{2}}$ be the weight operator Eq. 8 defined
based on the information of $\mathbf{X}^{(k)}\in\mathbb{R}^{n_{1}\times
n_{2}}$. Then the solution of the weighted least squares step Eq. 11 of
Algorithm 1

(22)
$\mathbf{X}^{(k+1)}=\operatorname*{arg\,min}\limits_{\mathcal{A}(\mathbf{X})=\mathbf{y}}\langle\mathbf{X},W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}(\mathbf{X})\rangle,$

is unique and solves Eq. 22 if and only if

(23) $\mathcal{A}(\mathbf{X}^{(k+1)})=\mathbf{y}\qquad\text{ and
}\qquad\langle
W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k+1)}),\boldsymbol{\Xi}\rangle=0\;\text{
for all }\;\boldsymbol{\Xi}\in\ker\mathcal{A}.$

For any iterate $\mathbf{X}^{(k)}$ of Algorithm 1, we furthermore abbreviate
the tangent space Eq. 21 of the fixed rank-$r$ manifold $\mathcal{M}_{r}$
restricted to $S$ at $\mathsf{H}_{s}(\mathbf{X}^{(k)})$ by

(24) $\displaystyle
T_{k}=T_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}},S},$

where $\widetilde{\mathbf{U}}\in\mathbb{R}^{n_{1}\times r}$ and
$\widetilde{\mathbf{V}}\in\mathbb{R}^{n_{2}\times r}$ are matrices with
leading444As $\widetilde{\mathbf{U}}$ and $\widetilde{\mathbf{V}}$ might not
be unique, any set of $r$ leading singular vectors can be chosen in this
definition. $r$ singular vectors of $\mathsf{H}_{s}(\mathbf{X}^{(k)})$ as
columns, and $S\in[n_{1}]$ is the support set of the $s$ rows of
$\mathbf{X}^{(k)}$ with largest $\ell_{2}$-norm.

The following lemma is the first crucial tool for showing local quadratic
convergence of Algorithm 1.

######  Lemma C.7.

Let $\mathbf{X}_{\star}\in\mathcal{M}_{r,s}$ and let $\mathbf{X}^{(k)}$ be the
$k$-th iterate of Algorithm 1 with rank and sparsity parameters
$\widetilde{r}=r$ and $\widetilde{s}=s$, let $\delta_{k},\varepsilon_{k}$ be
such that $s_{k}$ and $r_{k}$ from Definition 2.1 satisfy $s_{k}\geq s$ and
$r_{k}\geq r$. Assume that there exists a constant $c>1$ such that

(25) $\displaystyle\left\|\boldsymbol{\Xi}\right\|_{F}\leq
c\left\|\mathbb{P}_{T_{k}^{\perp}}(\boldsymbol{\Xi})\right\|_{F}\qquad\text{for
all }\quad\boldsymbol{\Xi}\in\ker(\mathcal{A}),$

where $T_{k}=T_{\widetilde{\mathbf{U}},\widetilde{\mathbf{V}},S}$ is as
defined in Eq. 24 for matrices
$\widetilde{\mathbf{U}}\in\mathbb{R}^{n_{1}\times r}$ and
$\widetilde{\mathbf{V}}\in\mathbb{R}^{n_{2}\times r}$ of leading $r$ left and
right singular vectors of $\mathsf{H}_{s}(\mathbf{X}^{(k)})$ and
$S\subset[n_{1}]$ is the support set of $\mathsf{H}_{s}(\mathbf{X}^{(k)})$.
Assume furthermore that

(26)
$\displaystyle\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|\leq\min\left\\{\frac{1}{2}\rho_{s}(\mathbf{X}_{\star}),\min\Big{\\{}\frac{1}{48},\frac{1}{19c}\Big{\\}}\sigma_{r}(\mathbf{X}_{\star})\right\\}.$

Then,

$\displaystyle\left\|\mathbf{X}^{(k+1)}-\mathbf{X}_{\star}\right\|$
$\displaystyle\leq
4c^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}$
$\displaystyle\qquad\cdot\left(\left\|W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{X}_{\star})\right\|_{*}+\left\|\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\cdot\mathbf{X}_{\star}\right\|_{1,2}\right),$

where
$\left\|\mathbf{M}\right\|_{1,2}=\sum_{i}\left\|\mathbf{M}_{i,:}\right\|_{2}$
denotes the row-sum norm of a matrix $\mathbf{M}$, and
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}$ and
$\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}$ are the weight
operators Eq. 9 and Eq. 10 from Definition 2.1.

The proof of Lemma C.7 is presented in Section D.4.

######  Remark C.8.

By revisiting the proof of Lemma C.7 (omit the bound in Eq. 48 and keep the
term $\langle\boldsymbol{\Xi},\bar{W}\boldsymbol{\Xi}\rangle$ until the end),
one can show under the same assumptions as in Lemma C.7 that

$\|\boldsymbol{\Xi}\|_{F}^{2}\leq
4c^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}\Big{\langle}\boldsymbol{\Xi},\big{(}\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}\mathbb{P}_{\mathbf{U},\mathbf{V}}^{\perp}+\mathbb{P}_{S^{c}}\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbb{P}_{S^{c}}\big{)}\boldsymbol{\Xi}\Big{\rangle},$

where $\mathbf{U}$ and $\mathbf{V}$ are containing the left and right singular
vectors of $\mathbf{X}^{(k)}$, see Definition 2.1.

The contribution of the norms of the weighted $\mathbf{X}_{\star}$ terms in
Lemma C.7 can be controlled by Lemmas C.9 and C.10 below.

######  Lemma C.9.

Let
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}\colon\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}^{n_{1}\times n_{2}}$ be the rank-based weight operator Eq.
9 that uses the spectral information of $\mathbf{X}^{(k)}$ and let
$\mathbf{X}_{\star}\in\mathbb{R}^{n_{1}\times n_{2}}$ be a rank-$r$ matrix.
Assume that there exists $0<\zeta<\frac{1}{2}$ such that

(27)
$\max\\{\varepsilon_{k},\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|\\}\leq\zeta\sigma_{r}(\mathbf{X}_{\star}).$

Then for each $1\leq q\leq\infty$,

$\left\|W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{X}_{\star})\right\|_{S_{q}}\leq\frac{r^{1/q}}{(1-\zeta)\sigma_{r}(\mathbf{X}_{\star})}\left(\frac{1}{1-\zeta}\varepsilon_{k}^{2}+\varepsilon_{k}K_{q}\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|+2\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|^{2}\right)$

and

$\left\|W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{X}_{\star})\right\|_{S_{q}}\leq\frac{1}{(1-\zeta)\sigma_{r}(\mathbf{X}_{\star})}\left(\frac{r^{1/q}}{1-\zeta}\varepsilon_{k}^{2}+2\left\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\right\|_{S_{q}}\left(\varepsilon_{k}+\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|\right)\right)$

where $K_{q}$ is such that $K_{q}=2^{1/q}$ for $1\leq q\leq 2$ and $4\leq q$,
$K_{q}=\sqrt{2}$ for $2<q\leq 4$ and $K_{q}=1$ for $q=\infty$.

######  Lemma C.10.

Let
$\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\in\mathbb{R}^{n_{1}\times
n_{1}}$ be the row-sparsity-based weight operator Eq. 10 that uses the current
iterate $\mathbf{X}^{(k)}$ with
$\delta_{k}=\min\left(\delta_{k-1},\rho_{s+1}(\mathbf{X}^{(k)})\right)$ and
let $\mathbf{X}_{\star}\in\mathbb{R}^{n_{1}\times n_{2}}$ be an $s$-row-sparse
matrix. Assume that there exists $0<\zeta<\frac{1}{2}$ such that

(28)
$\displaystyle\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|_{\infty,2}=\max_{i\in[n_{1}]}\|(\mathbf{X}^{(k)})_{i,:}-(\mathbf{X}_{\star})_{i,:}\|_{2}\leq\zeta\rho_{s}(\mathbf{X}_{\star}),$

where $\rho_{s}(\mathbf{M})$ denotes the $\ell_{2}$-norm of the in
$\ell_{2}$-norm $s$-largest row of $\mathbf{M}$. Then

$\displaystyle\|\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\leavevmode\nobreak\
\cdot\leavevmode\nobreak\
\mathbf{X}_{\star}\|_{1,2}\leq\frac{s\delta_{k}^{2}}{(1-\zeta)^{2}\rho_{s}(\mathbf{X}_{\star})}$

Lemma C.9 is a refined version of [KMV21, Lemma B.9] the proof of which we
omit here.555This result is a technical result of an unpublished paper. In
this paper, we only use that result as a tool. If the reviewers think that
adding the proof is relevant here, we are happy to provide it. The proof of
Lemma C.10 is provided in Section D.5. Finally, the following lemma will allow
us to control the decay of the IRLS parameters $\delta_{k}$ and
$\varepsilon_{k}$.

######  Lemma C.11 ([KMV21, Lemma B.5]).

Let $\mathbf{X}_{\star}\in\mathcal{M}_{r,s}$, assume that $\mathcal{A}$ has
the $(r,s)$-RIP with $\delta\in(0,1)$, and let us abbreviate
$n=\min\\{n_{1},n_{2}\\}$.

Assume that the $k$-th iterate $\mathbf{X}^{(k)}$ of Algorithm 1 with
$\widetilde{r}=r$ and $\widetilde{s}=r$ updates the smoothing parameters in
Eq. 12 such that one of the statements
$\varepsilon_{k}=\sigma_{r+1}(\mathbf{X}^{(k)})$ or
$\delta_{k}=\rho_{s+1}(\mathbf{X}^{(k)})$ is true, and that $r_{k}\geq r$ and
$s_{k}\geq s$. Furthermore, let

$\displaystyle\varepsilon_{k}\leq\frac{1}{48}\sigma_{r}(\mathbf{X}_{\star})$

with $c_{\left\|\mathcal{A}\right\|_{2\to
2}}=\sqrt{1+\tfrac{\left\|\mathcal{A}\right\|_{2\to 2}^{2}}{(1-\delta)}}$, let
$\boldsymbol{\Xi}^{(k)}:=\mathbf{X}^{(k)}-\mathbf{X}_{\star}$ satisfy

(29)
$\displaystyle\|\boldsymbol{\Xi}^{(k)}\|\leq\min\Big{\\{}\frac{1}{2}\rho_{s}(\mathbf{X}_{\star}),\min\Big{\\{}\frac{1}{48},\frac{1}{21c_{\left\|\mathcal{A}\right\|_{2\to
2}}}\Big{\\}}\sigma_{r}(\mathbf{X}_{\star})\Big{\\}}.$

Then

$\displaystyle\|\boldsymbol{\Xi}^{(k)}\|_{F}\leq
2\sqrt{2}\sqrt{n}c_{\left\|\mathcal{A}\right\|_{2\to
2}}\sqrt{4\varepsilon_{k}^{2}+\delta_{k}^{2}}.$

The proof of Lemma C.11 is provided in Section D.6. We finally have all the
tools to prove Theorem 2.4. Note that Eq. 15 implies

(30) $\displaystyle\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|$
$\displaystyle\leq\min\left\\{\frac{1}{48c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{2}}\min\left\\{\frac{\sigma_{r}(\mathbf{X}_{\star})}{r},\frac{\rho_{s}(\mathbf{X}_{\star})}{s}\right\\},\frac{1}{4\mu\sqrt{5n}c_{\left\|\mathcal{A}\right\|_{2\to
2}}}\right\\},$

which we will use in the proof below.

######  Proof C.12 (of Theorem 2.4).

First note, that by assumption $\tilde{r}=r$ and $\tilde{s}=s$. Furthermore,
since $\frac{1}{48c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{2}}\leq\frac{1}{2}$, the closeness assumption Eq. 30 implies that
$\mathsf{H}_{s}(\mathbf{X}^{(k)})$ and $\mathbf{X}_{\star}$ share the same
support due to Lemma C.2.

Let $\mathbf{X}^{(k)}$ be the $k$-th iterate of Algorithm 1. Since the
operator $\mathcal{A}\colon\mathbb{R}^{n_{1}\times n_{2}}\to\mathbb{R}^{m}$
has the $(r,s)$-RIP with $\delta\in(0,1)$, Lemma C.4 yields for all
$\mathbf{U}\in\mathbb{R}^{n_{1}\times
r}$,$\mathbf{V}\in\mathbb{R}^{n_{2}\times r}$ with
$\mathrm{supp}(\mathbf{U})=S$, $|S|\leq s$, that

$\displaystyle\left\|\boldsymbol{\Xi}\right\|_{F}\leq
c_{\left\|\mathcal{A}\right\|_{2\to
2}}\left\|\mathbb{P}_{T_{\mathbf{U},\mathbf{V},S}}^{\perp}(\boldsymbol{\Xi})\right\|_{F},$

for any $\boldsymbol{\Xi}\in\ker(\mathcal{A})$. Furthermore, due to our
assumption that $\tilde{s}=s$ and $\tilde{r}=r$, the smoothing parameter
update rules in Eq. 12, i.e.,
$\delta_{k}=\min\left(\delta_{k-1},\rho_{s+1}(\mathbf{X}^{(k)})\right)$ and
$\varepsilon_{k}=\min\left(\varepsilon_{k-1},\sigma_{r+1}(\mathbf{X}^{(k)})\right)$,
imply that $r_{k}\geq r$ and $s_{k}\geq s$ for all $k$. We can thus apply
Lemma C.7 for $\boldsymbol{\Xi}^{(k)}:=\mathbf{X}^{(k)}-\mathbf{X}_{\star}$
(note at this point that Eq. 30 implies the closeness assumption Eq. 26 of
Lemma C.7) and obtain

(31)
$\begin{split}\left\|\boldsymbol{\Xi}^{(k+1)}\right\|&=\left\|\mathbf{X}^{(k+1)}-\mathbf{X}_{\star}\right\|\\\
&\leq 4c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}\left(\left\|W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{X}_{\star})\right\|_{*}+\left\|\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\cdot\mathbf{X}_{\star}\right\|_{1,2}\right),\end{split}$

where
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}\colon\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}^{n_{1}\times n_{2}}$ is the low-rank promoting part Eq. 9
of the weight operator associated to $\mathbf{X}^{(k)}$ and
$\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\in\mathbb{R}^{n_{1}\times
n_{1}}$ the sparsity promoting part Eq. 10. Since by assumption

$\displaystyle\max(\varepsilon_{k},\|\boldsymbol{\Xi}^{(k)}\|)\leq\frac{1}{48}\sigma_{r}(\mathbf{X}_{\star}),$

Lemma C.9 yields

(32)
$\displaystyle\big{\|}W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{X}_{\star})\big{\|}_{*}\leq
0.995\frac{42}{40\sigma_{r}(\mathbf{X}_{\star})}\left(\varepsilon_{k}^{2}r+2\varepsilon_{k}\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|+2\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|^{2}\right).$

Similarly, by assumption

$\displaystyle\|\boldsymbol{\Xi}^{(k)}\|_{\infty,2}\leq\|\boldsymbol{\Xi}^{(k)}\|\leq\frac{1}{48s}\rho_{s}(\mathbf{X}_{\star})\leq\frac{1}{48}\rho_{s}(\mathbf{X}_{\star}),$

such that Lemma C.10 yields

(33)
$\displaystyle\|\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\leavevmode\nobreak\
\cdot\leavevmode\nobreak\ \mathbf{X}_{\star}\|_{1,2}\leq
0.995\frac{21s\delta_{k}^{2}}{20\rho_{s}(\mathbf{X}_{\star})}.$

Inserting Eq. 32 and Eq. 33 into Eq. 31 we obtain that

(34)
$\displaystyle\begin{split}\left\|\boldsymbol{\Xi}^{(k+1)}\right\|\leq&\;0.995\cdot
4.2c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{2}\min\left\\{\frac{\sigma_{r+1}(\mathbf{X}^{(k)})}{\varepsilon_{k}},\frac{\rho_{s+1}(\mathbf{X}^{(k)})}{\delta_{k}}\right\\}^{2}\\\
&\cdot\bigg{(}\frac{r}{\sigma_{r}(\mathbf{X}_{\star})}\left(\varepsilon_{k}^{2}+2\varepsilon_{k}\|\boldsymbol{\Xi}^{(k)}\|+2\|\boldsymbol{\Xi}^{(k)}\|^{2}\right)+\frac{2s}{\rho_{s}(\mathbf{X}_{\star})}\delta_{k}^{2}\bigg{)}.\end{split}$

Due to the assertion that $r_{k}\geq r$, it holds that
$\varepsilon_{k}\leq\sigma_{r+1}(\mathbf{X}^{(k)})$. Therefore, Lemma D.5
yields that

$\left(\varepsilon_{k}^{2}+2\varepsilon_{k}\|\boldsymbol{\Xi}^{(k)}\|+2\|\boldsymbol{\Xi}^{(k)}\|^{2}\right)\leq
5\|\boldsymbol{\Xi}^{(k)}\|^{2}.$

and, since $s_{k}\geq s$, also that

$\delta_{k}^{2}\leq\|\boldsymbol{\Xi}^{(k)}\|_{\infty,2}^{2}\leq\|\boldsymbol{\Xi}^{(k)}\|^{2},$

since $\delta_{k}\leq\rho_{s+1}(\mathbf{X}^{(k)})$ in this case.

Thus, using the assertion that one of the statements
$\varepsilon_{k}=\sigma_{r+1}(\mathbf{X}^{(k)})$ or
$\delta_{k}=\rho_{s+1}(\mathbf{X}^{(k)})$ is true, we obtain from Eq. 34 that

(35) $\displaystyle\begin{split}\left\|\boldsymbol{\Xi}^{(k+1)}\right\|\leq
0.995\cdot 4.2c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{2}\cdot\left(\frac{5r}{\sigma_{r}(\mathbf{X}_{\star})}+\frac{2s}{\rho_{s}(\mathbf{X}_{\star})}\right)\left\|\boldsymbol{\Xi}^{(k)}\right\|^{2}.\end{split}$

For
$\left\|\boldsymbol{\Xi}^{(k)}\right\|<\frac{1}{48}c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{-2}\min\\{\tfrac{\sigma_{r}(\mathbf{X}_{\star})}{r},\tfrac{\rho_{s}(\mathbf{X}_{\star})}{s}\\}$
(as implied by Eq. 30), this yields

(36) $\displaystyle\|\boldsymbol{\Xi}^{(k+1)}\|<0.9\|\boldsymbol{\Xi}^{(k)}\|$

and the quadratic error decay

$\displaystyle\|\boldsymbol{\Xi}^{(k+1)}\|\leq\mu\|\boldsymbol{\Xi}^{(k)}\|^{2}$

if we define $\mu=4.179c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{2}\big{(}\frac{5r}{\sigma_{r}(\mathbf{X}_{\star})}+\frac{2s}{\rho_{s}(\mathbf{X}_{\star})}\big{)}$.

To show the remaining statement, we need to argue that the assertions of
Theorem 2.4 are satisfied not only for $k$, but for any $k+\ell$ with
$\ell\geq 1$. For this, it is sufficient to show that

  1. 1.

$r_{k+1}\geq r$,

  2. 2.

$s_{k+1}\geq s$,

  3. 3.

$\varepsilon_{k+1}\leq\frac{1}{48}\sigma_{r}(\mathbf{X}_{\star})$,

  4. 4.

Eq. 30 holds for $\mathbf{X}^{(k+1)}$, and that

  5. 5.

one of the statements $\varepsilon_{k+1}=\sigma_{r+1}(\mathbf{X}^{(k+1)})$ or
$\delta_{k+1}=\rho_{s+1}(\mathbf{X}^{(k+1)})$ is true,

as in this case,
$\mathbf{X}^{(k+\ell)}\overset{\ell\to\infty}{\to}\mathbf{X}_{\star}$ follows
by induction due to successive application of Eq. 36.

For 1. and 2., we see that this follows from the smoothing parameter update
rules Eq. 12 which imply that
$\varepsilon_{k+1}\leq\sigma_{r+1}(\mathbf{X}^{(k+1)})$ and
$\delta_{k+1}\leq\rho_{s+1}(\mathbf{X}^{(k+1)})$.

3\. follows from Eq. 14 and the fact that due to Eq. 12,
$(\varepsilon_{k})_{k\geq 1}$ is non-increasing. 4. is satisfied due to Eq. 36
and Eq. 30.

To show 5., we note that due to Eq. 30, the assertion Eq. 29 is satisfied, and
therefore it follows from Eq. 35 and Lemma C.11 that

$\left\|\boldsymbol{\Xi}^{(k+1)}\right\|\leq
4.179c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{2}\left(\frac{5r}{\sigma_{r}(\mathbf{X}_{\star})}+\frac{2s}{\rho_{s}(\mathbf{X}_{\star})}\right)\left\|\boldsymbol{\Xi}^{(k)}\right\|\cdot
2\sqrt{2}\sqrt{n}c_{\left\|\mathcal{A}\right\|_{2\to
2}}\sqrt{4\varepsilon_{k}^{2}+\delta_{k}^{2}}.$

We now distinguish the case (i) $\delta_{k}<\varepsilon_{k}$ and the case (ii)
$\delta_{k}\geq\varepsilon_{k}$.

In case (i), it holds that

$\begin{split}\sigma_{r+1}(\mathbf{X}^{(k+1)})\leq\left\|\boldsymbol{\Xi}^{(k+1)}\right\|&\leq
4.179c_{\left\|\mathcal{A}\right\|_{2\to
2}}^{2}\left(\frac{5r}{\sigma_{r}(\mathbf{X}_{\star})}+\frac{2s}{\rho_{s}(\mathbf{X}_{\star})}\right)2\sqrt{10n}c_{\left\|\mathcal{A}\right\|_{2\to
2}}\left\|\boldsymbol{\Xi}^{(k)}\right\|\varepsilon_{k}\\\ &=\mu
2\sqrt{10n}c_{\left\|\mathcal{A}\right\|_{2\to
2}}\left\|\boldsymbol{\Xi}^{(k)}\right\|\varepsilon_{k}\\\
&<\varepsilon_{k},\end{split}$

where the last inequality holds since by (30) the $k$-th iterate
$\mathbf{X}^{(k)}$ additionally satisfies

(37)
$\|\mathbf{X}^{(k)}-\mathbf{X}_{\star}\|<\frac{1}{2\mu\sqrt{10}c_{\left\|\mathcal{A}\right\|_{2\to
2}}}.$

In this case, due to the smoothing parameter update rule Eq. 12, we have that
$\varepsilon_{k+1}=\sigma_{r+1}(\mathbf{X}^{(k+1)})$.

In case (ii), we have likewise that

$\rho_{s+1}(\mathbf{X}^{(k+1)})\leq\left\|\boldsymbol{\Xi}^{(k+1)}\right\|\leq\mu
2\sqrt{10n}c_{\left\|\mathcal{A}\right\|_{2\to
2}}\left\|\boldsymbol{\Xi}^{(k)}\right\|\delta_{k}<\delta_{k},$

due to Eq. 35, Lemma C.11, and Eq. 37. Hence,
$\delta_{k+1}=\rho_{s+1}(\mathbf{X}^{(k+1)})$ which shows the remaining
statement 5. and concludes the proof of Theorem 2.4.

###  C.2 Proof of Theorem 2.5

1.) Let $\varepsilon,\delta>0$ be arbitrary. Due to the additive structure of
$\mathcal{F}_{\varepsilon,\delta}(\cdot)$, cf. Eq. 4, it is sufficient to
establish that

(38)
$\mathcal{F}_{\textit{sp},\delta}(\mathbf{Z})\leq\mathcal{Q}_{\textit{sp},\delta}(\mathbf{Z}|\mathbf{X})=\mathcal{F}_{\textit{sp},\delta}(\mathbf{X})+\langle\nabla\mathcal{F}_{\textit{sp},\delta}(\mathbf{X}),\mathbf{Z}-\mathbf{X}\rangle+\frac{1}{2}\langle\mathbf{Z}-\mathbf{X},\mathbf{W}_{\mathbf{X},\delta}^{\textit{sp}}(\mathbf{Z}-\mathbf{X})\rangle$

for any $\mathbf{Z},\mathbf{X}\in\mathbb{R}^{n_{1}\times n_{2}}$, where
$\mathbf{W}_{\mathbf{X},\delta}^{\textit{sp}}:\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}^{n_{1}\times n_{2}}$ is defined analogously to Eq. 10 and

(39)
$\mathcal{F}_{\textit{lr},\varepsilon}(\mathbf{Z})\leq\mathcal{Q}_{\textit{lr},\varepsilon}(\mathbf{Z}|\mathbf{X})=\mathcal{F}_{\textit{lr},\varepsilon}(\mathbf{X})+\langle\nabla\mathcal{F}_{\textit{lr},\varepsilon}(\mathbf{X}),\mathbf{Z}-\mathbf{X}\rangle+\frac{1}{2}\langle\mathbf{Z}-\mathbf{X},W_{\mathbf{X},\varepsilon}^{\textit{lr}}(\mathbf{Z}-\mathbf{X})\rangle,$

for any $\mathbf{Z},\mathbf{X}\in\mathbb{R}^{n_{1}\times n_{2}}$, where
$W_{\mathbf{X},\varepsilon}^{\textit{lr}}:\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}^{n_{1}\times n_{2}}$ is defined analogously to Eq. 9.

The argument for Eq. 38 is standard in the IRLS literature [AH15, ODBP15,
PKV23] and is based on the facts that both
$\mathcal{Q}_{\textit{sp},\delta}(\mathbf{Z}|\mathbf{X})$ and
$\mathcal{F}_{\textit{sp},\delta}(\mathbf{Z})$ are row-wise separable, and
that $t\mapsto f_{\sqrt{\tau}}(\sqrt{t})$ is concave and therefore majorized
by its linearization: indeed, let $g_{\tau}:\mathbb{R}\to\mathbb{R}$ be such
that

$g_{\tau}(t):=\begin{cases}\frac{1}{2}\tau\log(e|t|/\tau),&\text{ if
}|t|>\tau,\\\ \frac{1}{2}|t|,&\text{ if }|t|\leq\tau.\end{cases}$

The function $g_{\tau}(\cdot)$ is continuously differentiable with derivative
$g^{\prime}_{\tau}(t)=\frac{\tau}{2\max(|t|,\tau)}\mathrm{sign}(t)$ and
furthermore, concave restricted to the non-negative domain $\mathbb{R}_{\geq
0}$.

Therefore, it holds for any $t,t^{\prime}\in\mathbb{R}_{\geq 0}$ that

$\displaystyle g_{\tau}(t)\leq
g_{\tau}(t^{\prime})+g^{\prime}_{\tau}(t^{\prime})(t-t^{\prime}).$

We recall the definition
$f_{\tau}(t)=\frac{1}{2}\tau^{2}\log(et^{2}/\tau^{2})$ for $|t|>\tau$ and
$f_{\tau}(t)=\frac{1}{2}t^{2}$ for $|t|\leq\tau$ from Eq. 2 with derivative
$f_{\tau}^{\prime}(t)=\frac{\max(t^{2},\tau^{2})t}{t^{2}}=\frac{\tau^{2}t}{\max(t^{2},\tau^{2})}$.
Thus, for any $x,z\in\mathbb{R}$, it follows that

$\begin{split}f_{\tau}(z)=g_{\tau^{2}}(z^{2})&\leq
g_{\tau^{2}}(x^{2})+g^{\prime}_{\tau^{2}}(x^{2})(z^{2}-x^{2})\\\
&=f_{\tau}(x)+\frac{\tau^{2}}{2\max(x^{2},\tau^{2})}(z^{2}-x^{2}),\end{split}$

and inserting $\tau=\delta$, $z=\|\mathbf{Z}_{i,:}\|_{2}$,
$x=\|\mathbf{X}_{i,:}\|_{2}$ and summing over $i=1,\ldots n_{1}$ implies that

$\begin{split}\mathcal{F}_{\textit{sp},\delta}(\mathbf{Z})&=\sum_{i=1}^{n_{1}}f_{\delta}(\|\mathbf{Z}_{i,:}\|_{2})\leq\mathcal{F}_{\textit{sp},\delta}(\mathbf{X})+\sum_{i=1}^{n_{1}}\frac{\delta^{2}}{2\max(\|\mathbf{X}_{i,:}\|_{2}^{2},\delta^{2})}(\|\mathbf{Z}_{i,:}\|_{2}^{2}-\|\mathbf{X}_{i,:}\|_{2}^{2})\\\
&=\mathcal{F}_{\textit{sp},\delta}(\mathbf{X})+\sum_{i=1}^{n_{1}}\frac{\delta^{2}}{\max(\|\mathbf{X}_{i,:}\|_{2}^{2},\delta^{2})}\langle\mathbf{X}_{i,:},\mathbf{Z}_{i,:}-\mathbf{X}_{i,:}\rangle+\frac{1}{2}\sum_{i=1}^{n_{1}}\frac{\big{\|}\mathbf{Z}_{i,:}-\mathbf{X}_{i,:}\big{\|}_{2}^{2}}{\max(\|\mathbf{X}_{i,:}\|_{2}^{2}/\delta^{2},1)}\end{split}$

From the chain rule, it follows that for all $i=1,\ldots,n_{1}$ for which
$\mathbf{X}_{i,:}\neq 0$,

(40)
$\frac{d}{d\mathbf{X}_{i,:}}f_{\delta}(\|\mathbf{X}_{i,:}\|_{2})=f_{\delta}^{\prime}(\|\mathbf{X}_{i,:}\|_{2})\frac{d\|\mathbf{X}_{i,:}\|_{2}}{d\mathbf{X}_{i,:}}=\frac{\delta^{2}\|\mathbf{X}_{i,}\|_{2}}{\max(\|\mathbf{X}_{i,:}\|_{2}^{2},\delta^{2})}\frac{\mathbf{X}_{i,:}}{\|\mathbf{X}_{i,:}\|_{2}},=\frac{\delta^{2}\mathbf{X}_{i,:}}{\max(\|\mathbf{X}_{i,:}\|_{2}^{2},\delta^{2})}$

and therefore

$\mathcal{F}_{\textit{sp},\delta}(\mathbf{Z})\leq\mathcal{F}_{\textit{sp},\delta}(\mathbf{X})+\langle\nabla\mathcal{F}_{\textit{sp},\delta}(\mathbf{X}),\mathbf{Z}-\mathbf{X}\rangle+\frac{1}{2}\langle\mathbf{Z}-\mathbf{X},\mathbf{W}_{\mathbf{X},\delta}^{\textit{sp}}(\mathbf{Z}-\mathbf{X})\rangle$

which shows the majorization of Eq. 38, recalling the definition

$\mathbf{W}_{\mathbf{X},\delta}^{\textit{sp}}=\mathrm{diag}\left(\max\big{(}\|\mathbf{X}_{i,:}\|_{2}^{2}/\delta^{2},1)_{i=1}^{d_{1}}\big{)}^{-1}\right)$

of Eq. 10.

The majorization of Eq. 39 is non-trivial but follows in a straightforward way
from [Küm19, Theorem 2.4] as the objective
$\mathcal{F}_{\textit{lr},\varepsilon}(\mathbf{Z})$ corresponds to the one of
[Küm19, Theorem 2.4] up to a multiplicative factor of $\varepsilon^{2}$ and
constant additive factors, and since the weight operator
$W_{\mathbf{X},\varepsilon}^{\textit{lr}}$ corresponds to the weight operator
used in [Küm19, Chapter 2] for $p=0$.

2.) Due to the definition Eq. 3 of
$\mathcal{F}_{\textit{sp},\delta_{k}}(\cdot)$ and the derivative computation
of Eq. 40, we observe that

$\nabla\mathcal{F}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k)})=\mathrm{diag}\left(\left(\max\left(\big{\|}(\mathbf{X}^{(k)})_{i,:}\big{\|}_{2}^{2}/\delta_{k}^{2},1\right)^{-1}\right)_{i=1}^{d_{1}}\right)\mathbf{X}^{(k)}=\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\cdot\mathbf{X}^{(k)},$

comparing the resulting term with the definition of Eq. 10 of
$\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}$. Furthermore, an
analogue equality follows from the the formula

$\nabla\mathcal{F}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)})=\begin{bmatrix}\mathbf{U}&\mathbf{U}_{\perp}\end{bmatrix}\mathrm{diag}\left(\left(\sigma_{i}^{(k)}\max\left((\sigma_{i}^{(k)})^{2}/\varepsilon_{k}^{2},1\right)^{-1}\right)_{i=1}^{d}\right)\begin{bmatrix}\mathbf{V}^{*}\\\
\mathbf{V}_{\perp}^{*}\end{bmatrix}$

with $\sigma_{i}^{(k)}=\sigma_{i}(\mathbf{X}^{(k)})$ for any $i\leq d$, which
is a direct consequence from the calculus of spectral functions Lemma D.1, and
inserting into the low-rank promoting weight operator formula Eq. 9

$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{X}^{(k)})=\begin{bmatrix}\mathbf{U}&\mathbf{U}_{\perp}\end{bmatrix}\boldsymbol{\Sigma}_{\varepsilon_{k}}^{-1}\mathrm{diag}\left(\left(\sigma_{i}^{(k)}\right)_{i=1}^{d}\right)\boldsymbol{\Sigma}_{\varepsilon_{k}}^{-1}\begin{bmatrix}\mathbf{V}^{*}\\\
\mathbf{V}_{\perp}^{*}\end{bmatrix}=\nabla\mathcal{F}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)})$

Inserting
$\nabla\mathcal{F}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k)})=\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\cdot\mathbf{X}^{(k)}$
and
$\nabla\mathcal{F}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)})=W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{X}^{(k)})$
into the definitions of
$\mathcal{Q}_{\textit{lr},\varepsilon_{k}}(\mathbf{Z}|\mathbf{X}^{(k)})$ and
$\mathcal{Q}_{\textit{sp},\delta_{k}}(\mathbf{Z}|\mathbf{X}^{(k)})$, we see
that it holds that

$\mathcal{Q}_{\textit{lr},\varepsilon_{k}}(\mathbf{Z}|\mathbf{X}^{(k)})=\mathcal{F}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)})+\frac{1}{2}\left(\langle\mathbf{Z},W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{Z})\rangle-\langle\mathbf{X}^{(k)},W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{X}^{(k)})\rangle\right)$

and

$\mathcal{Q}_{\textit{sp},\delta_{k}}(\mathbf{Z}|\mathbf{X}^{(k)})=\mathcal{F}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k)})+\frac{1}{2}\left(\langle\mathbf{Z},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbf{Z}\rangle-\langle\mathbf{X}^{(k)},\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\mathbf{X}^{(k)}\rangle\right).$

Therefore, we see that the weighted least squares solution
$\mathbf{X}^{(k+1)}$ of Eq. 11 for $k+1$ coincides with the minimizer of

(41)
$\begin{split}&\min\limits_{\mathbf{Z}:\mathcal{A}(\mathbf{Z})=\mathbf{y}}\left[\mathcal{Q}_{\textit{lr},\varepsilon_{k}}(\mathbf{Z}|\mathbf{X}^{(k)})+\mathcal{Q}_{\textit{sp},\delta_{k}}(\mathbf{Z}|\mathbf{X}^{(k)})\right]\\\
&=\min\limits_{\mathbf{Z}:\mathcal{A}(\mathbf{Z})=\mathbf{y}}\bigg{[}\mathcal{F}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)})+\mathcal{F}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k)})\\\
&\;\;\;\;\;\;\quad\quad\quad\quad+\frac{1}{2}\left(\langle\mathbf{Z},W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}(\mathbf{Z})\rangle-\langle\mathbf{X}^{(k)},W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k)})\rangle\right)\bigg{]}\end{split}$

with the weight operator $W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}$ of
Eq. 8, which implies that

(42)
$\mathcal{Q}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k+1)}|\mathbf{X}^{(k)})+\mathcal{Q}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k+1)}|\mathbf{X}^{(k)})\leq\mathcal{Q}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)}|\mathbf{X}^{(k)})+\mathcal{Q}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k)}|\mathbf{X}^{(k)}).$

Using the majorization Eq. 17 established in Statement 1 of Theorem 2.5 and
Eq. 42, it follows that

(43)
$\displaystyle\begin{split}\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k+1)})&\leq\mathcal{Q}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k+1)}|\mathbf{X}^{(k)})+\mathcal{Q}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k+1)}|\mathbf{X}^{(k)})\\\
&\leq\mathcal{Q}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)}|\mathbf{X}^{(k)})+\mathcal{Q}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k)}|\mathbf{X}^{(k)})\\\
&=\mathcal{F}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)})+\mathcal{F}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k)})=\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k)}),\end{split}$

using in the third line that
$\mathcal{Q}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)}|\mathbf{X}^{(k)})=\mathcal{F}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)})$
and
$\mathcal{Q}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k)}|\mathbf{X}^{(k)})=\mathcal{F}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k)})$.

To conclude, it suffices to show that
$\varepsilon\mapsto\mathcal{F}_{\varepsilon,\delta_{k}}(\mathbf{X}^{(k+1)})$
and $\delta\mapsto\mathcal{F}_{\varepsilon_{k},\delta}(\mathbf{X}^{(k+1)})$
are non-decreasing functions, since Eq. 43 then extends to

$\displaystyle\mathcal{F}_{\varepsilon_{k+1},\delta_{k+1}}(\mathbf{X}^{(k+1)})\leq\mathcal{F}_{\varepsilon_{k},\delta_{k+1}}(\mathbf{X}^{(k+1)})\leq\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k+1)})\leq\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k)}),$

where we used that the sequences $\varepsilon_{k}$ and $\delta_{k}$ defined in
Algorithm 1 are decreasing. So let us prove this last claim. We define for
$t\in\mathbb{R}$ the function $h_{t}:\mathbb{R}_{>0}\to\mathbb{R}$ such that
$h_{t}(\tau)=f_{\tau}(t)$, i.e.,

$h_{t}(\tau)=\begin{cases}\frac{1}{2}t^{2},&\text{ if }\tau\geq|t|,\\\
\frac{1}{2}\tau^{2}\log(et^{2}/\tau^{2}),&\text{ if }\tau<|t|.\end{cases}$

This function is continuously differentiable with $h_{t}^{\prime}(\tau)=0$ for
all $\tau>|t|$ and

$h_{t}^{\prime}(\tau)=\tau\left(\log(et^{2}/\tau^{2})-1\right)$

for $\tau<|t|$, which implies that $h_{t}^{\prime}(\tau)\geq 0$ for all
$\tau\geq 0$ and thus shows that
$\varepsilon\mapsto\mathcal{F}_{\varepsilon,\delta_{k}}(\mathbf{X}^{(k+1)})$
and $\delta\mapsto\mathcal{F}_{\varepsilon_{k},\delta}(\mathbf{X}^{(k+1)})$
are non-decreasing functions due to the additive structure of
$\mathcal{F}_{\varepsilon,\delta}(\mathbf{X}^{(k+1)})$ and Eq. 4.

3.) First, we argue that $(\mathbf{X}^{(k)})_{k\geq 1}$ is a bounded sequence:
Indeed, if $\overline{\varepsilon}:=\lim_{k\to\infty}\varepsilon_{k}>0$ and
$\overline{\delta}:=\lim_{k\to\infty}\delta_{k}>0$, we note that

$\begin{split}&\frac{1}{2}\overline{\varepsilon}^{2}\log(e\|\mathbf{X}^{(k)}\|^{2}/\overline{\varepsilon}^{2})+\frac{1}{2}\overline{\delta}^{2}\log(e\max_{i}\|\mathbf{X}^{(k)}\|_{\infty,2}/\overline{\delta}^{2})\\\
&=\frac{1}{2}\overline{\varepsilon}^{2}\log(e\sigma_{1}^{2}(\mathbf{X}^{(k)})/\overline{\varepsilon}^{2})+\frac{1}{2}\overline{\delta}^{2}\log(e\max_{i}\|\mathbf{X}^{(k)}_{i,:}\|_{2}/\overline{\delta}^{2})\\\
&\leq\frac{1}{2}\varepsilon_{k}^{2}\log(e\sigma_{1}^{2}(\mathbf{X}^{(k)})/\varepsilon_{k}^{2})+\frac{1}{2}\delta_{k}^{2}\log(e\max_{i}\|\mathbf{X}^{(k)}_{i,:}\|_{2}/\delta_{k}^{2})\\\
&\leq\mathcal{F}_{\textit{lr},\varepsilon_{k}}(\mathbf{X}^{(k)})+\mathcal{F}_{\textit{sp},\delta_{k}}(\mathbf{X}^{(k)})=\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k)})\leq\mathcal{F}_{\varepsilon_{1},\delta_{1}}(\mathbf{X}^{(1)})\\\
&\leq\frac{1}{2}\min(d_{1},d_{2})\sigma_{1}^{2}(\mathbf{X}^{(1)})+\frac{1}{2}d_{1}\max_{i}\|\mathbf{X}^{(1)}_{i,:}\|_{2}^{2}=:C_{\mathbf{X}^{(1)}},\end{split}$

which implies that $\\{\|\mathbf{X}^{(k)}\|\\}_{k\geq 1}$ is bounded by a
constant that depends on $C_{\mathbf{X}^{(1)}}$.

Furthermore, we note that the optimality condition of Eq. 11 (see Lemma C.6)
implies that $\mathbf{X}^{(k+1)}$ satisfies

$\langle
W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k+1)}),\boldsymbol{\Xi}\rangle=0\text{
for all }\boldsymbol{\Xi}\in\ker\mathcal{A}\text{ and
}\mathcal{A}(\mathbf{X}^{(k+1)})=\mathbf{y}.$

Choosing $\boldsymbol{\Xi}=\mathbf{X}^{(k+1)}-\mathbf{X}^{(k)}$ and using the
notation $W^{(k)}=W_{\mathbf{X}^{(k)},\varepsilon_{k},\delta_{k}}$ we see that

(44)
$\begin{split}&\langle\mathbf{X}^{(k+1)},W^{(k)}(\mathbf{X}^{(k+1)})\rangle-\langle\mathbf{X}^{(k)},W^{(k)}(\mathbf{X}^{(k)})\rangle\\\
&=\langle\mathbf{X}^{(k+1)},W^{(k)}(\mathbf{X}^{(k+1)})\rangle-\langle\mathbf{X}^{(k)},W^{(k)}(\mathbf{X}^{(k)})\rangle-2\langle
W^{(k)}(\mathbf{X}^{(k+1)}),\mathbf{X}^{(k+1)}-\mathbf{X}^{(k)}\rangle\\\
&=-\left(\langle\mathbf{X}^{(k+1)},W^{(k)}(\mathbf{X}^{(k+1)})\rangle-2\langle
W^{(k)}(\mathbf{X}^{(k)}),\mathbf{X}^{(k+1)}\rangle+\langle\mathbf{X}^{(k)},W^{(k)}(\mathbf{X}^{(k)})\rangle\right)\\\
&=-\big{\langle}(\mathbf{X}^{(k+1)}-\mathbf{X}^{(k)}),W^{(k)}(\mathbf{X}^{(k+1)}-\mathbf{X}^{(k)})\big{\rangle}.\end{split}$

Due to the definition of $W^{(k)}$, we note that its smallest singular value
(interpreted as matrix operator) can be lower bounded by

$\begin{split}\sigma_{\min}(W^{(k)})&\geq\sigma_{\min}\big{(}W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}\big{)}+\sigma_{\min}\big{(}\mathbf{W}_{\mathbf{X}^{(k)},\delta_{k}}^{\textit{sp}}\big{)}\geq\delta_{k}^{2}/\max_{i}\big{\|}\mathbf{X}^{(k)}_{i,:}\big{\|}_{2}^{2}+\varepsilon_{k}^{2}/\sigma_{1}^{2}(\mathbf{X}^{(k)})\\\
&\geq\overline{\delta}^{2}/c_{\text{sp},\mathbf{X}^{(1)}}+\overline{\varepsilon}^{2}/c_{\text{lr},\mathbf{X}^{(1)}},\end{split}$

where $c_{\text{sp},\mathbf{X}^{(1)}}$ and $c_{\text{lr},\mathbf{X}^{(1)}}$
are constants that satisfy
$c_{\text{sp},\mathbf{X}^{(1)}}\leq\overline{\delta}^{2}\exp(C_{\mathbf{X}^{(1)}}/\overline{\delta}^{2}-1)$
and
$c_{\text{lr},\mathbf{X}^{(1)}}\leq\overline{\varepsilon}^{2}\exp(C_{\mathbf{X}^{(1)}}/\overline{\varepsilon}^{2}-1)$.

Combining this with Eq. 44, the monotonicity according to Statement 2 of
Theorem 2.5, and Eq. 41, it follows that

$\begin{split}\mathcal{F}_{\varepsilon_{k},\delta_{k}}(\mathbf{X}^{(k)})-\mathcal{F}_{\varepsilon_{k+1},\delta_{k+1}}(\mathbf{X}^{(k+1)})&\geq\frac{1}{2}\big{\langle}(\mathbf{X}^{(k)}-\mathbf{X}^{(k+1)}),W^{(k)}(\mathbf{X}^{(k)}-\mathbf{X}^{(k+1)})\big{\rangle}\\\
&\geq\frac{1}{2}\left(\overline{\delta}^{2}/c_{\text{sp},\mathbf{X}^{(1)}}+\overline{\varepsilon}^{2}/c_{\text{lr},\mathbf{X}^{(1)}}\right)\|\mathbf{X}^{(k+1)}-\mathbf{X}^{(k)}\|_{F}^{2}.\end{split}$

Summing over all $k$, this implies that
$\lim_{k\to\infty}\|\mathbf{X}^{(k+1)}-\mathbf{X}^{(k)}\|_{F}=0$.

Since $(\mathbf{X}^{(k)})_{k\geq 1}$ is bounded, each subsequence of
$(\mathbf{X}^{(k)})_{k\geq 1}$ has a convergent subsequence. Let
$(\mathbf{X}^{(k_{\ell})})_{\ell\geq 1}$ be such a sequence with
$\lim_{\ell\to\infty}\mathbf{X}^{(k_{\ell})}=\bar{\mathbf{X}}$, i.e.,
$\bar{\mathbf{X}}$ is an accumulation point of the sequence. As the weight
operator $W^{(k_{\ell})}$ depends continuously on $\mathbf{X}^{(k_{\ell})}$,
there exists a weight operator $\bar{W}:\mathbb{R}^{n_{1}\times
n_{2}}\to\mathbb{R}^{n_{1}\times n_{2}}$ such that
$\bar{W}=\lim_{\ell\to\infty}W^{(k_{\ell})}$.

Since $\lim_{k\to\infty}\|\mathbf{X}^{(k+1)}-\mathbf{X}^{(k)}\|_{F}=0$, it
also holds that $\mathbf{X}^{(k_{\ell}+1)}\to\bar{\mathbf{X}}$ and therefore

$\langle\nabla\mathcal{F}_{\overline{\varepsilon},\overline{\delta}}(\bar{\mathbf{X}}),\boldsymbol{\Xi}\rangle=\langle\bar{W}(\bar{\mathbf{X}}),\boldsymbol{\Xi}\rangle=\lim_{\ell\to\infty}\langle
W^{(k_{\ell})}(\mathbf{X}^{(k_{\ell}+1)}),\boldsymbol{\Xi}\rangle=0$

for all $\boldsymbol{\Xi}\in\ker\mathcal{A}$. The statement is shown as this
is equivalent to $\bar{\mathbf{X}}$ being a stationary point of
$\mathcal{F}_{\overline{\varepsilon},\overline{\delta}}(\cdot)$ subject to the
linear constraint $\\{\mathbf{Z}\in\mathbb{R}^{n_{1}\times
n_{2}}:\mathcal{A}(\mathbf{Z})=\mathbf{y}\\}$.

##  Appendix D Technical addendum

###  D.1 Auxiliary Results

In the proof of Theorem 2.5, we use the following result about the calculus of
spectral functions.

######  Lemma D.1 ([LS05],[FRW11, Proposition 7.4]).

Let $F:\mathbb{R}^{d_{1}\times d_{2}}\to\mathbb{R}$ be a spectral function
$F=f\circ\sigma$ with an associated function $f:\mathbb{R}^{d}\to\mathbb{R}$
that is absolutely permutation symmetric. Then, $F$ is differentiable at
$\mathbf{X}\in\mathbb{R}^{d_{1}\times d_{2}}$ if and only if $f$ is
differentiable at $\sigma(\mathbf{X})\in\mathbb{R}^{d}$.

In this case, the gradient $\nabla F$ of $F$ at $\mathbf{X}$ is given by

$\nabla F(\mathbf{X})=\mathbf{U}\mathrm{diag}\big{(}\nabla
f(\sigma(\mathbf{X})\big{)}\mathbf{V}^{*}$

if
$\mathbf{X}=\mathbf{U}\mathrm{diag}\big{(}\sigma(\mathbf{X})\big{)}\mathbf{V}^{*}$
for unitary matrices $\mathbf{U}\in\mathbb{R}^{d_{1}\times d_{1}}$ and
$\mathbf{V}\in\mathbb{R}^{d_{2}\times d_{2}}$.666Here, for
$\mathbf{v}\in\mathbb{R}^{\min(d_{1},d_{2})}$,
$\mathrm{diag}(\mathbf{v})\in\mathbb{R}^{d_{1}\times d_{2}}$ refers to the
matrix with diagonal elements $\mathbf{v}_{i}$ on its main diagonal and zeros
elsewhere.

###  D.2 Proof of Lemma C.1

The projection operators for $\mathcal{M}_{r}^{n_{1},n_{2}}$,
$\mathcal{N}_{s}^{n_{1},n_{2}}$, and $T_{\mathbf{U},\mathbf{V}}$ are well-
known, see e.g. [Bou23]. To see the final statement assume that $\mathbf{U}$
has row-support $S$ and note that $\mathbb{P}_{\mathbf{U},\mathbf{V},S}$ is
idempotent, i.e.,

$\displaystyle\mathbb{P}_{\mathbf{U},\mathbf{V},S}\mathbb{P}_{\mathbf{U},\mathbf{V},S}\mathbf{Z}$
$\displaystyle=\mathbf{U}\mathbf{U}^{*}(\mathbf{U}\mathbf{U}^{*}\mathbf{Z}+\mathbb{P}_{S}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}-\mathbf{U}\mathbf{U}^{*}\mathbf{Z}\mathbf{V}\mathbf{V}^{*})+\mathbb{P}_{S}(\mathbf{U}\mathbf{U}^{*}\mathbf{Z}+\mathbb{P}_{S}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}-\mathbf{U}\mathbf{U}^{*}\mathbf{Z}\mathbf{V}\mathbf{V}^{*})\mathbf{V}\mathbf{V}^{*}$
$\displaystyle-\mathbf{U}\mathbf{U}^{*}(\mathbf{U}\mathbf{U}^{*}\mathbf{Z}+\mathbb{P}_{S}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}-\mathbf{U}\mathbf{U}^{*}\mathbf{Z}\mathbf{V}\mathbf{V}^{*})\mathbf{V}\mathbf{V}^{*}$
$\displaystyle=\mathbf{U}\mathbf{U}^{*}\mathbf{Z}+\mathbf{U}\mathbf{U}^{*}\mathbb{P}_{S}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}-\mathbf{U}\mathbf{U}^{*}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}+\mathbf{U}\mathbf{U}^{*}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}+\mathbb{P}_{S}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}-\mathbf{U}\mathbf{U}^{*}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}$
$\displaystyle-\mathbf{U}\mathbf{U}^{*}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}-\mathbf{U}\mathbf{U}^{*}\mathbb{P}_{S}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}+\mathbf{U}\mathbf{U}^{*}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}$
$\displaystyle=\mathbf{U}\mathbf{U}^{*}\mathbf{Z}+\mathbb{P}_{S}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}-\mathbf{U}\mathbf{U}^{*}\mathbf{Z}\mathbf{V}\mathbf{V}^{*}=\mathbb{P}_{\mathbf{U},\mathbf{V},S}\mathbf{Z}.$

One can easily check that $\mathbb{P}_{\mathbf{U},\mathbf{V},S}$ acts as
identity when applied to matrices in $T_{\mathbf{U},\mathbf{V},S}$ and that
$\mathbb{P}_{\mathbf{U},\mathbf{V},S}=\mathbb{P}_{\mathbf{U},\mathbf{V},S}^{*}$
since
$\langle\mathbf{Z}^{\prime},\mathbb{P}_{\mathbf{U},\mathbf{V},S}\mathbf{Z}\rangle_{F}=\langle\mathbb{P}_{\mathbf{U},\mathbf{V},S}\mathbf{Z}^{\prime},\mathbf{Z}\rangle_{F}$,
for any $\mathbf{Z},\mathbf{Z}^{\prime}$. This proves the claim.

###  D.3 Proof of Lemma C.4

Let $\boldsymbol{\Xi}\in\ker(\mathcal{A})$. Note that

$\displaystyle
0=\left\|\mathcal{A}(\boldsymbol{\Xi})\right\|_{2}=\left\|\mathcal{A}(\mathbb{P}_{T_{\mathbf{U},\mathbf{V},S}}(\boldsymbol{\Xi})+\mathbb{P}_{T_{\mathbf{U},\mathbf{V},S}}^{\perp}(\boldsymbol{\Xi}))\right\|_{2}\geq\left\|\mathcal{A}(\mathbb{P}_{T_{\mathbf{U},\mathbf{V},S}}(\boldsymbol{\Xi}))\right\|_{2}-\left\|\mathcal{A}(\mathbb{P}_{T_{\mathbf{U},\mathbf{V},S}}^{\perp}(\boldsymbol{\Xi}))\right\|_{2}$

By the RIP we hence get that

$\displaystyle\left\|\mathbb{P}_{T_{\mathbf{U},\mathbf{V},S}}(\boldsymbol{\Xi})\right\|_{F}^{2}$
$\displaystyle\leq\frac{1}{1-\delta}\left\|\mathcal{A}(\mathbb{P}_{T_{\mathbf{U},\mathbf{V},S}}(\boldsymbol{\Xi}))\right\|_{2}^{2}\leq\frac{1}{1-\delta}\left\|\mathcal{A}(\mathbb{P}_{T_{\mathbf{U},\mathbf{V},S}}^{\perp}(\boldsymbol{\Xi}))\right\|_{2}^{2}$
$\displaystyle\leq\frac{\left\|\mathcal{A}\right\|_{2\to
2}^{2}}{(1-\delta)}\left\|\mathbb{P}_{T_{\mathbf{U},\mathbf{V},S}}^{\perp}(\boldsymbol{\Xi})\right\|_{F}^{2}.$

Consequently,

$\displaystyle\left\|\boldsymbol{\Xi}\right\|_{F}^{2}=\left\|\mathbb{P}_{T_{\mathbf{U},\mathbf{V},S}}(\boldsymbol{\Xi})\right\|_{F}^{2}+\left\|\mathbb{P}_{T_{\mathbf{U},\mathbf{V},S}}^{\perp}(\boldsymbol{\Xi})\right\|_{F}^{2}\leq\left(1+\frac{\left\|\mathcal{A}\right\|_{2\to
2}^{2}}{(1-\delta)}\right)\left\|\mathbb{P}_{T_{\mathbf{U},\mathbf{V},S}}^{\perp}(\boldsymbol{\Xi})\right\|_{F}^{2}.$

###  D.4 Proof of Lemma C.7

In the proof of Lemma C.7 we will use the following fact.

######  Lemma D.2.

Let $W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}$ be the weight
operator defined in Eq. 9, which is based on the matrices
$\mathbf{U}\in\mathbb{R}^{n_{1}\times r_{k}}$ and
$\mathbf{V}\in\mathbb{R}^{n_{2}\times r_{k}}$ of leading $r_{k}$ left and
right singular vectors of $\mathbf{X}^{(k)}$. If $\mathbf{M}\in
T_{\mathbf{U},\mathbf{V}}$, then
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{M})\in
T_{\mathbf{U},\mathbf{V}}$. If $\mathbf{M}\in
T_{\mathbf{U},\mathbf{V}}^{\perp}$, then
$W_{\mathbf{X}^{(k)},\varepsilon_{k}}^{\textit{lr}}(\mathbf{M})\in
T_{\mathbf{U},\mathbf{V}}^{\perp}$.

######  Proof D.3.

If $\mathbf{M}\in T_{\mathbf{U},\mathbf{V}}$, there exist
$\mathbf{M}_{1}\in\mathbb{R}^{r_{k}\times r_{k}}$,
$\mathbf{M}_{2}\in\mathbb{R}^{r_{k}\times(n_{2}-r_{k})}$,
$\mathbf{M}_{3}\in\mathbb{R}^{(n_{1}-r_{k})\times r_{k}}$ such that

$\mathbf{M}=\begin{bmatrix}\mathbf{U}&\mathbf{U}_{\perp}\end{bmatrix}\begin{bmatrix}\mathbf{M}_{1}&\mathbf{M}_{2}\\\
\mathbf{M}_{3}&0\end{bmatrix}\begin{bmatrix}\mathbf{V}^{*}\\\