Recovering Simultaneously Structured Data via Non-Convex Iteratively Reweighted Least Squares
Authors: Christian Kümmerle
Johannes Maly
Venue: NeurIPS 2023 poster (2023)
Keywords: low-rank models
sparsity
iteratively reweighted least squares
non-convex optimization
quadratic convergence
simultaneously structured data
TL;DR: We propose a novel iteratively reweighted least squares method for the recovery of simultaneously structured data from measurements, and show its local quadratic convergence in a minimal sample complexity regime.

Abstract: We propose a new algorithm for the problem of recovering data that adheres to multiple, heterogenous low-dimensional structures from linear observations. Focussing 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.

Recovering Simultaneously Structured Data via Non-Convex Iteratively Reweighted Least Squares
Christian Kümmerle
Johannes Maly
Department of Computer Science University of North Carolina Charlotte
Department of Mathematics Ludwig-Maxings Maximilians- -nniiiisaannniiieiitttt München
Charlotte, NC 28223, USA kuemmerle@uncc.eer
80799 Munich, Germany and Munich Center for Machine Learning (MCML) maly@math 1mu.lmu..... de
Abstract
We propose new algorithm for the problem of recovering data that adheres to multiple, heterogeneous low-dimensional structures from linear observations. Focusing 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 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 simultaneously structured data matrix in a regime minimal sample complexity (up constants and logarithmic factor), which is known to be impossible for 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-ar methods.
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 [7]. If the latent image is n-dimensional, for N. it is well-known that ((n) observations are required for robust identification in general. In practice, imaging problems are however often ill- posed, i.e., the number of observations smaller than n or the operator creating the observations is defective [88, 93]. In such situations, the fundamental lower bound of (2)) 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, sparsity 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 [10, 68]. For instance, consider the recovery of an nı ni-dieensional image X, g xa2 from linear observations
y A(X,)
R.
where A: R x - RTh is linear operator modeling the impulse response the sensing device and 77 E m models additive noise. Whereas this problem is ill-posed for n 122, it has been established [38, 15, 78] that it becomes well-posed X, is sparse or of low rank. The aforementioned works prove that
[(31sg) observations suffice for robust reconstruction if X, is -oooo-spe
vector x E R" is called s-sparse if x has at most s non-zero entries. For a matrix X E R Xna there are various ways to count the level of sparsity. In this work, we use the most common definition and call X s-row-sparse (resp. -column-sparse) if most rows (resp. columns) non-zero.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

and bg-columm-spea sparse, and that m [(((m1 + n2)) observations suffice if X, is rank-r matrix. These bounds, which relax the general lower bound of m (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 he formulated as instances of (1) with X. being sinultaneously structured, i.c., X, is both of rank r and Examples encompass sparse phase retrieval [51, 46, 11, 47, 85], sparse blind deconvolution [59, 83]. hyperspectral imaging [43, 41, 89, 102], sparse reduced-rank regression [22, 100], and graph denoising and refinement [80, 103]. In these settings, the hope is that due to leveraging the simultaneous structure, [(r(s1 + $2)) observations suffice to identify the data matrix. For r' < $1, $2 < n1. n2., these bounds are significantly smaller than the bounds for single-structured data matrices.
8i heop--ppareesy--couun---ps
From an algorithmic point of view, however, the simultaneously structured recovery problem poses obstacles that are not present for problems where X, is only of low-rank or (group) sparse: In the latter case, variational methods [6] that formulate the reconstruction method in terms of optimizing a suitable objective with a structural regularization term involving Schatten (quasi-)norms have been well-understood, leading tractable algorithms in the information theoretically optimal regime [28, 17].
and and Sp-
For simultaneously structured problems, on the other hand, Oymak al. showed in [74] that mere linear combination of convex regularizers for different sparsity structures in our case, nuclear and E2,1 12,1-nomss -nomms 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 single prior; was also shown in [74] that minimizing a linear combination of rank and row sparsity can indeed lead to guaranteed recovery from D(r(51 + 82)) measurements. The downside is that the combination of these non-conve 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 [74] 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) [26, 35, 70, 5, 1, 54, 55], the convergence of which we analyze in the information theoretically (near-)optimal regime. Specifically, our main contributions are threefold:
(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 (12). Notably, the algorithm automatically balances between its low-rank and its sparsity promoting terms, leading to a reliable identification of s--row-sparse and rank-r ground truths?.
(ii) Under the assumption that A behaves almost isometrically on the set of row-sparse and low-rank matrices, we show in Theorem 2.5 that locally Algorithm 1 exhibits quadratic convergence towards X.. Note that if A is, e.g., a Gaussian operator, the isometry as- sumption is (up log factors) fulfilled the information theoretic (near- hoptimal regime
= =[(r(s) + v22)) [60].
(iii) Finally, in Section 2.3 we identify the underlying family of objectives that are minimized by Algorithm 1. To make this precise, define for and denoting Fuler's number the real-valued function fr R + I such that
1 7 log(et?
if i > T, if f S T,
50(1)
(2)
"To enhance readability of the presented proof' and results, we restrict ourselves to row-sparsity of X. here. It is straight- forward to generalize the arguments to column-sparsity well.

---

which is quadratic around the origin and otherwise logarithmic in its argument. Using this definition, we define for E > the (E-)smoothed log-determnnant objective Fire IR - R and for o > 0 the -smmototh sum of logarithmic row-wise fn-norms objective Fsp.8 : R XR2 + R such that
min(nimme)
11 Σ fo(()x
(x) =
Σ
(x(),
F
(3)
Combining the above, we further define the (e, smoothed logarithmic surrogate objective Fe.8 : R X + - as
FX) == Fir, (X) ( Fop (X))
(4)
In Theorem 2.6, we prove for any A that the iterates of Algorithm 1 minimize quadratic majorizations of Fe.s and form a non-increasing sequence on Fe,s. 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 rigorous variational interpretation.
The numerical simulations in Section 4 support our theoretical findings and provide empirical evidence for the efficacy of the proposed method.
1.2 Related Work
Sparse and Low-Rank Recovery, Whereas leveraging a single matrix structure like sparsity or low-rankness in the reconstruction process easily be obtained by convex regularizers [78, 16], Oymak et al. [74] showed that, if one is interested in near-optimal sampling rates, one cannot expect comparably simple solutions for identifying simultaneously structured objects; minimization of (4) with convex terms (-) and F.).((-) would be only as good as using the one structure that is information theoretically more favorable. A closely related problem' that appears in statistical literature under the name Sparse Principal Component Analysis (SPCA) [104, 25] is known to be NP-hard in general [67]. Despite the the intrinsic hardness of simultaneously structured recovery problems, promising empirical results for hyperspectral image demixing were shown [40], where minimization problems involving the sum of reweighted convex surrogates are solved by proximal scheme based on ADMM for the case of simultaneously sparse, low-rank and non-negative matrices. For the problem of simultaneously sparse and low-rank matrix recovery, there exist only a handful approaches that come with rigorous theoretical analysis. The first line of works [4, 37] aims to overcome the aforementioned limitations of purely convex methods in neat way. They assume that the operator A has a nested structure such that basic solvers for low-rank resp. row/column-sparc recovery can be applied in two consecutive steps. Despite being an elegant idea, this approach clearly restricts possible choices for A and is of hardly any practical use.
In a second line of work, Lee et al. [60] consider general impulse response operators that satisfy a suitable restricted isometry property for -oow and 8 s2-column-sparss 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 [48] and uses hard thresholding pursuit [36] to enforce sparsity in addition. In particular, they show that if X, is rank-R, has s1-sparse columns and sp-sparse rows, then m ~ R(s1 + 82) (oggmax((en) (exxx(en///// $1, enz/ s2 }) 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 X, 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 [59]. In [100], related approach that combines gradient descent of a smooth objective with hard thresholding considered, for which the authors show linear convergence from suitable intialization if the measurement operator satisfies a restricted strong convexity and smoothness assumption.
third line of work, approaches the problem from a variational point of view. In [33, 69] the
*Since 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.
3

---

authors aim at enhancing robustness of recovery by alternating minimization of an 4 no-norm based multi-penally 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 [60].
The works [37, 29I, which build upon generalized projection operators to modify iterative hard thresholding to the simultaneous setting, share the lack of global convergence guarantees.
In [79], the authors examine the use of atomic norms to perform recovery of simultaneously sparse and low-rank matrices, which uses a related, but different sparsity assumption compared to the row or column sparsity studied here. From a practical point of view, the such norms are hard to compute and the paper only proposes a heuristic polynomial time algorithm for the problem. Finally, the alternative approach of using optimally weighted sums maxima of convex regularizers [52] 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 long line of research on IRLS going back to Weiszfeld's ss algorithm proposed in the 1930s for a facility location problem [95, 5J. 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 [42, 18, 26, 57, 34, 55], robust statistics [45, 2, 72], computer vision [19, 61, 84], low-rank matrix recovery and completion [35, 70, 56, 54], and in inverse problems involving group sparsity [21, 101, 20]. Recently, has been shown [62] 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 [58], optimal performance is achieved if intertwined with a smoothing strategy for the original objective, in which case globally linear (for convex objectives) [26, I, 72, 55, 75] and locally superlinear (for non-convex objectives) [26, 56, 54, 75] convergence rates have been shown under suitable conditions on the linear operator A.
However, there has only been little work on IRLS optimizing sum of heterogenous objectives [81] including the combination of low-rank promoting and sparsity-prone promoting objectives nor on the convergence analysis of any such methods. The sole algorithmically related approach for our setting has been studied in [23], where a method has been derived in sparse Bayesian learning framework, the main step of which amounts to the minimization of weighted least squares problems. Whereas the algorithm of [23] 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 cnarmation-thhortticn
optimal regime. Furthermore, [23] 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 of matrix Z E R by Z;,; and the j-th column of Z by Z..g- We abbreviate n = min(n), [11, n22)- We denote the r-th singular value of a matrix Z E m X by (Z). Likewise, we denote the in nn-mrmm s-largest row of Z by (Z). (To determine Pa (Z), we form the in fn-morm non-increasing rearrangement of the rows of Z and, by convention, sort rows with equal norm according to the row-index.) We use denote the Hadamard (or Schur) product, i.e., the entry-wise product of two vectors/matrices. We denote the Euclidean lg- fo-norm of a vector z R" by
||z For Z E Rm1x"2, the matrix norms we use encompass the operator norm || |Z|| := sup/|w|[a= 1= /Zw|22 the row-sum p-quasinorm
|Z|pp2 (E IZ. 12) 1/p the row-max norm
| || xx. = maxte[e.) Z2, and the
Schatten-p quasinorm [[ls, (Er= o, (z)))))'. Note that two special cases of Schatten quasinorms are the nuclear norm ||]||
:= ||]]|, and the Frobenius norm ||]| := |/]||s.

---

IRLS for Sparse and Low-Rank Reconstruction
Recall that we are interested in recovering rank-r and s-row-sparse matrix X, [2 linear observations
X 1 from m
y A(X)) ER"
(5)
ie., A: R X + T is lincar. We write X, - == M /W where M C R denotes the set of matrices with rank most and N C R denotes the set of matrices with at most s non-zero rows. For convenience, we suppress the indices n1 and 122 whenever the ambient dimension is clear from the context. In particular, know that X, U.E.V; where U. E NNr Σ, E Rexr, and V. & R22xr denote the reduced SVD of X.. Furthermore, the row supports of X, and U, (the index sets of non-zero of X, resp. U.) are identical, i.e., (X.)(X.) = supp(U.) .)= = S. C [mi] = {l....,... R1).
2.1 How to Combine Sparse and Low-Rank Weighting
As discussed in Section 1, the challenge in designing a reconstruction method for (5) lies in simul- taneously leveraging both structures X. in order to achieve the optimal sample complexity of r(s + n2). To this end,
propose novel (RLS-based approach in Algorithm 1. The key in- gredient for computing the (k: + 1)-st iterate ((+1) E [ X in Algorithm is the multi strti-structural weight operator W Wx.+).e, R Xx2 + R x of Definition 2.1, which depends on the current
iterate (A)
Definition 2.1. For db, En 0, of = (x()), and ((k) E g Xna, let
Th = ii [n]] o; (k) Ex}
(6)
denote the number of singular values o (o, (k) T4 ) of X (k) larger than ek and furthermore be
Sk := [mi] | 112 5R)
the number of rows of X with E2-norm larger than dy- Define the left and right singular vectors of (k) as columns of E Hmix XTa and E 1R=2XT respectively, corresponding to the
leading singular values o (k) = (o,")) = (x))) )=
We define the weight operator W (1) 84.88 8x 1Ra1X=2 + guixna at iteration k of Algorithm I as
W ( 5u (Z) = W (Z)
W X(0
Z
(8)
where W
R Xna
+ R Xna is its low-rank promoting part
U
V
Whrr X( (Z) = U UIE
Z
V1
(9)
U
Ek
Vi
with Eg, (a) /EKs , / 1) and Wa
X(0, E TRRixns xn1 its pparsiy-prooootinn part, which is
diagonal with
( x
X (2) max (I.
for all [ni].
(10)
The matrices U. and V are arbitrary complementary orthogonal bases for and V that do do not need to be computed in Algorithm 1.
Remark 2.2. For the sake of conciseness, we only consider row-sparsity here. Algorithm 1 and its analysis can however be modified to cover row- and column-sparse matrices well. For instance, in the symmetric setting X, = =X! (naturally occurring in applications like sparse phase retrieval) one would define the weight operator W Wxt))
sn.52 as (8), but with an additional term that multiplies
wap X &v. from the right Z, which corresponds to minimizing the sum of three smoothed logarithmic surrogates. In this case, the solving modified weighted least squares problem (11) will have similar complexity (potentially smaller complexity, as additional symmetries can be exploited).
5

---

Algorithm IRLS for simultaneously low-rank rand row-sparse matrices
1: Input: Linear operator A: R=1x" - R" data y E R", rank and sparsity estimates 7 and s.
2: Initialize k = 0, W Wx(().ook Id and set 8z. Ex = c0.
3: for k == K do
4: Weighted Least Squares: Update iterate x by
X (k) arg min (X, Wxa Maot--1,cc 1. (X)) X:A(x)=y
(11)
5: Update Smoothing: Compute F + 1-st singular value 0,+1 == 0j+1 00+1 (x)) and (5 + 1)-st
largest row 22-norm (x+1(X(k)) of X (k), update
EK = min Ek-1, 0+11
6k min (0k-1, P5+1
(12)
6: Update Weight Operator: For тı and 8 as in (6) and (7),
compute first rh singular triplets (*)) R U E Rmix* and V E R compute Wx0,,,. in (8) via Wike X in (9) and W X in (10).
7: end for
8: Output: X(K)
Recall Fer,ber Fire and Fsp.ó from (3)-(4). The high-level idea of Algorithm 1, for other IRLS methods, is to minimize quadratic functionals, which call Qir.a (: (x) R x R and
2pp.6 (+) 12 x + R and define them by
Qir.a (z)x = Fire (x) (VFine (x())), x(x)) +
2 7z (k)
W
(Z x())),
(1k
2.p.6. (Z)x()) = Fop.6r (x()) (VF (x(n), (x 2- --x()) + 2 (Z x
w ( on (Z x
(13)
that majorize Fec, 62 (-) (see Theorem 2.6 below) for any iteration k. This minimization leads the weighted least squares problem (11) in Algorithm 1. This step can be implemented by standard numerical linear algebra (see the supplementary material for discussion of its computational complexity). As a second ingredient of the method, the smoothing parameters & and 5 of Fex are updated (i.e., decreased) in step (12) before the weight operator updated according the current iterate information. For the weight operator update, it is only necessary to compute row norms and leading singular triplets of (k)
Remark 2.3. The particular form of the low-rank promoting part of the weight operator W X(t) eew in
(9) is due to [53, 54] and captures optimally spectral information both in the column and row space, unlike prior work on low-rank IRLS [35, 70], while retaining the property that the induced quadratic model Qir.e (x( majorizes Fure, (-) (see proof of Theorem 2.6). This choice is critical to enable fast local rate as established in Theorem 2.5.
2.2 Local Quadratic Convergence of IRLS
Our first main result states that Algorithm 1 exhibits quadratic convergence in a local neighborhood of X.. a property Algorithm 1 shares with several methods from the IRLS family. We only need to assume that A acts almost isometrically on the set My,s-
Definition 2.4. We say that a linear operator A: Ra x - T satisfies the rank-r and row-s- sparse restricted isometry property (or (r, R)-R with RIP-constant 8 E (0, 1) if
(1 (1- 5) | 2
||((Z)] < (1
for all Z Mr,.-
is worth highlighting that Gaussian operators satisfy the above RIP with high-probability vn 2 cr(s + n12) log(en1/s) for some absolute constant e > 0, see for instance [60]. Up to log- factors, this at the information theoretic limit which we discussed in the beginning. The convergence result for Algorithm 1 now reads as follows.

---

Theorem 2.5 (Local Quadratic Convergence). Let X. E E,,,, be a fixed ground-true truth matrix that s-row-sparse and of rank r. Let linear observations y == A(X.) be given and assume that A has the (r, (r.8)-RIP with E (0, 1). Assume that the k-th iterate X of Algorithm with F r and s updates the smoothing parameters in (12) such that one of the statements EK 0r+1 (X()) or
5g Ps+1 (x)))) is true, and that rk 2 T and Sk 2
X (R) satisfies
min or (x((X.) Ps (X.))
1x X.II
(14)
48 niy All2_
where 2-
1
AA (1-6)
and n = min fn), fni n2), then the local convergence rate is quadratic
in the sense that
X (++1) - X.Il S min
0.11
-x.11,
for
5r
28
#1 4.19cc TAll
(15)
(x,
Ps Px(
and X (k++) 2000 X..
The proof of Theorem 2.5 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-arr competitor methods [60, 37, 69, 29] reach a local linear error decay at best.
On the other hand, (14) is rather pessimistic since for Gaussian A the constant C|AAll2 .2 scales like ([a(nn (2)) )/2)) whereas experiments within an empirically much larger convergence radius. Closing this gap between theory and practical performance is future work.
v /mnnaa)/m, which means that the Eight-hand side of (14) behaves like m (/3///(n(nn we observe quadratic convergence
noteworthy that the theory in [60] - to our knowledge the only other related work explicitly characterizing the convergence radius - holds
neighborhood of X. that is independent of the ambient dimension. The authors of [60] however assume that the RIP-constant decays with the conditioning number ε of A, quantity that might be large in applications. Hence, ignoring log-factors the sufficient number of measurements [60] scales like m = 2(x2r(s1 + n2)). In contrast, Theorem 2.5 works for any RIP-constant less than one which means for m [[((s1 + n2)).
2.3 IRLS as Quadratic Majorize-Miniminn Algorithm
With Theorem 2.5, we have provided a local convergence theorem that quantifies the behavior of Algorithm in a small neighbourhood of the simultaneously row-sparse and low-rank ground-truth X. "R"1Xmg Xm2, The result based on sufficient regularity of the measurement operator A, which in turn is satisfied with high probability if A consists of sufficiently generic random linear observations that concentrate around their mean.
In this section we establish that, for any measurement operator A, Algorithm 1 can be interpreted within the framework of iteratively reweighted least squares (IRLS) algorithms [26, 70, 75], which implies strong connection to the minimization of a suitable smoothened objective function. In our case, the objective Fa,s in (4) is a linear combination of sum-of-logarithms terms penalizing both non-zero singular values [31, 70, 14] as well as non-zero rows of matrix X [50].
We show in Theorem 2.6 below that the IRLS algorithm Algorithm studied in this paper is based on minimizing each iteration quadratic models that majorize Fe.s, and furthermore, that the iterates (X 1 of Algorithm define non-increasing sequence (Fe,. (X())) k21 with respect the objective Fs.s of (4). The proof combines the fact that for fixed smoothing parameters en and òx. the weighted least squares and weight update steps Algorithm can be interpreted as step of Majorize-Minimias algorithm [86, 58], with decrease the underlying objective (4) for updated smoothing parameters.
Theorem 2.6. Let E Rm,, let the linear operator A: R - Rm be arbitrary. If (X (k) 1221 a sequence of iterates of Algorithm and (ōx )22 and (E1) are the sequences of smoothing parameters defined therein, then the following statements hold.
7

---

1. The quadratic model functions Qir.es (- (x() and 2ap ()) (xx()) defined in (13) globally majorize the (Ex, (onx:)))smoohhee logarithmic surrogate objective Ferby for any Z R , it holds that
Fen.be (2) S Qir,e (Z)x(k)) + ( (z)x
(16)
2. The sequence (F (X (k)) is non-increasing.
3. I == lim Ek 0 and o lim 6% O, then limg- [X Furthermore, in this case, every accumulation point of (x is stationary point of
x(++)) O.
min X:A(X)=y
(x)
The proof of Theorem 2.6 is presented in the supplementary material.
Discussion of Computational Complexity
It is well-known that the solution of the linearly constrained weighted least squares problem (11) can be written as
( W (AW
(17)
where Wk-1 == W 1. is the weight operator (8) of iteration 1: [26, 54]. In [54, Theorem 3.1 and Supplementary Material), it was shown that in the case of low-rank matrix com- pletion without the presence of row-sitrs inducing term, this weighted least squares problem can be solved by solving an equivalent, well-conditiones linear system via an iterative solver that uses the application of system matrix whose matrix-vector products have time complexity of O(mr + 2 max(n1, 122))
In the of Algorithm 1, the situations is slightly more involved as cannot provide explicit formula for the inverse of the weight operator Wi-1 as it amounts to the sum of the weight operators Wir () and 'P x(k--))6 that are diagonalized by different, mutually incompatible bases.
However, computing this inverse facilitated by the Sherman-
Woodbury formula [97]
(ECF" -)- B B B- + FBB F "F'B'
for suitable matrices of compatible size E, F and invertible C, B and the fact that both W x(k-1),
and W xx*-1),60-1 exhibit "low-rank plus (scaled) identity" or a "sparse diagonal plus (scaled) identity" structure. After simple application of the SMW formula, (17) can be rewritten such that the computational bottleneck becomes the assembly and inversion of a Orr max(n), n2)) linear system. We note that in general, this can be done exactly in time complexity of O(rgmaz(n) n2)²) using standard linear algebra. A crucial factor in the computational cost of the method is also the structure of the measurement operator A defining the problem, as the application of itself and its adjoint can significantly influence the per-iteration cost of IRLS; for dense Gaussian measurements, just processing the information of AA amounts to mning flops. If rank-one Fourier-type measurements are taken, this cost can significantly be reduced, see [29, Table 1 and Section 3] for an analogous discussion.
We refer to the MATLAB implementation available in the repository sttps://githie com/ ckuemmerle/ simirls for further details. While our implementation is not optimized for large-scale problems, the computational cost of Algorithm was observed to be comparable to the implementa- tions of SPF RiemAdaIHT provided by the authors [69, 29]. We leave further improvements and adaptations to large-scale settings to future work.
4 Numerical Evaluation
In this section, we explore the empirical performance of IRLS in view the theoretical results of Theorems 2.5 and 2.6, 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) [60] and Rimmannian adaptive iterative hard thresholding (RiemAdaIHT) [29], which are

---

among the methods with the best empirical performance reported in the literature. The method ATLAS [33] and its successor [69] are not used in our empirical studies since they are tailored to robust recovery and yield suboptimal performance when seeking righ-precision reconstruction in low noise scenarios. We use spectral initialization for SPF and RiemAdaIHT. The weight operator IRLS initialized by the identity as described in Algorithm 1, solving unweighted least squares problem in the first iteration. A detailed description of the experimental setup can be found in Appendix A.
Performance in www-Measurement Regime. Figures and show the empirical probability of successful recovery when recovering s-row sparse ground-truths X, E R 2556x40 of rank 1 (resp. r = 5) from Gaussian measurements under oracle knowledge on and The results are averaged over random trials. As both figures illustrate, the region of success of IRLS comes closest to the information theoretic limit of r(s n2 r) which is highlighted by red line, requiring significantly lower oversampling factor than the baseline methods.
Nr. messinmettss
200
200
200
0.8
150
150
150
0.6
0.4
100
100
100
0.2
50
50
0
10
20
30
40
10
20
30
40
10
20
30
40
sparsity
Figure Left column: RiemAdaIHT, center: SPF, right: IRLS. Phase transition experiments with M1 256, n22 = 40, 1, Gaussian measurements. Algorithmic hyperparameters informed by model order knowledge (i.e., r and s for IRLS). White corresponds empirical success rate of 1, black to 0.
Number musorrctentts
1,000
1,000
1,000
800
800
800
600
600
600
400
400
400
60
20
40
60
20
40
60
1,000
1,000
1 1,000
800
800
800
600
600
600
400
400
400
20
40
60
20
40
60
20
40
60
sparsity
Figure Left column: RiemAdaIHT, center: SPF, right: IRLS. First row: As in Figure 1, but for data matrix X, of rank
5. Second row: As first row, but hyper-parameter, and are overestimated
as 10, 115ss
In Appendix A.2 and Appendix A.3 in the supplementary material, we report on similar experiments conducted for other measurement operators than dense Gaussians, which cases the empirical relative behavior of the methods is comparable.
Sensitivity to Parameter Choice. In applications our setting, the quantities and 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 algorithms with rank and sparsity estimates of 7 2r and 3 [1.5x]. Whereas all considered methods suffer a deterioration of performance, we observe that IRLS deteriorates relatively the least by large margin. Furthermore, we observe that even if IRLS does not recovery X,, it converges typically to matrix that is still low-rank and row-sparse (with

---

larger and s) satisfying the data constraint, while the other methods fail to convergence to such matrix.
Convergence Behavior. Finally, we examine the convergence rate of the iterates to validate the theoretical prediction of Theorem 2.5 in the setting of Figure Figure 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.5 could possibly be improved by future investigations.
10-1
IRLS- RRLS-Lowe oow
Rel Frehenin enn
Sparse Power Factorization [60]
Riemannian Adaptive IHT [29]
10
10
7
14
21
28
35
42
49
56
63
70
77
84
91
98
105
112
Figure 3: Comparison of convergence rate. Setting in Figure 2 with 40 and m 1125.
Further experiments. In Appendix A.4, we provide additional experiments investigating the self- balancing property of the objective (4), well as the experiments on the noise robustness of the method in Appendix A.5.
5 Conclusion, Limitations and Future Work
Conclusion. In this paper, we adapted the IRLS framework the problem of recovering simulta- neously structured matrices from linear observations focusing on the special case of row sparsity and low-rankness. Our convergence guarantee Theorem 2.5 hereby the first one for any method minimizing combinations of structural surrogate objectives that holds in the information-thhoree near-optimal regime and exhibits local quadratic convergence. The numerical experiments we con- ducted for synthetic data suggest that, due to its weak dependence on the choice of hyperparameters, IRLS in the form of Algorithm can be practical method for identifying simultancously structured data even in difficult problem instances.
As in the case of established IRLS methods that optimize non- convex surrogate objectives representing single structure [26, 56, 54, 75], the radius of guaranteed quadratic convergence Theorem 2.5 the most restrictive assumption. Beyond the interpretation in terms of surrogate minimization presented in Theorem 2.6, which holds without any assumptions the initialization, our method shares a lack of global convergence guarantees with other non-convex IRLS algorithms [70, 54, 75].
Limitations and Future Work.
The generalization and application of our framework to combinations of structures beyond rank and (or column-)sparsite lies outside the scope of the present paper, but could involve subspace- structured low-rankness [31, 24, 99, 87] or analysis sparsity [30, 77]. A generalization of the presented IRLS framework to higher-order objects such as low-rank tensors is of future interest as convexifications structure- structure-promoting objectives face similar challenges [71, 74, 98] in this case.
In parameter-efficieen deep learning, both sparse [39, 44, 32, 91] and low-rank [96, 94, 82] weight parameter models have gained considerable attention due to the challenges training and storing, e.g., in large transformer-based models [92, 9]. It will be of interest 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.
10

---

Acknowledgements
The authors thank Max Pfeffer for providing their implementation of the algorithm RiemAdaIHT as considered in [29].
References
[1] Deeksha Adil, Richard Peng, and Sushant Sachdeva. Fast, provably convergent IRLS algorithm for p-norm linear regression. Advances in Neural Information Processing Systems (NeurIPS 20/9), 2019.
[2] Khurrum Aftab and Richard Hartley. Convergence iteratively re-weighted least squares to robust m-estimators. In IEEE Winter Conference on Applications of Computer Vision (WACV 2015), pages 480-487, 2015.
[3] Ali Ahmed, Benjamin Recht, and Justin Romberg, Blind deconvolution using convex program- ming. IEEE Transactions on Information Theory, 60(3):1111-1777 2014.
[4] Sohail Bahmani and Justin Romberg. Near-optimal estimation of simultancously sparse and low-rank matrices from nested linear measurements. Information and Inference: Journal of the IMA, 5(3):331-351, 2016.
[5] Amir Beck and Shoham Sabach. Weiszfeld's method: Old and new results. Journal of Optimization Theory and Applications, 164(1):1-40, 2015.
[6] Martin Benning and Martin Burger. Modern regularization methods for inverse problems. Acta Numerica, 27:1-111, 2018.
[7] Mario Bertero, Patrizia Boccacci, and Christine De Mol. Introduction to inverse problems in imaging. CRC Press, 2021.
[8] Nicolas Boumal. An introduction to optimization on smooth manifolds. Cambridge University Press, 2023.
[9] Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhari- wal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al. Language models are few-shot learners. Advances in Neural Information Processing Systems (NeurIPS 2020), 33: 877-11991,, 2020.
[10] Alfred M. Bruckstein, David L. Donoho, and Michael Elad. From sparse solutions of systems of equations to sparse modeling of signals and images. SIAM Review, 51(1):34 81, 2009.
[11] T Tony Cai, Xiaodong Li, and Zongming Ma. Optimal rates of convergence for noisy sparse phase retrieval via thresholded wirtinger flow. The Annals of Statistics, pages 2221-2251. 2016.
[12] T Tony Cai and Anru Zhang. Rate-optimal perturbation bounds for singular subspaces with applications high-dimensional statistics. The Annals of Statistics, 46(1):60-89, 2018.
[13] Patrizio Campisi and Karen Egiazarian. Blind Image Deconvolution: Theory and Applications. CRC Press, 2017.
[14] F. J. Candès, Y. Eldar, T. Strohmer, and V. Voroninski. Phase Retrieval via Matrix Completion. 2013.
SIAM J. Imag. Sci., 6(1):199-223 2013.
[15] 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):2333 5744)2342-23355 2011.
[16] Venkat Chandrasekaran, Benjamin Recht, Pablo A Parrilo, and Alan S Willsky. The convex geometry of linear inverse problems. Foundations of Computational Mathematics, 12(6):805-
849, 2012.

---

[17] Rick Chartrand. Exact reconstruction of sparse signals via nonconvex minimization. IEEE Signal Processing Letters, 14(10):70-770, 2007.
[18] Rick Chartrand and Wotao Yin. Iteratively reweighted algorithms for compressive sensing, In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3869-3872, 2008.
[19] Avishek Chatterjee and Venu Madhav Govindu. Robust relative rotation averaging. IEEE Transactions on Pattern Analysis and Machine Intelligence, 40(4):958-972, 2017.
[20] Chen Chen, Lei He, Hongsheng Li, and Junzhou Huang. Fast iteratively reweighted least squares algorithms for analysis-based sparse reconstruction. Medical Image Analysis, 49:141- 152, 2018.
[21] Chen Chen, Junzhou Huang, Lei He, and Hongsheng Li. Preconditioning for accelerated iteratively reweighted least squares in structured sparsity reconstruction. In Proceedings of the IEEE Conference on Computer Vision and Puttern Recognition (CVPR 2014), pages 2713-2720, 2014.
[22] Lisha Chen and Jianhua Z Huang. Sparse reduced-rank regression for simultaneous di- mension reduction and variable selection. Journal of the American Statistical Association, (77550):1533-1544 2012.
[23] Wei Chen. Simultancously sparse and low-rank matrix reconstruction via nonconvex and nonseparable regularization. IEEE Transactions on Signal Processing, 66(20)55333-5333, 2018.
[24] Yuxin Chen and Yuejie Chi. Robust spectral compressed sensing via structured matrix completion. IEEE Transactions on Information Theory, 60(10))6576-6611, 2014.
[25] Alexandre d' Aspremont, Laurent Ghaoui, Michael Jordan, and Gert R Lanckriet. A direct formulation for sparse PCA using semidefinite programming. In Advances in Neural Information Processing Systems (NIPS 2005), pages 41 41-48, 2005.
[26] I. Daubechies, R. Vere, M. Fornasier, and C.S. Güntürk. Iteratively reweighted least squares minimization for sparse recovery. Commun. Pure Appl. Marh., 63:1-38, 2010.
[27] Achiya Dax. On extremum properties of orthogonal quotients matrices. Linear Algebra and Applications, 433 2010.
[28] David L. Donoho. Compressed sensing. IEEE Transactions on Information Theory, 5284):1289-
1306, 2006.
[29] Henrik Eisenmann, Felix Krahmer, Max Pfeffer, and André Uschmajew. Riemannian thresh- olding methods for row-sparse and low-rank matrix recovery. Numerical Algorithums, pages 1-25, 2022.
[30] Michael Flad and Michal Aharon. Image denoising via sparse and redundant representations over learned dictionaries. IEEE Transactions on Image Processing, 15(12):336-33755, 2006.
[31] Maryam Fazel, Haitham Hindi, and Stephen 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.
[32] William Fedus, Barret Zoph, and Noam Shazeer. Switch transformers: Scaling to trillion parameter models with simple and efficient sparsity. The Journal of Machine Learning Research, (000))2222-277000 2022.
[33] Massimo Fornasier, Johannes Maly, and Valeriya Naumova. Robust recovery of low-rank matrices with not-orthogonal sparse decomposition from incomplete measurements. Applied Mathematics and Computation, 392, 2021.
[34] Massimo Fornasier, Steffen Peter, Holger Rauhut, and Stephan Worm. Conjugate gradient acceleration of iteratively re- reightet least squares methods. Comput. Optim. Appl., 65(1):205.
259, 2016.
12

---

[35] Massimo Fornasier, Holger Rauhut, and Rachel Ward. Low-rank matrix recovery via iteratively reweighted least squares minimization. SIAM Journal on Optimization, 2144);61141660, 2011.
[36] Simon Foucart. Hard thresholding pursuit: algorithm for compressive sensing. SIAM Journal on Numerical Analysis, 49(6):2533-22 25663,3333,3,,,, 2011.
[37] Simon Foucart, Rémi Gribonval, Laurent Jacques, and Holger Rauhut. Jointly low-rank and bisparse recovery: Questions and partial answers. Analysis and Applications, 18(01):25-48, 2020.
[38] Simon Foucart and Holger Rauhut. A Mathematical Introduction to Compressive Sensing. Birkhäuser Basel, 2013.
[39] Jonathan Frankle and Michael Carbin. The lottery ticket hypothesis: Finding sparse, trainable neural networks. In International Conference on Learning Representations, 2019.
[40] Paris V Giampouras, Konstantinos E Themelis, Athamasios A Rontogiannis, and Konstanti- nos D Koutroumbas. Simultaneously sparse and low-rank abundance matrix estimation for hyperspectral image unmixing. IEEE Transactions Geoscience and Remote Sensing, 54(8):4775-4789, 2016.
[41] Mohammad Golbabaee and Pierre Vandergheynst. Hyperspectral image compressed sensing via low-rank and joint-sparse matrix recovery. Proceedings of ICASSP 2012, 2012.
[42] Irina Gorodnitsky and Bhaskar Rao. Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm. IEEE Trans. Signal Process. 45(3):600-616, 1997.
[43] Benjamin D Haeffele and René Vidal. Structured low-rank matrix factorization: Global optimality, algorithms, and applications. IEEE Transactions Pattern Analysis and Machine Intelligence, 4466):1448-B8 1482, 2019.
[44] Torsten Hoefler, Dan Alistarh, Tal Ben-Nun, Nikoli Dryden, and Alexandra Peste. Sparsity in deep learning: Pruning and growth for efficient inference and training neural networks. The Journal of Machine Learning Research,
2021.
[45] Paul W. Holland and Roy E. Welsch. Robust regression using iteratively reweighted least- squares. Communications Statistics Theory and Methods, 6(9):813 827, 1977.
[46] Kishore Jaganathan, Samet Oymak, and Babak Hassibi. Sparse phase retrieval: Uniqueness guarantees and recovery algorithms. IEEE Transactions on Signal Processing, 65(9):2402- 2410, 2017.
[47] Gauri Jagatap and Chinmay Hegde. Fast, sample-efficient algorithms for structured phase retrieval. Advances in Neural Information Processing Systems (NeurIPS 2017), 2017.
[48] Prateek Jain, Praneeth Netrapalli, and Sujay Sanghavi. Low-rank matrix completion using alternating minimization. Proceedings of the forty-fifth Annual ACM Symposium on Theory of Computing, pages 665 674, 2013.
[49] Stuart M Jefferies and Julian C Christou. Restoration of astronomical images by iterative blind deconvolution. The Astrophysical Journal, 415:862, 1993.
[50] Chengyu Ke, Miju Ahn, Sunyoung Shin, and Yifei Lou. Iteratively reweighted group lasso based on log-composite regularization. SIAM Journal on Scientific Computing, 5365): S678, 2021.
[51] Michael V. Klibanov, Paul E. Sacks, and Alexander V. Tikhonravov. The phase retrieval problem. Inverse problems, 11(1):1, 1995.
[52] Martin Kliesch, Stanislaw J Szarek, and Peter Jung. Simultaneous structures in convex signal recovery- revisiting the convex combination of norms. Frontiers in Applied Mathematics and Statistics, 5:23, 2019.
13

---

[53] Christian Kümmerle. Understanding and Enhancing Data Recovery Algorithms: From Noise- Blind Sparse Recovery to Reweighted Methods for Row-Rank Matrix Optimization. Dissertation, Technical University of Munich, 2019. Available at https: bittps://nmtt //tp://tmmmmammmmmmmmmmmmmmm .b.
tum. wtt.dd//oo//214366666333 1121/122211111 pdf.
[54] Christian Kümmerle and Claudio Mayrink Verdun. A scalable second order method for ill-conditioned matrix completion from few samples. In International Conference on Machine Learning (ICML 2021), 2021.
[55] Christian Kümmerle, Claudio Mayrink Verdun, and Dominik Stöger. Iteratively reweighted least squares for basis pursuit with global linear convergence rate. Advances in Neural Information Processing Systems (NeurIPS 2021), 2021.
[56] Christian Kümmerle and Juliane Sigl. Harmonic mean iteratively reweighted least squares for low-rank matrix recovery. The Journal of Machine Learning Research, 19(1):1815-186; 2018.
[57] Ming-Jun Lai, Yangyang Xu, and Wotao Yin. Improved iteratively reweighted least squares for unconstrained smoothed Eg my-miiimization. SIAM Journal on Numerical Analysis, 51())9227 957, 2013.
[58] Kenneth Lange. MM optimization algorithms. SIAM, 2016.
[59] Kiryung Lee. Yanjun Li. Marius Junge, and Yoram Bresler. Blind recovery of sparse signals from subsampled convolution. IEEE Transactions on Information Theory, 63(2):802-821, 2016.
[60] Kiryung Lee, Yihong Wu, and Yoram Bresler. Near-optimal compressed sensing a class of sparse low-rank matrices via sparse power factorization. IEEE Transuctions Information Theory, 64(3):1666-1698 2018.
[61] Seong Hun Lee and Javier Civera. Hara: A hierarchical approach for robust rotation averaging. In Proceedings of the IEEE/CVF Conference on Computer Vision and Puttern Recognition (CVPR 2022), pages 15777-15786, 2022.
[62] Stamatios Lefkimmiatis and Iaroslav Sergeevich 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.
[63] Adrian S. Lewis and Hristo S. Sendov. Nonsmooth Analysis of Singular Values. Part I: Theory, Set- Valued Analysis, 2132221,, 2005.
[64] Xiaodong Li, Shuyang Ling, Thomas Strohmer, and Ke Wei. Rapid, robust, and reliable blind deconvolution via nonconvex optimization. Applied and Computational Harmonic Analysis, 47(3):893-9344 2019.
[65] He Lyu and Rongrong Wang. An exact sin 0 formula for matrix perturbation analysis and its applications. arXiv preprint arXiv:2011.07669, 2020.
[66] C Ma, K Wang, Y Chi, and Yuxin Chen. Implicit regularization in nonconvex statistical estimation: Gradient descent converges linearly for phase retrieval, matrix completion, and blind deconvolution. Foundations of Computational Mathematics, 2020.
[67] Malik Magdon-Ismail Np-hardness and inapproximability of sparse pea. Information Pro- cessing Letters, 126:35-35 30, 2017.
[68] Julien Mairal, Francis Bach, Jean Ponce, al. Sparse modeling for image and vision processing. Foundations and Trends@ in Computer Graphics and Vision, 8(2-3):85-283, 2014.
[69] Johannes Maly. Robust sensing of low-rank matrices with non-orthogonal sparse decomposi- tion. Applied and Computational Harmonic Analysis, 2023.
[70] Karthik Mohan and Maryam Fazel. Iterative reweighted algorithms for matrix rank minimiza- tion. J. Mach. Learn. Res., 13(1):44413443, 2012.
14

---

[71] Cun Mu, Bo Huang, John Wright, and Donald Goldfarb. Square deal: Lower bounds and improved relaxations for tensor recovery. In Proceedings of the 31st International Conference Machine Learning (ICMI. 2014), volume 32 of Proceedings of Machine Learning Research, pages 73-81, 2014.
[72] Bhaskar Mukhoty. Govind Gopakumar, Prateck Jain, and Purushottam Kar. Globally- convergent iteratively reweighted least squares for robust regression problems. In International Conference on Artificial Intelligence and Statistics, pages 313-322 2019.
[73] Peter Ochs, Alexey Dosovitskiy, Thomas Brox, and Thomas Pock. On iteratively reweighted algorithms for nonsmooth nonconvex optimization in computer vision. SIAM Journal on Imaging Sciences, 8(1):331-372, 2015.
[74] Samet Oymak, Amin Jalali, Maryam Fazel, Yonina C Eldar, and Babak Hassibi. Simultane- ously structured models with application to sparse and low-rank matrices. IEEE Transactions on Information Theory, 61(5):2886 2908, 2015.
[75] Liangzu Peng, Christian Kümmerle, and René 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.
[76] Liangzu Peng. Christian Kümmerle, and René Vidal. On the Convergence of IRLS and Its Variants nutlier-Robur Estimation. In IEEE Conference on Computer Vision and Pattern Recognition (CVPR 2023), 2023.
[77] Saiprasad Ravishankar, Jong Chul Ye, and Jeffrey A Fessler. Image reconstruction: From sparsity to data-adaptive methods and machine learning. Proceedings of the IEEE, 108(1):86- 109, 2019.
[78] Benjamin Recht, Maryam Fazel, and Pablo A. Parrilo. Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization. SIAM Review, 52(3):471-5011 2010.
[79] Emile Richard, Guillaume R Obozinski, and Jean- Dean-Philippe Vert. Tight convex relaxations for sparse matrix factorization. Advances in Neural Information Processing Systems (NIPS 2014), 27, 2014.
[80] Emile Richard, Pierre-André Savalle, and Nicolas 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.
[81] Masaki Samejima and Yasuyuki Matsushita. Fast general norm approximation via itera- tively reweighted least squares. In Computer Vision-ACCV 2016 Workshops: ACCV 2016 International Workshops, Taipei, Taiwan, November 20-24, 2016, pages 207-221. Springer, 2017.
[82] Steffen Schotthöfer, Emanuele Zangrando. Jonas Kusch, Gianluca Ceruti, and Francesco 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.
[83] Laixi Shi and Yuejie Chi. Manifold gradient descent solves multi-channel sparse blind deconvolution provably and efficiently. IEEE Transactions on Information Theory, 67(7):4784- 4811, 2021.
[84] Yunpeng Shi, Cole M Wyeth, and Gilad Lerman. Robust group synchronization via quadratic programming. In International Conference on Machine Learning (ICML 2022), pages 20095- 20105. PMLR. 2022.
[85] Mahdi Soltanolkotabi. Structured signal recovery from quadratic measurements: Breaking sample complexity barriers via nonconvex optimization. IEEE Transactions on Information Theory, 65(4):2374-2400, 2019.
15

---

[86] Ying Sun, Prabhu Babu, and Daniel P Palomar. inporizationmmimmmzzaaiii algorithms in signal processing, communications, and machine learning. IEEE Transactions on Signal Processing, 65(3):794-816, 2016.
[87] Yue Sun, Samet Oymak, and Maryam 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 1 Jun 2020.
[88] Andrei Nikolajevits Tikhonov. Solution of incorrectly formulated problems and the regulariza- tion method. Soviet Math., 1011003-1030 1038, 1963.
[89] Christos Tsisoss, Athanasios A Rontogiannis, and Kostas Berberidis. Distributed blind hyper- spectral unmixing via joint sparsity and low-rank constrained non-negative matrix factorization. IEEE Transactions on Computational Imaging, 3(2):160-174, 2017.
[90] Vandereycken. Low-Rank Matrix Completion by Riemannian Optimization. SIAM J. Optim., 23(2), 2013.
[91] Mukund Varma T, Xuxi Chen, Zhenyu Zhang, Tianlong Chen, Subhashini Venugopalan, and Zhangyang Wang. Sparse winning tickets are data-efficient image recognizers. Advances in Neural Information Processing Systems (NeurIPS 2022), 35:4652-444 26666, 2022.
[92] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin, Attention is all you need. Advances in Neural Information Processing Systems (NeurIPS 2017), 30, 2017.
[93] Curtis R Vogel. Computational methods for inverse problems. SIAM, 2002.
[94] Sinong Wang, Belinda Li,, Madian Khabsa, Han Fang, and Hao Ma. Linformer: Self-attention with linear complexity. arXiv preprini arXiv:2006. 04768, 2020.
[95] Endre Weiszfeld. Sur le point pour lequel somme des distances de points donnés est minimum. Tohoku Mathematical Journal, First Series, 43:355-386, 1937.
[96] Wei Wen, Cong Xu, Chunpeng Wu, Yandan Wang, Yiran Chen, and Hai Li. Coordinating filters for faster deep neural networks. In Proceedings of the IEEE International Conference on Computer Vision (CVPR 2017), pages 658-666, 2017.
[97] Max A. Woodbury. Inverting modified matrices. Memorandum report, 42(106):336, 1950.
[98] Quanming Yao, James Tin- Yau Kwok, and Bo Han. Efficient nonconvex regularized tensor completion with structure- proximal iterations. In International Conference on Machine Learning (ICML 2019), pages 7035-7044. PMLR, 2019.
[99] Jong Chul Ye, Yoseob Han, and Eunju Cha. Deep convolutional framelets: A general deep learning framework for inverse problems. SIAM Journal on Imaging Sciences, 11(2)9991-1048 2018.
[100] Ming Yu, Varun Gupta, and Mladen Kolar. Recovery of simultaneous low rank and two-way sparse coefficient matrices, a nonconvex approach. Electronic Journal of Statistics, 14:413-457, 2020.
[101] Zeinab Zeinalkhani and Amir H Banihashemi. Iterative reweighted Eg/2, recovery algorithms for compressed sensing of block sparse signals. IEEE Transactions on Signal Processing, (3(17):451-4553) 2015.
[102] Yanbin Zhang, Long-Tine Huang, Yangqing Li, Kai Zhang, and Changchuan Yin. Low-rank and sparse matrix recovery for hyperspectral image reconstruction using bayesian learning. Sensors, 22(1):343, 2022.
[103] Zhenyue Zhang, Zheng Zhai, and Limin Li. Graph refinement via simultaneously low-rank and sparse approximation. SIAM Journal on Scientific Computing,
2022.
[104] Hui Zou, Trevor Hastie, and Robert Tibshirani. Sparse principal component analysis. Journal of Computational and Graphical Statistics, 15(2):265-286. 2006.
16

---

Supplementary material for Recovering Simultaneuusly Structured Data via Non-Convex Iteratively Reweighted Least Squares
This supplement is structured as follows.
Appendix A presents some details about the experimental setup well additional numerical experiments.
Appendix B.1 presents the proof of Theorem 2.3.
Appendix B.2 presents the proof of Theorem 2.6.
Appendix details some technical results that are used in Appendices B. and B.2.
A Experimental Setup and Supplementary Experiments
In this section, we elaborate on the detailed experimental setup that was used in Section 4 of the main paper. Furthermore, we provide additional experiments comparing the behavior of the three methods studied in Section 4 for linear measurement operators A that are closer to operators that can be encountered in applications of simultaneous low-rank and troup-sparse recovery. Finally, we shed light on the evolution of the objective function (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 4 were conducted using MATLAB implementations of the three al- gorithms on different Linux machines using MATLAB versions R2019b or R2022b. total, the preparation and execution of the experiments used approximately 1200 CPU hours. The CPU models used in the simulations are Dual 18-Cort 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) [60], we used our custom imple- mentation of [60, Algorithm 4 thresholding (RiemAdaIHT) [29], we used an implementation provided to us by Max Pfeffer in private communications. We refer to Section 3 for implementation details for the IRLS method Algorithm 1.
and for Riemannian adaptive iterative hard
In all phase transition experiments, we define successful recovery such that the relative Frobenius emrr ( () -x
|x of the iterate (K) returned by the algorithm relative the simultaneously low-rank and row-sparse ground truth matrix X, is smaller than the threshold 10-1, As stopping criteria, we
[[x[(t -x(-1)
used the criterion that the relative change of Frobenius norm satisfies
tol for IRLS, the change in the matrix factors norms satisfy || - Uk-1 < tol and V V <
X
tol for SPF, and the norm of the Riemannian gradient in Ri emAdaIHT being smaller than tol for tol = 10 10-16 or if a maximal number of iterations is reached. This iteration threshold was chosen as max iter = 250 for IRLS and SPF and as max_iter iter- = 2000 for Ri RimAAaIHT, reflecting the fact that RiemAdaIHT is gradient-type method which might need many iterations to reach 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 X, eventually.
In the experiments, chose random ground truths X, E E. x of rank and row-sparsity s such that X, X./ X. | where X, U,diag (d.) V;, and where U. E g is matrix with non-zero rows whose location is chosen uniformly at random and whose entries are drawn from i.idd standard Gaussian random variables, d has i.idd standard Gaussian entries and V. g xr has likewise i.i.d. standard Gaussian entries.
A.2 Random Rank-One Measurements
In Section 4, we considered only measurement operator R X - R whose matrix repre- sentation consists of i.i.d. Gaussian entries, i.e., operators such that there are independent matrices
17

---

Ai.... Am with i..d. standard Gaussian entries such that
A(X), = X),
any Xn2, While known that such Gaussian measurement operators satisfy the (r, s)PP of Section 4, which is the basis of our convergence theorem Theorem 2.5, in regime of near-optimal number measurements with high probability, practically relevant measurement operators often more structured; another downside of dense Gaussian mcasurements 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 [46, 11, or blind deconvolution [59, 83], the measurement operator consists of rank-one measurements. For this reason, we now conduct experiments in settings related to the ones depicted in Figure and Figure 2 Section 4, but for random rank-one measurements where the action of A R2 Xniz + E on X can be written
A(X), = (a, ba,
for each j = 1, ..... m, where aj, by are independent random standard Gaussian vectors. In Figure we report the phase transition performance of RiemAda SPF and IRLS for (256 x x40) 2256x40)-diimnsiie dimeeninnla ground truths of different row-sparsities and different ranks if we are given such random rank-one measurements.
Number messonments
200
200
200
150
150
150
100
100
100
50
50
50
10
20
30
40
10
20
30
40
10
20
30
40
1
1,000
1,000
1,000
0.8
0.6
800
800
800
0.4
600
600
600
0.2
400
400
400
0
20
40
60
20
40
60
20
40
60
Row sparsity
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 X, (cf. Figure 1. Second row: Rank-ō ground truth X, (cf. Figure 2).
We observe in Figure 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 performs clearly best for both r = 1 and r = 5, whereas SPF has the second best performance for 1, it is somewhat unclear whether RiemAdaIHT or SPF performs better.
success and high empirical success rate is extends across larger area. IRLS
5. For
A.3 Discrete Fourier Rank-One Measurements
We now revisit the experiments of Appendix A.2 for third measurement setup motivated from blind deconvolution problems [3, 59, 64, 66, 83, 29], which prevalent in astronomy, medical imaging and communications engineering [49, 13]. particular, these settings, if z I is an (unknown) signal and _ is an (unknown) convolution kernel, assume convolution ÿ = * w. If we know that z Au for some known matrix A E Roxm and an s-sparse vector v E "2 applying the discrete Fourier transform (represented via the DFT matrix F E Cm we write the coordinates of
given the entries of their
T and w Bv for some known matrix B IR mmx and arbitrary vector
y Fÿ diag Fz Fw (laagFAA)FBy
18

---

as
In this setting, Theorem 2.6 guarantees that the behavior of the IRLS methods still benign as the y, A(=v"), (FA)), uv"rr
for each. 1,..., m, which allows to write the problem simultaneously rank-1 and s-row sparse recovery problem from Fourier-type measurements.
Number inssorcennnt
200
200
200
150
150
150
100
100
100
50
50
50
10
20
30
40
10
20
30
40
10
20
30
40
1
1,000
1,000
1,000
0.8
0.6
800
800
800
0.4
600
600
600
0.2
400
400
400
0
20
20
40
60
20
40
60
Row spursity
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 X, Second row: Rank-5 ground truth X, (cf. Figure 2).
In Figure
we report the results simulations with and 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 increasing number of measurement (with fixed dimension parameters 71 and n2) happens earlier than for the random Gaussian rank-one measurements of Appendix A.2, but slightly later than for dense Gaussian measurements. Again, IRLS exhibits the best empirical data-efficiency with sharpest phase transition curves.
As summary, we observe that IRLS is able to recovery simultaneously low-rank and row-sparse matrices empirically from fewer measurements than state-of-the methods for variety of linear measurement operators, including in cases where the RIP assumption of Definition 2.4 is not satisfied and in cases that are relevant for applications.
A.4 Evolution of Objective Values
While Theorem 2.5 guarantees local convergence if the measurement operator AA is generic enough and contains enough measurements (RIP-assumption), is instructive to study the behavior of Algorithm in situations where there are not enough measurements available to identify a specific low-rank and row-sparse ground truth X, which respect to which the measurements have been taken.
sequence of E- and 5-smoothed log-objectives (Fse. (x)) k21 from (4) is non-increasing. In Figure 6, we illustrate the evolution of the relative Frobenius error of an iterate to the ground truth X,, the (Ek, 6x) -smoothedd logarithmic surrogates Fe.. (X(k) as well as of the rank and row-sparsity parts Firss (x)) and Fsp (X (x)) of the objective, respectively, in two typical situations.
In particular,
the evolution of these four quantities in the setting of data of dimensionality created as in the other experiments, where a number of oversampling factor of 3.0 and 1.0. respectively) dense Gaussian measurements are provided to Algorithm 1.
71 = 128, T22 E {20, 40). } = 20 and 875 and 175 (corresponding
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 match the order of the (unsquared) relative Frobenius error),
19

---

converging to values between 10 10-13 13 and 10 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 (k) does not converge to X.. However, can be seen that Algorithm still finds very meaningful solutions: It can be seen that within 86 iterations, Fe..,h, (X()) converges
10-12 (since FF (X(A)) 10-6) in manner that partially "stairease-like": After 20 initial iterations where Fs (X (x) decreases significantly at each iteration, its decrease is dominated by relatively sudden, alternating drops of the (blue) sparsity objective Fop.b, (x)) and the (red) rank objective Fir.z (xk),, which typically do not occur simultaneously.
This illustrates the self-balancing property of the two objective terms in the IRLS objective FeL.b, (X(+)): while the final iterate at iteration 86 is not of the target row-sparsity 20 and r = 5, it is still 20-row sparse and has essentially rank 6. This means that Algorithm has found an alternative parsimonious solution 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 Ri RimAAdIHH, which typically breakdown in such a regime.
Rank Objective V Fir.ey (x())
Spariity Objective v Fap..6 (X())) Rei. Frob. error (x( X. |F/ X.
IRI i Chijective / Fee s, (x()))
102
10
10
10-2
10
10
10-6
10
10-5
10-10
10-6
10-14
10-7
10
12
0
20
40
60
80
iteration &
Figure Objective/ error quantities iterates (k) for iterations k. Left: Typical result for ni = 128.
40,
875. Right: Typical result for nı 128,
20,
175.
A.5 Robustness under Noisy Measurements
The convergence theory for the IRLS method Algorithm established in Theorem 2.5 assume that exact linear measurements y = AA(X.) of a row-sparse and low-rank ground truth X, are provided to the algorithm. However, in practice, one would expect that the linear measurement model is only approximately accurate. For IRLS for sparse vector recovery, theoretical guarantees have been established for this in [26, 57]. We do not extend such results to the simultaneously structured case, but we provide numerical evidence that IRLS as defined in Algorithm can be used directly also for noisy measurements.
To this end, we conduct an experiment the problem setup of Figure in Section for a fixed row-sparsity of 40, in which the measurements provided to the algorithms IRLS, RiemAdaIHT
20

---

and SPF are such that
y= A(X.) + w,
where w a Gaussian vector (i.i.d, entries) with standard deviation of V SNR is a varying signal-to-noks ratio. We consider SNRs between 10 and 10 and report the resulting relative Frobenius error statistics in Figure 7.
mn SNR and where
10
X XollF/ Koll//XXo XolF
IRLS- LowRankRow.Sprree
10
Sparse Power Factorization [60]
Riemannian Adaptive IHT [29]
10-2
10-3
10
10-5
relutvv Prohenius
10 10-6
10-7
und
10
102
103
104
105
10
107
10
10
10
10
10
Sigoal-to-noise
Figure 7: Median relative Frobenius reconstruction errors of different algorithms given noisy Gaussian measurements, n, = 256, n2 = 40, row-sparsinn sow-sparstty 40 and rank r 1, oversampling factor of 3. Error bars correspond 25% and 75% percentiles.
We observe that the reconstruction error is consistently roughly proportional to the inverse square root of the signal-to-noise ratio, for all three algorithms considered. This suggests that IRLS is as noise robust as comparable algorithms, and expected to return estimates of the original ground truth that has reconstruction error that is of the order of the norm of the noise.
B Proofs
The following two sections contain the proofs of our main results. Let begin with some helpful observations.
First note that the low-rank promoting part W can be re-written as
R -> B x of our weight operator
Wh X(0) ssk (Z)
U U1] (H(
Ex)
Uº [
zi
VH
V+
(18)
Vi
where
(k)((), (k) ,xx) = min (ex/o((t) (+)
min (Ek/OG(1) (k) 1)
iij=
rk
c (2) o
(3)
i,j=1
i,j=1
R
d1-Pk
1
1,j=1
Consequently, all weight operators in Definition 2.1 are self-adjoint and positive. Whereas for W ( ,,, this is obvious, for Wir xx.). it followw from the matrix representation
Wh (()) (IU U DHOw,e
IU
U o Vil')
21

---

where Di( .oy ] R" is diagonal matrix with the entries of H(o((),,))) Exx) which are all positive, on diagonal.
B.1 Proof of Theorem 2.5
Before approaching the proof of Theorem 2.5, let us collect various important observations. order to keep the presentation concise, we defer part of the proofs to Appendix C.
For rank-r matrix UEV* we define the tangent space of the manifold of rank-r matrices at Z as
Tu.v = {UZ; + Z2 Z1 E E Z2 E ER
(19)
In a similar manner, we can define for Z UEV* E M,,, and S = supp (z) supp(U) C [m]] the tangent space of M, restricted to S
Tu.v,s {UZ] Z₂ V* Z₁ E _ Ze E R with supp (Z)) S}.
(20)
As the following lemma shows, orthogonal projections onto the sets M No.. Tu.v. and Tu. can be efficiently computed.
Lemma B.1. We denote the projection operators onto Ma and N by T, and H₂- T truncates a matrix to the dominant singular values; H, sets all but the in E2-norm largest rows to zero. In case of ambiguities (multiple singular values/rows of same magnitude), by convention we choose the (respectively 8) with smallest index.
For U and V fixed, the orthogonal projection onto Tu.v is given by
Fu,v == Pru. Z UU + ZVV* ZU
For S C. [n1] and U, V fixed with supp(U)- = S, the orthogonal projection onto Tu, is given by Pu,v,s == Pru.v.s Z = (ps( ZVV* - UU"ZVV")
UU P3ZVV* UU"ZVV",
where Pss projects to the row support i.e., sets all rows to zero which are not indexed by
The proof Lemma B.1 is provided Appendix C.2. In contrast to the above named projections, the projection onto M is not tractable. However, [29, Lemma 2.4] shows that locally IP P,M,, can be replaced by the concatenation of T, and H., i.e., for Z, E &,,,. and Z., one has that
P P.M.. (Z) = (Z)..
For a matrix X E Ra X and i E [n,]. we set p; (X) = (X)x,ll2 where is row index corre- sponding to the i-th largest row of X in f2-norm. More precisely, if X is a decreasing rearrangement of X with rows ordered by magnitude ε2-norm, then p. (X) (I2 As the following lemma shows, the quantity Ps (X) determines local neighborhood of X on which H,, preserves the row-support.
Lemma B.2. Let X x be matrix with row-support S [n]] and [S] Then, for any Z E R2ixx2 with X o0.2 == max,e(nu] [X Zz.:122 I po (X) the matrix (Z) has row-support S.
Proof: Note that
1 2 (x(
max [Z)).. - (X)a.lla S (E[ns]
implies that any non-zero row of X corresponds to non-zero row of H. (Z) and hence yields the claim.
A first important observation is that if A has the (r, s)-IP, then the norm of kernel elements of .A is bounded in the following way. Lemma B.3. If A has the (r, RR-IIP with (0, and U Raixy V Reaxe with supp(UU = S.
|S| < 8, then
|All
2-
|
V
(1
8)
for all = ker(A))
22

---

The proof of Lemma B.3 is presented in Appendix C.3.
Remark B.4. If A is a Gaussian operator with standard deviation m. one has with high probability that || ~ 112222
m
We use of the following lemma to characterize the solution of the weighted least squares problem (11). Its proof analogous [54, Lemma B.7] and [26, Lemma 5.2].
Lemma B.5. Let A : IR x - Rh' and y E R Let W ,x,), R R be the weight operator (8) defined based on the information of X(k) E Ruixma, Then the solution of the weighted least squares step (11) of Algorithm I
(
arg mimmi...... W 52,55 (X)),, A(X)=y
(21)
unique and solves (21) and only if
A(X ((+))+)) = and
(W e0, (x E) for all E e ker A.
(22)
For any iterate (k) of Algorithm I, we furthermore abbreviate the tangent space (20) of the fixed rank-r manifold M, restricted to at (x()) by
Th To V.s.
(23)
where U IR and V E R are matrices with leading" r singular vectors of Hx(((()) as columns, and S E [mi] is the support set of the s rows of (4) with largest Ey-norm.
The following lemma the first crucial tool for showing local quadratic convergence of Algorithm 1. Lemma B.6. Let X, E M,,,, and let (() be the k-th iterate of Algorithm I with rank and sparsity parameters F = r and 3 = s, let ōg, Ek be such that sę and r from Definition 2.1 satisfy 2 s and Th 2 r. Assume that there exists constant
such that
E F S PT (=)
for all E ker(
(24)
where Ts = TÜN is as defined in (23) for matrices Ü E Raixr and V E R of leading r
left and right singular vectors of (x(()))) and S [mi] is the support set of (X()). Assume furthermore that
|X X.I S min
min {
(48 19c S
1
Jo,
(X.)
(25)
Then,
| xx, s 2 min { (x(1)(((())
Ps+ Pa+1
(X(()) 2
EN
5x
(1 (x.)l.
X. x.112)
where |IM], = E, |M,il11 denotes the row-sum norm of a matrix M, and W x(00 and ap are the weight operators (9) and (10) from Definition 2.1.
( 8a.
The proof of Lemma B.6 presented in Appendix C.4.
Remark B.7. By revisiting the proof of Lemma B.6 (omit the bound in (48) and keep the term [3, WE) until the end), one can show under the same assumptions as in Lemma 6 that
(x(x)(((k)) Ek
p--1(X())) -
}
S 2 min
X EN P (P..
Ps. WP X
Ps- Ps:)=)
8%
AN
where U and V are containing the left and right singular vectors of ((), see Definition 2.
As U and V might not be unique, any of leading singular vectors be chosen in this definition.
23

---

The contribution of the norms of the weighted X. terms in Lemma B.6 can be controlled by Lemmas B.8 and B.9 below.
Lemma B.8. Let W X Eo
R Rnixna be the rank-based weight operator (9) that
the spectral information of ( and let X. E IRaixna be rank-r matrix. Assume that there exists Ç < such that
max |X X.I> S (X,)
(26)
Then for each 5500
12/9
1
Wi X(x, (X.)
$ -11- (X,,
kk + Ek K (X -x,1 xx 2 -x,11
11-5
5,
and
Wile X(0),s (x+) xx,01ss,
///
=k + x -x. -x.1s
s. (Ex + + (k) -x(1)xx11)
$ (1- (X+) (-
where Kg is such that Kg 22/2221///// for S S 2 and 4 S Kg = for 2 < q S 4 and K₂ = for q OC.
Lemma B.9. Let W X 5y E R be the row-sparsity-basee weight operator (10) that uses the current iterate (1) with dx min (8x-1, Pe+1 ( (x)) and let X. E & x be s-row-sparse
matrix, Assume that there exists 0 < < < 2 such that
X,11,,22 = max ( (b)
|| (xx())
-X,,))llll S (X,)
(27)
sE|ml
where (M( denotes the fn-norm of the in ly-norm s-largest row M. Then
s6
X00,8 X. S (1 -1- (1=)))p((X.)
Lemma B.8 a refined version of [54, Lemma B.9] the proof of which we omit here. The proof of Lemma B.9 is provided in Appendix C.5. Finally, the following lemma will allow us to control the decay of the IRLS parameters dū and Ek-
Lemma B.10 ([54, Lemma B.5]). Let X, E Mrrs, assume that A has the (r, R)IP) with 5 E (0, 1), and let us abbreviate n = min(ns n1, n2)-
Assume that the k-th iterate (k) of Algorithm I with r and updates the smoothing parameters in (12) such that one of the statements and that rk 2 T and 85 2 Furthermore, let
(X(1) or 8r. (x(()) is true,
E3 < 1
(X.)
48
|A 12 lel (k) with 91A22-12 V 1 (1-6)
X X, satisfy
[=() S min { (X.). min
148' 21c1 jo,(
1
(28)
Then
E ( 222 2/2/ncc/All
V 46 + 82.
The proof of Lemma B.10 10 is provided in Appendix C.6. We finally have all the tools prove Theorem 2.5. Note that (14) implies
min (x.) (x())) r s
x (1) X.I < min
}
(29)
48c2 Alla-2
4 5nici
'This result is technical result of an unpublished paper. In this paper, we only use that result a tool. If the reviewers think that adding the proof is relevant here, are happy to provide it.
24

---

and
1
Ek S
48 (x..)
(30)
which we will use in the proof below. The latter follows from the fact that for i =
min (Ex-1, & 0r+1 (k)) S X (k)
x. S (x,(x)448.
Proof Theorem 2.5: First note, that by assumption r and s. Furthermore, since 48c All2-22
S 1
the closeness assumption (29) implies that H2 (x()) and X. share the same support due Lemma B.2. Let X (K) be the k-th iterate Algorithm 1. Since the operator A: Rmixug R has the (r, s)-P with
5 E (0, 1), Lemma B.3 yields for all E " . E R with supp(U) S. IS S 8, that
r 1 91412-2 Plo.v.s
for any E E ker(.A), Furthermore, due to our assumption that and r, the smoothing parameter update rules in (12), i.c., 6g = min ( Pot1 (x(+)) and Ex min (6 imply that 2 and
su 2 s for all k. We can thus apply Lemma B.6 for =(h) == x X. (note at this point that (29) implies the closeness assumption (25) of Lemma B.6) and obtain
c (((+1) - X..
S Acfi 12
min (o-1(((())
(
Wecss. (x.)l. +
| x
X.
Ex
82
(31)
where With x((),
Xm2 + R X the low-rank promoting part (9) of the weight operator associated to E BRiXX the sparsity promoting part (10), Since by assumption
(k) and W
max(E1,
[=(+ S
48 o. (X.),,
Lemma B.8 yields
[ x (x.)
S 0.995
42
cer 25 | x.1I 211x - x.113).
(32)
(X.))
Similarly, by assumption
10.2 S
(K)
(X.)
18 15 (X.),
488
such that Lemma B.9 yields
21
X, 0.995 20p.(X.)
(33)
000.8 8
Inserting (32) and (33) into (31) we obtain that
0r (x(1
Pati (x(1)
(k+1)
0.995
min
(34)
xx(X) 22 2=x
(+) + 21)=
Pe (x.
Due to the assertion that 2 r, it holds that En S 5o-1 (X (.)). Therefore, Lemma C.3 yiclds that
(o) 2ca
| (4) ( 5 133
and, since 8 2 8, also that
g < kk)
(4) 12
since S Ps+1 (x) in this
Thus, using the assertion that one of the statements EN = a (4) or 8n = (x(1( is true, we obtain from (34) that
(++1)
0.995
4.20f1Al1a->
(oxxx.) + o.2x.) ||
(35)
25

---

For (
min( =x(x.),
ex(X.)] (as implicd by (29)), this yields
18 1All2-2
0
(36)
the quadratic error decay
(++1)
1
if we define = 4.179ci.All2-> (x()
2x
To show the remaining statement, need argue that the assertions of Theorem 2.5 satisfied not only for but any + 8 with 2 For this, is sufficient show that
1. 2
2. 8k+1 2 8,
3. Eb+1 S 18 (X.))
4. (29) holds for ((++1), and that
5. one of the statements Ei-1 (x(1(X((+1)) ox+1 Pa+1 Po+1 (X(++1) is true,
as in this case, ((++)) 2 X. follows by induction due to successive application of (36).
For 1. and 2., we see that this follows from the smoothing parameter update rules (12) which imply that Ex+1 S 1 and 6 I P (X (++1)).
follows from (30) and the fact that due (12), (E#) is non-increasing. is satistied due (36) and (29). To show 5., we note that due to (29), the assertion (28) is satisfied, and therefore it follows from (35) and Lemma B.10 that
(K++))
< 4.179cfiAllsss
(x,(X))
5r
2s
|
2 / 45 + o2.
Pa (x(
We now distinguish the case (i) dx < Eu and the case (ii) de 2 Ek.
case (i), holds that
o (x (x(++)) $
[ S 4.179c/iAA122. 2-2
5r (x.).)
2v 10nc All2-12 10
p.(X.)
= p2v 10me
ER.
where the last inequality holds since by (29) the k-th iterate X additionally satisfies
x.I |l
(37)
2pv Luv/DeeA22-2
In this case, due to the smoothing parameter update rule (12), we have that Ek+1 ar+t (1))))). In case (ii), we have likewise that
Pa+1 (x (++1) < = (B(ttli/| S. n2v 10nc 1012m [
|
8n < ou.
due to (35), Lemma B. 10, and (37). Henee, 5et1 Po+1 (X((+1) which shows the remaining statement 5. and concludes the proof of Thcorem 2.5.
B.2 Proof of Theorem 2.6
1.) Ect > 0 be arbitrary. Due to the additive structure of FF, ef. (4), is sufficient to establish
that
F 5p.5 (2) S (Z(x) = Fxp.s (X) + (VF sp.5 (X)), X) z- W X)) (38)
for any Z. E R Xn2 where Wh X.5
g22 R mnix. xn=
defined analogously to (10) and
Fine (Z) < Qir.s (Z(X) = Fine + (VFir X), Z x z
W X)), (39)
for any Z.X E 1 where W
: 121xn + R np. is defined analogously to (9).
26

---

The argument for (38) is standard in the IRLS literature [2, 73, 76] and is based on the facts that both 2sp. (Z) and (Z) are row-wise separable, and that t + fvT(vt) is concave and therefore majorized by its linearization: indeed, let 9r IR be such that
log(e)t//T). if [t] -r(1) |t)
if it S.
The function 9- (-) is continuously differentiable with derivative g1) = ? sign(t) and furthermore, concave restricted to the non-negative domain Rzo- Therefore, it holds for any 1, R2z0 that
(+(1) S 97 (t())))))))))))) (t)(
We recall the definition (t) T log(et (T) for 7 max x2) Thus, for any 2, I., it follows that
and (t) 2/2²² for |t] S from (2)
with derivative (t( max((2 3r
f+(2) = =()))=)((²))) S g72 ( + gra (z (=² x²)
= (x +
2 max -2) (2" x²))
and inserting 7 7= -, = = Zz. 12, x | and summing over = 1,... my implies that
71
8
F (Z)
Sa(1 Z Zi. (2) S F w. (X)
|X,,: 1
2 max X 12, 82)
Z X 1|
= F X)
1 Z x,,) 2
max X,, 11, 82)
max 1,,11//00, 2//2,11
i=1
From the chain rule, it follows that for all i = 1,..., n1 for which X,: * 0,
d) X1:112
8² X
X1,
82X,
16 /6(1)x,, f6(1I,,,
dx,
max X, 82)
12,
max X
12 52)
(40)
and therefore
(Z) S Fxp.((X) + (VF sp,6 (X) 1-
X (Z -2-
X, W X.5
(Z X))
which shows the majorization of (10).
(38), recalling the definition
W X,5
diag [max (II 02 1)
The majorization of (39) is non-trivial but follows in straightforward way from [53, Theorem 2.4] as the objective Fire (Z) corresponds to the one of [53, Theorem 2.4] up to multiplicative factor of e² and constant additive factors, and since the weight operator W operator used in [53, Chapter 2] for 0.
X.s corresponds to the weight
2.) Due to the definition (3) of Fup.6, (-) and the derivative computation of (40), we observe that d,
VF sp (x()) = diag ( (((((())) 12/6//
(a)
(b)
W
(k)
x Sa.
i=1
comparing the resulting term with the definition of (10) of Wap "(1))6," equality follows from the the formula
Furthermore, an analogue
((050) max (1o(8))/((,11))
d
V- V
Firce (X (()) = [U U1] diag
i=1
with oi (x))) for any i < d, which is a direct consequence from the calculus of spectral functions Lemma C. and inserting into the low-rank promoting weight operator formula (9)
Wele xe (x (1)) 1x))= = U U1 Eo diag (o
(k)
d
Σー1
V
V (x))
V
i=1
27

---

Inserting V (x(1) W
x( (x) and V.rre (X(4) = W
X( (x(*)
into the defini-
tions of Dire: (2|X (k)) and 2v..d. (2)x(1)), we see that it holds that
Cires ( = FFriee (x(1))
(Z, W x(h,e (2)) - (A), Wx X(4) (X(1)))
+ 2
and
2x (x) Fop.6 (x(1)
1
(Z, Wap
Z
[X
(k)
W
X x
2
X .6.
Therefore, we see that the weighted least squares solution X of (11) for k + 1 coincides with the minimizer of
Qinic (x) 2xp.6 (x))
min Z:A(Z)
min 2:A(Z)=y
Firca (X(1)) + F 5p (x(())
(41)
1 (Z, (x()) 52. (Z)) (X (K), W x 5e (x
2
with the weight operator W ,x, of (8), which implies that
Cirer (X ( (x)) + 6. (X (1+1) X (k)) S Qir.es (x(*) (1) + Qup.bu (x(*) (x())) (42)
Using the majorization (16) established in Statement of Theorem 2.6 and (42), follows that
Fen,ba (x S Ores (X(++1) (()) + 2ap (x (x(++)xx(())
I Qinor (x(1) (x + 2xp,8e (x (x()) (8))
(43)
= Firex (x()) + Furde (x()) = Fex (x
using in the third line that Qreee (x(x) = Firen (X(1)) and Qxp. (X())
=
Fsp.b. (x).
To conclude, it suffices to show that e + Fe,hy (X(++1) and 8 + Fer (X(++1) are non-decreasing functions, since (43) then extends to
Je- (X (1)) S Fea (x (x(+))) S Fen.b (x(++)) S FF (x()),
where we used that the sequences EY and dg defined in Algorithm I are decreasing, So let us prove this last claim. We define for é IR the function h RR50 + IR such that he(T) = fr (t),
if 2 |t|
(T(T)) T log(et if < [t].
This function continuously differentiable with (T(T) = for all r > |t and h((T) = (log(et? ///)//) -
for r < [t], which implies that (T(r) 2 0 for all 7 2 and thus shows that & +> Fe,d, (X(++)) and
8 + Fe,, (x (x(++)) are non-decreasing functions due to the additive structure of Fe, (X (k+tl) and
(4).
3.) First, we argue that (x (x))) )221 is a bounded sequence: Indeed, if == tom Ek > and
5 : lim 5% > 0, we note that
2 log(el|X 12/// +
1 2
log(e max ||X (k) 08)
1 = log(eof X (1) //2)
1
log(emax X (k) 112/0²)
2
2 & log(eoi (4) /
1 2
+ 02
1
log(e max ||x (k:) 112/02)
2
$ Finee (x) + Fyp.i (x = F (x( S Fe (x(1))
1 2 min(di, do)of(X())) (1)) =dh di max [IX Cxa,,
28

---

which implies that {|| (8) [1211 is bounded by a constant that depends on Cxo.
Furthermore, we note that the optimality condition of (11) (see Lemma B.5) implies that X (k+1)
satisfies
( (Wx,)) 16 (x =) = 0 for all = E ker A and A(X (k+1)) = y.
Choosing E = ((++1) - (*) and using the notation W = W Wx())e 61,66 we see that
X (k+))) (()) x (k- () (x) (x())
(x((tt)) w(()) (x ( (x(1) w(() (x - 2(W((k) ((( ())+ x(++1) x(()) ((x(+t2), w (X (k+1))) 2(W(()) (x()),, x(++2)) + (x(*), w(( ((()(
(44)
= -(((((+t) w (k) (x (k+ )))
Due the definition of (k),, we note that its smallest singular value (interpreted as matrix operator) can be lower bounded by
miin (W(k))= min (W X(0,) + o min ( x 6 max | +
x()) + + /
where X()) and (()()) are constants that satisfy X(.) S exp(Cxux/0 - and (1,((1) S E² (xx 1).
Combining this with (44), the monotonicity according Statement of Theorem 2.6, and (41), follows that
Fen F Fereribu (x(++))
2 [x (8) (++1)),, (x())(((()) (k++))))
2
882 /0spxx1)) + /(1x).))) || (+(1)
13.
2
Summing over all this implies that limg [X+++] - xx() ||= =0. =.
Since (x(k) is bounded, each subsequence of (k)) has convergent subsequence. Let (X(*)) 122 be such sequence with lime+ (ks) X, i.e., X is accumulation point of the ( there exists weight
sequence. As the weight operator ((k) depends continuously operator W Xma - RB1xan such that W lime- ((k)),
Since lim |X (x+1) x (1:) || = 0, it also holds that X (k++1) + X and therefore
VVF (X) E) (W (X). E
lim (W(ke) (x(k))(X((t)) ())+ E)
E-300
for all E ker A. The statement is shown as this is equivalent to X being stationary point of (.)(()) subject to the linear constraint (Z E 1 x A(Z) yy.
C Technical addendum
C.1 Auxiliary Results
In the proof of Theorem 2.6, we use the following result about the calculus of spectral functions. Lemma C.1 ([63],35, Proposition 7.4]). Let F: IRd xdz + R. be spectral function F' = f f a with an associated function f: d IR that is absolutely permutation symmetric. Then, F is differentiable at X Rdh xda and only if f is differentiable at oXX) Idd.
In this case, the gradient F of F at X is given by
VF(X) = Udiag (V)((())
if X (o V* for unitary matrices Rd xd: and V 6
'Here, for v E gmin(d;.dy diag(v) E gd xd2 refers to the matrix with diagonal elements v, on its main diagonal and zeros elsewhere.
29

---

C.2 Proof of Lemma B.1
The projection operators for M, " N/ and Tu.v are well-known, see e.g. [8]. To see the final statement assume that U has row-stpport S and note that Pu,v.s is idempotent, i.e.
Pu.v.sPu.v.s2
-UUZZ - P3ZVV' +
-U" PsZVV UU'ZVV"
UU UUMP,ZVV' UU + UU + VV
UU - -U' UU'ZVV'
= UU + Pg UU' = Fu,v
One can easily check that Pu. acts as identity when applied to matrices in Vu.v and that Fu,v U.v,s since (Z', s = (Pu,v.sZ' Z) F. for any Z, Z' This proves the claim.
C.3 Proof of Lemma B.3
Let E E ker (A). Note that
0 ||A(E)) 1| A(P ru.v.s [ P Tu. {
||A(( (=))ll2
||
|| (E)
By the RIP we hence get that
Pru.s (B)]l -
1 [((Prrvv. (
1
A(T
Tu.v. [3)
1 5
1 1
- 5)
||Al|
[P] Tu.sss
3 |
Consequently,
[ | | Pruvvs
[|
+
P p
Vu.
2 $
1+ (1+
|| ||A1
(1 (1-6)
|| | [
C.4 Proof of Lemma B.6
In the proof of Lemma B.6 we will use the following fact. Lemma C.2. Let Wir xxt,een be the weight operator defined in (9), which is based on the matrices
U E 1R and V E R of leading To left and right singular vectors of (k), IfM E TuN, then Wi (M) E Tu.v. IfM E Tu.y. then Wi x00,
(M) E TU.
X00) .EN
Proof: If E Tu.v. there exist Mi E R Mr E 1 rk), M3 E (g( such that
Mi
M.
V
M U
U.]
M
V v
0
e.g., see 190, Proposition 2.1]. We thus observe that the weight operator W
B Xng
R x from
(18) satisfics
W
Xw., .ek
(M)
U
U1]
H(( 2x
(E)))
U+
U
U.]
M
M2
V Vi
V
V* vul))
U
M
0
Vi
M
M2
V
UU
U1]
H(o
Ma
Vi
H
(k)
M
H (k)
M2
= U
U.]
E Tu.v
H
M3
Vi
Similarly, if M E T..y there exists M, E (m (ma mw)
such that M = U. ULM,(V_)" and
Wh
X(1)
(M)
U
U
H
M. Vi
V
U
U
Tt...
M
30

---

Proof of Lemma B.6: Let = 12 First note that by Lemma B.I, T2 Ts
x be arbitrary. We start with some simple but technical observations.
Fr E (Id Pr,,) (E)
= (Ü (Ü -U
(45)
(Id -o,v) (ld Ps\EEV'
P v = EVV"
with
(F.. (Ū - V V.VV =0,
(46)
where used that U. P - U for U1 E - and V1
R (na-(( being the complementary orthonormal bases of U and V, and that ViV=0.
Second, let now U E R and V R= be matrices with leading left and right singular vectors of (
their columns which coincide with the matrices and V from Definition 2.1 in their first columns. Then follows from (45) and (46) that
(47)
= -| + + Pú vt. I (VV
(V--vv
$2 < | Pt + IVs VVV
By an argument analogous to (46), we observe that
P |I 12 P P Ps (B,E),
where E (M) is an element of the subspace M Mi M2 R21 that is the direct sum of the subspaces Mi = [Pt,v : Z E R=1*"3] and M2 == ZVV Z E R"1*=2]
Let now 4 (e),,e, be the rank promoting part the weight operator from (9) and W x(0)5, be the row-sparsity promoting part from (10). Note that restriction of
W = Pú W ,,, LU. + Ps W
(()), Ps=
M invertible as first summand is invertible on M₁ Tú, second summand invertible on M22 (recall that the weight operators are positive definite). and M1 1 .M2. Therefore it holds that
|b.v= 1 |IPs l = - = (WME,, W Wind
Sai (wiz)) (a, W.ME) omin (Wim)
W
E., (E.WE) wE
Omin (Pb.y xe P(.))
(PeSo xi Pss.)
+ Tmii
IM
IM
E WE)
P7 (X
In the first inequality, we used that WIM is positive definite. In the second inequality, we used that Gmin (A +
B) > (A) + aann (B), for any positive semidefinite operators A and B, and and in the third incquality that
[ WE) < (=, WE). The latter ohservation can he deduced as follows: Note that, by the self-adjointness of
W.
+ (PM( (P) (E, WE) = PP ) WP (B) (PM))
E w E + (E),
WP)
31

---

Since W is positive semi-definite (due to the fact that both W W *x00_ .co
and W x are positive definite). all
that remains is to argue that the mixed term on the right-hand side vanishes. To this end, note that PU v Z Ps-ZVV = 0, for any Z E Mr and compute
( (P)((E)) = PL) e (P.. (3)I)) + (PM Ps W
X Ps PM.
= (Pú.v(PL())) W X (PU.v(PM, + (Fs- W X00, P)==E)
= 0.
We now continue by estimating
(=, w=) = (s.Pb. x( Pú.vE) + (=, Ps W x Ps-
(Ft.v(E), W a Pú. Pú.v(E)) PP X P30
(48)
I E, W =) (=, (=. w (1, =)
=(E, (E.W
using the positive semidefiniteness of Wit x(,) and W in the last inequality. To be precise, the last inequality can be argued as follows: Due to complimentary supports and 5", we see that
x(K),5 &
(E, W
E
S W
Ps
(Ps E, W
x Ps
x(k) 82
(
+ sP W
Ps. + SP W
x PoE)
(49)
S. ( P..
Similarly, note that W ((()_ cc acts diagonally on Tu,v and Tú.v. Indeed, we have by Lemma C.2 that if
M Tu.v, then W W (0) (M) E Tu,v and Tú,y, then W x(t)) (M) E Tú,y, which implies
(Fu, W x()),, (PU.vE) = (Pu.v,, W , (Pu. =
due to the orthogonality of elements in TŪ v and Tu.v, respectively, and therefore it follows from = TE) (E) that
( W (E)) (Pu. M (Pu, (Pt. W , (Fo,vE))
+ (Pu.v,, W (Pú,v)) + (PU. W (Pu.vE))
(50)
2 (Plyv W ( (P).
Combining the previous estimates with (47) and noticing that
pi (x(1)) (x(())
min
2
0%
( 1 (X(TT
we obtain
[/ S min (x(())
(x(1)
|
E. Wy (x)))) ,ct)
E
o2
.cc
(51)
+ 2
+2| (Po. - ) +2 EV
Next, control the last summands in (51) using matrix perturbation results. Recall that U, R"1 and V. E x are the singular vector matrices the reduced singular value decomposition of X.. First observe that
Pivv - Pu..v.)
|00 U.U, - -0 v.V
=
jūū u.u: |
[Id v.v.ll +
Id عي
< -ūū0~u.u V - -v.v:) |
32

---

Now note that, by [12. Lemma 1] and [65, Theorem 3.5], we obtain
= - u.u: -vv: s (IU.011+ +I..+IIIVVII VIll 54 ox(),) 0r+i (X(()))))
(X x.
| X. < 0, (X.) 0x+1(((6)
I (1 - |x((] - X.I
1/48)),(X..)
(52)
where we used some small observations in the last two incqualities: First, 0++1 O (H.(X()))) I o (X )
which follows from the rectangular Cauchy interlacing theorem [27, Theorem 23]. Second, according to
Lemma B.2 and (25), the row-support S of Ha 4x( (*) coincides with the row-support S. = {i E [nn]
Ixxx)):: o 0} of X. and hence H. (x
( X. is submatrix of (() X.. Finally, (x
[T ()) x [|) s IX. - x [1 S 48 (X.) due to (25). Consequently,
|IPUN - Po..v.
S 1x
x. (1-1/48)o, (X.)
E
and, by similar argument,
[X(() X.
I(Pu..v. - Pl.v) $ (1 - 1/18)0,, (X.)
E|
such that it follows from U Id PU.v and P. U.v Id Pü,v that
| - Pú
S -
(Pü.
Pu..v. |= + V
=
|X X-I|
I 8 (1 1/18)o..X.) (X.)
|ا||
(53)
To estimate the fourth summand in (51), we argue analogously that
| -
< E(VV vvr)) - (IVV - v. + IIV. V. vvi) |
22 V.. || II.. |5||
(x) X.|
[ x, ||
22 (x.)) 0r+1 (H, [X (x))
|
(x,) 0r+1 (X()))
(x X.
(+) X.
< 4 or (X.) - or+1 1 (x (
| < 4 (1 148) (x.)
| [
(54)
using again that 0r+1 (x(1)) S X - x| S 15 ,x, due (25).
Let now E ((++1)-
x (k+1) X.. Combining (24) and (51)-(52) can proceed to estimate that
| 12
S o² | (
2 F
(
(x (=tt))), W o.. ((t1))
$ 26² min
Ek
de
2
17718) )²+²²²
8
2 +(1+7)
(17/18)
4
1 (x X.I1
xX().
$ 26² min
{
(x (A))
(X()) (1)
(=(++)) W x (=(+))))) 167c? 1 (++1)) 12
(196)²
(x()
(x))(
< 26 min
X c. ((+t))))
1 1
E#
(55)
where the third inequality follows from (8) and (25). Hence, rearranging (55) yields
1 < 46d min { (x()(((())
(x
}
2
E
[ (k+1) W Wx()) 6, (E(kt2))).
Ek
6n
33

---

By Lemma B.5, we know that X (k+1) fulfills
0 = (= (k+1) W 56. (x ((+)))) = -=((+) W Wxim (E(t)t))) + (=(t11), W Wxin (X.))
such that we conclude that
/ = 1 1|
0r+1(X(()) (
p=+1(X(6)) 2
1c2 min
(= g(kt)) W x 2e [= 1)
EK
(r+1((((h) Pati(X(())
=4 min
(= (=(++) W (k) (X.))
de
(x(1(((h))
(x)))
2
4o min
(=(t)) W Wim (x.)) (=(+t1)
w x x.))
8k
4c min
(x(1)
P (et1(((()) (x(1))
W
x 184
(x.)
= {(+
+
|
why
X, xx1112)
x(x) 50
Ek
46 min
(x(1)((()))
Dat1 (x(1))
( (x.)
W
x.
(1..2)
|
Ek
dx
x
which completes the proof. We used the penultimate line Hölder's inequality and that
F) = Σ
An, Bas
$ Σ
Aa,lll [B...ll2 I ( maxx [|A.:ll)) Σ B..lll &
|All
|
Ad
for all matrices A, B.
C.5 Proof of Lemma B.9
Note that by Lemma B.2 and (27), S = supp(H,(X(*))) supp(X.). Since by assumption
5% S (x()) we have by definition of W X( that Z = W X(1)6 X. is matrix with row-support S and rows
Z min
(X. (X.)
(X.)..
(X(())
(k) )x. 11
(1(x(1))
(4)) )1.1
for i E S. Now note that if (27) holds, then
(2.1( ( (x..))p(Xx))
52
I2..112
<
where we used the last estimate that with (27) and || )t.112 2 (X.)). for S, we have
I|xx (A) Ja. 12 N
(1(x(()) (X.) (x.))..1122
N (a1 (11- (1-))
for all i E S. The claim easily follows since Z has only s non-zero rows.
C.6 Proof of Lemma B.10
In the proof of Lemma B.10, we use simple technical observation. Lemma C.3. Let X, E Mr.s, let ((k) be the k-th iterate of Algorithm 1, and abbreviate (k)
X (k) X.. Then the following two statements hold true:
1. 1f S S0++1 (x(1), then En < [E()
2. If 0₁ S Pot1 (x (k)), then 8r. 10.2
34

---

Proof: By defining [x(+), to be the best rank-r approximation of X in any unitarily invariant norm, we bound
S ()++1(((*))
[Xx
[ [-11 - I|x(
x.
==
where the incquality follows the fact that X, is rank-r matrix. Similarly, for the second statement, we have that
S Ps+1 ( [|x( (x(x())) 1oo 51 1x X. using that X, is s-TOw sparse.
18 100.21
Proof of Lemma B.10: First. note that, using Lemma B.3, the observation in Remark B.7 yields
(r+1(X(()) (
| f I 4cjjAlla-2
min
EK
de
(56)
= (P W x(+ se Pú.v
Ps W
P
x 6
for all R as the assumption (28) implies [X
X.II S min
(x.(
Thus, this holds
48
Toc
particular also for ( (*)
X.. (Recall that and V contain the leading singular vectors of (k)
see Definition 2.1.) We estimate that
( Wit x(t))ee ( E(t) = || -| x ( ("UV
(1 x ) (PO.
+ ( )) (1,. (x)))
m2n(n1
$ <|((W (22()))
max S o0 + vn rek
o
Furthermore, since max(e:, 1I use variant of Lemma B.8 obtain
35 (X.) ) by assumption and Ek
[1 by Lemma C.3,
[ x x.)
48 47
(VEEL + 2EL-
|
$
ox(X.)
(X.)
< 1004/TE + 1.16 (k) (X))
On the other hand, note that ( restricted to S coincides with the restriction of (k) to S under assumption (28), cf. Lemma B.2, and therefore
( X(1 Ps-=(M) (())))))))))))))))
W
X (4) W
(x(()) )
(N) ),, 1 17 S (n1 - 8)0
max
With the estimate of above, this implies that
E, (Pt. W
Pú.v + Ps. W
x(8),6v Pso) )=
X00
< ( x(
1/2 (x
+ /n - TEX +(11 8)02
- rek +
13
52 11
| (k) } 3(nz - -ree + (1 - -388
(X.)
Inserting these estimates into (56), we obtain
< 4 min
r.+1 (1))
(x (+)
13
1 nEE + 21
52 (4) 112
)
|
Es
8c
of (
If now either one of the two equations
1001111XX (k))
85 Ps+1 (x(1) (8) is true, it follows that
11 S AiAl (13ne + 4n1
+ ciAll-z
208
(2)
(X.)
< < (13ne + 4n
(k)
if the proximity condition || (x()) X. S 2ll or (X.) is satisfied. Rearranging the latter inequality yields the conclusion of Lemma B.10.
35


Reviews:
[{"review_id": "C63I17wpw0", "reviewer": "Reviewer_bSFz", "summary": "This work studies the problem of recovering a low-rank and row-sparse matrix from its compressed linear observations. A method based on iteratively re-weighted least squares is proposed, in which the sparsity inducing function is a non-convex log function. For theoretical contributions, this work provides a local quadratic convergence analysis of the iterates of the proposed algorithm under a restricted isometry property assumption, and it also studies the convergence of function value for any linear measurement operator using a Majorize-Minimize algorithm framework. Numerical experiments on synthetic datasets show that the proposed method can recover the ground truth with fewer measurements than some state-of-the-art methods.", "questions": "(1) In Theorem 2.4, will the assumptions $\\varepsilon_k = \\sigma_{r+1} (X^{(k)})$, $\\delta_k = \\rho_{s+1}(X^{k})$, and $\\varepsilon_k \\leq \\sigma_r(X_\\star) / 48$ ever be satisfied? \n\n(2) The inequality at the end of page 6 only holds for certain $k$ satisfying the assumptions, so how do we arrive at the conclusion $X^{(k+l)} \\to X_\\star$ from such an inequality for one step?\n\n(3) In Theorem 2.5, how strong/restricted is the assumption that $\\delta_k$ and $\\varepsilon_k$ have limits, and the limits are positive?", "limitations": "As far as I see, there is no potential negative societal impact of this work.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "1. This work provides convergence analyses not only for linear operators satisfying restricted isometry property but also for more general linear operators. Furthermore, it mentions when the restricted isometry property can be satisfied, and it also discusses the radius of the neighborhood in the local quadratic convergence.\n\n2. The proposed algorithm is well explained by identifying that the underlying family of objectives that are minimized during the iterations is (2).\n\n3. There are experiments for cases where the rank and sparsity are not accurately known as a prior, and thus the sensitivity to hyper-parameter choice is discussed, and also it provides a better simulation of potential real-world scenarios.", "weaknesses": "(1) The algorithm requires estimates of the rank and row sparsity, and the theoretical conclusions require that these estimates are accurate, i.e., $\\tilde r = r$ and $\\tilde s = s$. These priors and assumptions can limit the application of this method and its theoretical convergence guarantees.\n\n(2) There are some assumptions on the iterates $\\varepsilon_k$ and $\\delta_k$ in Theorem 2.4 and Theorem 2.5, but it is not clear whether these assumptions will ever be satisfied. Please refer to question (1) and (3) for details.\n\n(3) In numerical experiments, it would be better to report the time cost of different methods, including both the per iteration time cost and the total time cost.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your constructive and very detailed feedback to our submission. A point-by-point\nresponse and clarification with regards to your comments follows below:\n \n> (1) The algorithm requires estimates of the rank and row sparsity, and the theoretical conclusions require that these estimates are accurate, i.e., $\\widetilde{r} = r$ and $\\widetilde{s}=s$. These priors and assumptions can limit the application of this method and its theoretical convergence guarantees.\n\nWe agree that our theoretical convergence analysis is limited to this choice of hyperparameters. We, however, would like to point out that competing methods we are aware of, including the ones of papers [24] and [53], likewise need rank and row sparsity estimates. Furthermore, the experiments we performed in Fig.\\ 2, in which the rank parameter was overestimated by $100\\%$ and the sparsity parameter was overestimated by $50\\%$, suggest that IRLS is more robust to poor tuning of the hyperparameters than other methods such as RiemAdaIHT and SPF of papers [24] and [53]. Finally, we would like to point out that the IRLS method we propose does not require any additional hyperparameters except rank and row sparsity estimates such as a step size choice, unlike, e.g., RiemAdaIHT.\n\n> (3) In numerical experiments, it would be better to report the time cost of different methods, including both the per iteration time cost and the total time cost.\n\nWe agree that this would be insightful. We decided against including direct timing comparisons with other methods as a fair such comparison would require efficient implementations of all algorithms involved; since we were more interested in the _quality_ of the returned solutions, we focussed less on efficiency. We refer to our general answer above for a discussion of the per iteration cost.\n\n> (2) There are some assumptions on the iterates $\\varepsilon_k$ and $\\delta_k$ in Theorem 2.4 and Theorem 2.5, but it is not clear whether these assumptions will ever be satisfied. Please refer to question (1) and (3) for details.\n\n**Questions**:\n\n> (1) In Theorem 2.4, will the assumptions $\\varepsilon_k = \\sigma_{r+1}(\\mathbf{X}^{(k)})$, $\\delta_k=\\rho_{s+1}(\\mathbf{X}^{(k)})$, and $\\varepsilon_k \\leq \\sigma_r(\\mathbf{X}_*)$ ever be satisfied?\n\nWe note that it suffices if $\\varepsilon_k = \\sigma_{k+1}(\\mathbf{X}^{(k)})$ \\emph{or} $\\delta_k = \\rho_{s+1}(\\mathbf{X}^{(k)})$ at some step $k$, i.e., we do  not require that both of these statements are satisfied at the same time. \n\n> (2) The inequality at the end of page 6 only holds for certain $k$ satisfying the assumptions, so how do we arrive at the conclusion $\\mathbf{X}^{(k+\\ell)} \\to \\mathbf{X}_*$ from such an inequality for one step?\n\nNote that this condition is always fulfilled in the first step of Algorithm 1 since both parameters are initialized with $\\infty$. Whether the condition still holds in later steps depends on the quality of initialization. The proof of Theorem 2.4 works via induction and relies on initializing sufficiently close to the ground truth. Given such an initialization, the proof shows that $\\varepsilon_k = \\sigma_{k+1}(\\mathbf{X}^{(k)})$ or $\\delta_k = \\rho_{s+1}(\\mathbf{X}^{(k)})$ holds for all consequent steps.\n\nRegarding the condition $\\varepsilon_k \\leq \\sigma_{r}(X_{\\star}) /48 $, we noticed that this condition is actually superfluous since the locality assumption on $\\mathbf{X}^{(k)}$ in (15) automatically implies that $\\varepsilon_k \\leq \\sigma_r(\\mathbf{X}_*)/48$. Just note that \n    \n $\\varepsilon_k = \\min\\left(\\varepsilon_{k-1}, \\sigma_{r+1}(\\mathbf{X}^{(k)})\\right) \\leq \\sigma_{r+1}(\\mathbf{X}^{(k)})$\n        $\\leq || \\mathbf{X}^{(k)} - X_*|| \\leq \\sigma_r(X_*)/48$\n    if $\\widetilde{r} = r$ as in the assumptions of Theorem 2.4. In the second inequality we used that $\\mathbf{X}_*$ is of rank $r$ and in the last inequality we used (15). We will remove (14) in the revised version. Thanks for catching this! \n\n> (3) In Theorem 2.5, how strong/restricted is the assumption that $\\delta_k$ and $\\varepsilon_k$ have limits, and the limits are positive?\n        \nRegarding Question 3, note that the sequences $\\delta_k$ and $\\varepsilon_k$ are, by definition, decreasing and bounded from below, i.e., by the monotone convergence theorem they converge to a limit in any case.\n\n\nPlease let us know if you have further questions or concerns regarding our submission. If we\ncould clarify your questions, we would appreciate it very much if you considered an adjustment of\nyour rating. Thank you!"}, {"review_id": "Ak6ErKLJf4", "reviewer": "Reviewer_W2en", "summary": "Paper proposes to solve inverse problems on matrices subject to multiple types of sparsity (e.g. low-rank and element-wise sparsity) using an algorithm designed for the non-convex objective functions involved. The algorithm is a re-weighted least squares, which despite not solving a convex problem, authors prove a local convergence and a global consistency theorem for. They illustrate the performance of their algorithm on synthetic data. ", "questions": "1. what is the sensitivity of the algorithm to performing approximate SVDs\n2. Figure 3: fewer iterates are needed, OK, but what's the FLOP in each iterate? what's the overall computational cost of achieving epsilon close results as compared to others? is anything parallelizable?\n3. A related stream of research (see Tight Convex Relaxations for Sparse Matrix Factorization, NIPS 2014) characterizes the atoms that form the basis in which a doubly sparse matrix is sparse. How does this work, especially the lower bound Omega(r(s1+s2)), compare to those methods which directly minimize for the desired structure to recover?", "limitations": "because the computational complexity is not deeply discussed, it's hard to know how applicable the methods are to larger data", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The algorithm proposed solves a difficult non-convex problem, and the theoretical results by authors are original and interesting. ", "weaknesses": "1. The paper contains statistical results of the type \"if at least m measurements are available then the algorithm achieves good performance\", but is missing a more thorough discussion on algorithmic complexity. In particular, I would be curious to know if the truncated SVD in the Update Smoothing step is the heaviest computational piece of the algorithm or if the computational bottleneck is somewhere else. \n2. Numerical experiments are run on very small data. To really show the benefit of two sparsity types I would believe that bigger values of n1 n2 result in more convincing evidence that there is benefit in mixing the two prior knowledges\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your constructive and very detailed feedback to our submission. A point-by-point\nresponse to your comments follows below:\n\n> The paper contains statistical results of the type ”if at least m measurements are available\nthen the algorithm achieves good performance”, but is missing a more thorough discussion\non algorithmic complexity. In particular, I would be curious to know if the truncated SVD\nin the Update Smoothing step is the heaviest computational piece of the algorithm or if the\ncomputational bottleneck is somewhere else.\n\n> Numerical experiments are run on very small data. To really show the benefit of two sparsity types I would believe that bigger values of $n_1$ and $n_2$ result in more convincing evidence that there is benefit in mixing the two prior knowledges.\n\nOn the one hand, we agree with the reviewer that a larger scale of experiments is always preferable. On the other hand, we would like to emphasize that the theoretical results on locally quadratic convergence rates and monotonic decrease of the objective function are guaranteed to hold in higher dimensions as well. In evaluating these, the scale of the specific simulation does not matter. Finally, note that the scales of the problem instances we use are comparable to the ones used in the papers [24] and [53] describing the state-of-the-art\nmethods RiemAdaIHT and SPF, respectively.\n\n**Questions**:\n\n> (1) What is the sensitivity of the algorithm to performing approximate SVDs?\n\nThis is an interesting point. In fact, we did not test the present IRLS method for robustness against using approximate SVDs. However, in [47], Kümmerle et al. used a related IRLS method for low-rank recovery (i.e., a setting with only a single parsimonuous structure) in which case high-accuracy low-rank solutions were obtained using a truncated SVD implemented by a randomized block Krylov method [NIPS 2015, \"Randomized Block Krylov Methods for Stronger and Faster Approximate Singular Value Decomposition\", Musco et al.]. We therefore expect that this behavior extends to the simultaneously structured case considered here.\n\n> (2) Figure 3: fewer iterates are needed, OK, but what's the FLOP in each iterate? What's the overall computational cost of achieving epsilon close results as compared to others? Is anything parallelizable?\n\nWe refer to our general response above for a discussion of the per-iteration computational cost.\nRegarding a potential way to parallelize the main computational steps, it would be possible to use a conjugate gradient method to solve the main linear system of the IRLS step whose complexity is dominated by matrix-vector and matrix-matrix multiplications (similar to previous work on low-rank IRLS of [47]), which parallelize well. Second, the matrix multiplications needed to obtained the truncated SVD via block Krylov could also be parallelized.\n\n> (3) A related stream of research (see Tight Convex Relaxations for Sparse Matrix Factorization, NIPS 2014) characterizes the atoms that form the basis in which a doubly sparse matrix is sparse. How does this work, especially the lower bound Omega(r(s1+s2)), compare to those methods which directly minimize for the desired structure to recover?\n\nThanks for pointing us to this interesting work! We will be happy to include a reference to the NIPS 2014 paper by Richard et al. to the literature review in the final version. From a practical point of view, however, the paper does not propose a polynomial time algorithm for the problem, and the proposed heuristic algorithm focuses exclusively on Sparse PCA, which is a different problem setting than the underdetermined matrix recovery we consider. \n\nTherefore, it is hard to compare the numerical performance to our work right away. Regarding the theory, the estimates on the statistical dimension of the proposed (k, q)-trace norm provided by Richard et al. are restricted to (atomic) rank-1 matrices, whereas our local convergence theory for IRLS hold for arbitrary rank.\n\nPlease let us know if you have further questions or concerns regarding our submission. If we could clarify your questions, we would appreciate it very much if you considered an adjustment of your rating. Thank you!"}, {"review_id": "t2KQOqkX6r", "reviewer": "Reviewer_3C1j", "summary": "This article introduces an algorithm for recovering jointly row-sparse and low-rank matrices from (underdetermined) linear measurements. The algorithm is based on iteratively reweighted least squares for a non-convex objective. The method is theoretically analysed, establishing local quadratic convergence rates under the restricted isometry property. The second theoretical result relates the algorithm more directly to IRLS to show that the objective is non-increasing and proves convergence to a stationary point.", "questions": "- can you mention the per-iteration complexity of your method and contrast it with existing methods?\n- the introduction mentions simultaneously low rank and both column and row sparse. Can your results be extended to this setting?", "limitations": "Yes.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "I found the paper well-written, the problem of simultaneously recovering sparse and low rank matrices is challenging one cannot simply rely on standard convex regularisation. This work provides a nice solution to this problem via IRLS and the theoretical results are interesting: the authors establish both 'compressed sensing' type results that guarantee recovery under near optimal sampling complexity, as well as results that provide an understanding of the optimisation properties of the algorithm.", "weaknesses": "- The numerical results are somewhat limited, it would have been nice to have a discussion on the applications of the proposed method and more realistic numerical examples.\n- the result of Theorem 2.4 is restricted to the setting of exact measurements and does not cover robustness results, i.e. what if y has been corrupted with noise? Also, what if $X_*$ is only approximately sparse and low rank? Both of these are the more realistic settings, so it would have been nice to see a more complete result. Also, theory aside, equation (1) mentions additive noise, but it is not clear from the algorithm and results that the proposed method can handle additive noise.\n- Given that the title mentions \"simultaneously structured data\", I was expecting more than just low rank + sparse. It would be interesting to see the method extended to other kind of structures such as sparsity with respect to some dictionary. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your constructive and very detailed feedback to our submission. A point-by-point\nresponse to your comments follows below:\n\n> The numerical results are somewhat limited, it would have been nice to have a discussion on the applications of the proposed method and more realistic numerical examples.\n\nWe agree that additional experiments on some of the real-world applications mentioned in our introduction (like sparse phase retrieval, blind deconvolution, or hyperspectral imaging) could enrich the paper. Since the focus of this work was on algorithm development and analysis, and an efficient implementation of IRLS will depend on the specific problem instance (i.e., type of measurement operator), we refrained from heading into this direction to keep the exposition concise.\n\n> The result of Theorem 2.4 is restricted to the setting of exact measurements and does not cover robustness results, i.e. what if $\\mathbf{y}$ has been corrupted with noise? Also, what if $\\mathbf{X}_*$ is only approximately sparse and low rank? Both of these are the more realistic settings, so it would have been nice to see a more complete result. Also, theory aside, equation (1) mentions additive noise, but it is not clear from the algorithm and results that the proposed method can handle additive noise.\n\nCertainly, the noise robustness of the proposed methods if of interest both theoretically and empirically. We restricted ourselves to setup of exact measurements in the theory as even for the case of only low-rank structures, IRLS has less developed theory in this case. However, from an empirical point of view, IRLS (without modification of the optimization problem) returns a solution whose reconstruction error is of the order of a constant times the magnitude of the additive perturbation. We add a plot exploring this behavior in the rank-$1$ case of Figure 1 of the paper in the attached PDF above.\n\n> Given that the title mentions \"simultaneously structured data\", I was expecting more than just low rank + sparse. It would be interesting to see the method extended to other kind of structures such as sparsity with respect to some dictionary.\n\nThis is a good point. Indeed, one could add a fixed dictionary to the problem formulation which can be absorbed into the measurement operator in the analysis. The only point where we would expect some additional work to be necessary would be in proving that the concatenation of measurement operator and fixed dictionary still satisfy the RIP. For Gaussian $\\mathcal{A}$ this should be possible. Also note that the paper [55, Lefkimmiatis et al.], which we mentioned in our related work discussion and which restricts itself to single structures, combines IRLS with dictionary learning techniques, i.e., it addresses the question of how to learn a dictionary on the fly.\n\n> The introduction mentions simultaneously low rank and both column and row sparse. Can your results be extended to this setting?\n \n Thanks for mentioning this important aspect of our setting! While we restricted the presentation to the one-sided sparsity case for simplicity, it is indeed possible to adapt our IRLS method to handle sparsity on both sides of the matrix, i.e., column- and row-sparsity, by adding a second sparsity weighting term to the combined weight operator $W_{\\mathbf{X}^{(k)},\\varepsilon_{k},\\delta_{k}}$. \n   \nFor instance, in the symmetric setting $\\mathbf X_\\star = \\mathbf X_\\star^T$ (naturally occurring in sparse phase retrieval) one would define the weight operator $W_{X^{(k)},\\varepsilon_{k},\\delta_{k}}$ as in (8), but with an additional term that multiplies $W_{X^{(k)},\\delta_{k}}^{sp}$ _from the right_ to $Z$, which corresponds to minimizing the sum of _three_ smoothed logarithmic surrogates.\n\n"}, {"review_id": "YCIA36iiw9", "reviewer": "Reviewer_P5td", "summary": "This paper proposes an IRLS method for recovering data with multiple, heterogeneous low-dimensional structures from linear observations. It combines non-convex surrogates for row-sparsity and rank, to identify simultaneously row-sparse and low-rank matrices from limited measurements. Theoretical results are provided that show locally quadratic convergence. The experiments demonstrate favorable empirical convergence and prove its efficacy in handling challenging data recovery scenarios.", "questions": "It would be good to remind the reader what \\cal{F}_ functions are in (13).\n\nAre there other matrices besides Gaussian that are known to satisfy the RIP given in Definition 2.3?", "limitations": "The authors discuss the non-optimality of the bounds in Theorem 2.4, arguing that experimentally the convergence radius appears to be much larger. The authors also discuss future work about generalizing to other combined structures, which the reviewer agrees is interesting future work but certainly not a limitation to the results in this paper.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper is well-written overall and the contributions are clear. The challenge of the combination of structures is made apparent. The related work discussion is thorough. The IRLS is a popular approach and is practical even for large-scale systems. Experiments suggest the method is also robust to the choice of parameters r and s.", "weaknesses": "Although it is also true of other results in this area, only local convergence is guaranteed and practically it may be challenging to guarantee an initialization within the required radius. The results would be strengthened if other matrices besides Gaussian were shown to satisfy this RIP, and if others were even just used empirically and shown to still offer convergence.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your constructive and very detailed feedback to our submission. A point-by-point response to your comments follows below:\n\n> Although it is also true of other results in this area, only local convergence is guaranteed and practically it may be challenging to guarantee an initialization within the required radius. The results would be strengthened if other matrices besides Gaussian were shown to satisfy this RIP, and if others were even just used empirically and shown to still offer convergence.\n\nFirst, we would like to mention that our supplementary material includes experiments with rank-1 and subsampled Fourier measurements, which show similar outcomes as in the Gaussian case. Extending the theory to these measurement operators would be attractive. As long as the RIP holds for such measurements our theory directly applies. Extending the RIP, e.g., from Gaussian to subgaussian measurements is straight-forward. \n\nFurthermore, in [52, ``Blind recovery of sparse signals from subsampled convolution''], Lee et al. used that, for rank one matrices, subsampled Fourier measurements of a certain shape satisfy the RIP we use (requiring additional log-factors in the sample complexity). Since extending such results to more general settings requires substantial additional effort, this is however beyond the scope of our present work.\n\n**Questions**:\n> It would be good to remind the reader what $\\mathcal{F}$ functions are in (13).\n\nThanks for this good suggestion. In the final version of the manuscript, we will be happy to remind the reader that Equation (13), which defines the quadratic models minimized in the weighted least squares steps, uses as functions $\\mathcal{F}$ the $\\varepsilon_k$-smoothed log-determinant function and the $\\delta_k$-smoothed sum of logarithmic row-wise $\\ell_2$-norm, which had been previously defined in Equation (2)-(3).\n\n> Are there other matrices besides Gaussian that are known to satisfy the RIP given in Definition 2.3?\n\nApart from the straight-forward generalization to subgaussian measurement matrices and the subsampled Fourier measurements by [52, Lee et al.] (where the RIP is only needed for rank one matrices), we are not aware of further matrix types which provably satisfy our RIP with near-optimal rates.\n\nPlease let us know if you have further questions or concerns regarding our submission. If we could clarify your questions, we would appreciate it very much if you considered an adjustment of your rating. Thank you!"}]

Author rebuttal:

Meta-review: The paper proposed an iteratively reweighted least squares (IRLS) algorithm to recover data matrices that are simultaneously row-sparse and low-rank, which achieves a locally quadratic convergence under minimal sample complexity, outperforming methods that use a combination of convex surrogates. The paper received unanimous support from the reviewers.
Decision: Accept (poster)