Efficient Alternating Minimization with Applications to Weighted Low Rank Approximation
Authors: Zhao Song
Mingquan Ye
Junze Yin
Lichen Zhang
Venue: ICLR 2025 Poster (2025)
Keywords: weighted low rank approximation
sketching
alternating minimization
Primary area: optimization
TL;DR: We provide a fast alternating minimization algorithm for weighted low rank approximation that runs in nearly-verification time.

Abstract: Weighted low rank approximation is a fundamental problem in numerical linear algebra, and it has many applications in machine learning. Given a matrix $M \in \mathbb{R}^{n \times n}$, a non-negative weight matrix $W \in \mathbb{R}_{\geq 0}^{n \times n}$, a parameter $k$, the goal is to output two matrices $X,Y\in \mathbb{R}^{n \times k}$ such that $\\| W \circ (M - X Y^\top) \\|_F$ is minimized, where $\circ$ denotes the Hadamard product. It naturally generalizes the well-studied low rank matrix completion problem. Such a problem is known to be NP-hard and even hard to approximate assuming the Exponential Time Hypothesis. Meanwhile, alternating minimization is a good heuristic solution for weighted low rank approximation. In particular, [Li, Liang and Risteski, ICML'16] shows that, under mild assumptions, alternating minimization does provide provable guarantees. In this work, we develop an efficient and robust framework for alternating minimization that allows the alternating updates to be computed approximately. For weighted low rank approximation, this improves the runtime of [Li, Liang and Risteski, ICML'16] from $\\|W\\|_0k^2$ to $\\|W\\|_0 k$ where $\\|W\\|_0$ denotes the number of nonzero entries of the weight matrix. At the heart of our framework is a high-accuracy multiple response regression solver together with a robust analysis of alternating minimization.

Published as a conference paper at ICLR 2025
EFFICIENT ALTERNATING MINIMIZATION WITH APPLI- CATIONS TO WEIGHTED Low RANK APPROXIMATION
Zhao Song
Mingquan Ye
Simons Institute for the Theory of Computing, UC Berkeley magic. 1 inuukde@@gaa com
University of Illinois Chicago mye9Quic.edu
Junze Yin
Lichen Zhang MIT CSAIL 1i lichenzecsail wwt.::essaai..........
Rice University jy158erice.edu
ABSTRACT
Weighted low rank approximation is fundamental problem in numerical linear algebra, and it has many applications in machine learning. Given matrix M I="x" a non-negative weight matrix W E 2200 a parameter k, the goal is to output matrices X, Y E R such that ||W (M N0m-xY--- >l is minimized, where denotes the Hadamard product. It naturally generalizes the well-studied low rank matrix completion problem. Such problem is known to be NP-hard and even hard to approximate assuming the Exponential Time Hypothesis (Gillis & Glineur, 2011; I; Razenshteyn et al., 2016). Meanwhile, alternating minimization is a good heuristic solution for weighted low rank approximation. In particular, Li al. (2016) shows that, under mild assumptions, alternating minimization does provide provable guarantees. In this work, we develop an efficient and robust framework for alternating minimization that allows the alternating updates be computed approximately. For weighted low rank approximation, this improves the runtime of et al. (2016) from W llak2 to ||W]lok where |[W denotes the number of nonzero entries of the weight matrix. At the heart of our framework is a high-accuracy multiple response regression solver together with a robust analysis of alternating minimization.
INTRODUCTION
Given matrix M E _2xa the low rank approximation problem with rank it asks to find pair of matrices X, y E Ra** such that | -XYTI || is minimized over all rank k matrices X and Y, where || F is the Frobenius norm of matrix. Finding a low rank approximation efficiently is algorithmic problem that is well studied in machine learning, numerical linear algebra, and theoretical computer science. The exact solution follows directly from singular value decomposition (SVD): let M = UEVT and set X = Vg. Y V i.e., picking the space spanned by the top-k singular values and corresponding singular vectors. Faster algorithms utilizing linear sketches can run in input sparsity time (Clarkson & Woodruff, 2013). In addition to the standard model and Frobenius norm, low rank approximation has also been investigated in distributed setting (Boutsidis et al., 2016), for entrywise & norm (Song et al., 2017) and for tensors (Song et al., 2019c).
In practice, it is often the case that some entries of M are more important than others and some entries can be completely ignored, so it's natural to look for weighted low rank approximation. More specifically, given target matrix M E xn and a non-negative weight matrix W E 20 the goal
to find X, Y E gu with || (M XYT)|| minimized, where is the Hadamard product of two matrices. The formulation of weighted low rank approximation covers many interesting matrix problems, for example, the classic low rank approximation can be recovered by setting W 1, TT and the matrix completion problem (Jain al., 2013) is by observing a subset of entries of M, equivalent to picking W as a Boolean matrix. In addition to its theoretical importance, weighted low rank approximation also has significant practical impact in many fields, such natural language

---

Published as a conference paper at ICLR 2025
processing (Pennington et al., 2014; Arora et al., 2016; Hsu al., 2022), collaborative filtering (Srebro & Jaakkola, 2003; Koren et al., 2009; Lee et al., 2013; Chen et al., 2015), ecology (Robin et al., 2019; Kidzinski et al., 2022), chromatin conformation reconstruction Tuzhilina et al. (2022) and statistics (Wentzell et al., 1997; Markovsky et al., 2006).
Algorithmic study for weighted low rank approximation dates back to Young (1941). On the com- putational hardness front, Gillis & Glineur (2011) has shown that the general weighted low rank approximation is NP-hard even if the ground truth matrix is rank 1. The hardness is further en- hanced by Razenshteyn et al. (2016) by showing that assuming the Random Exponential Time Hypothesis. the problem hard to approximate beyond a constant factor. Despite its hardness, many heuristic approaches have been proposed and witnessed many successes. For example, Shpak (1990) implements gradient-based algorithms, while Lu et al. (1997); Lu & Antoniou (2003) use the alternating minimization framework. Srebro & Jaakkola (2003) develops algorithm based on expectation-maxinne
(EM). Unfortunately, all these approaches are without provable guaran- tees. Razenshteyn et al. (2016) is the first to provide algorithms with theoretical guarantees. They propose algorithms with parameterized complexity on different parameters of W, such as the number of distinct columns or low rank. In general, these algorithms are not polynomial which is also indicated by their lower bound results. Ban et al. (2019) subsequently studies the weighted low rank approximation problem with regularization, and they manage to obtain an improved running time depending on the statistical dimension of the input, rather than the rank. When one relaxes to a bi-criteria solution with additive error guarantees, Bhaskara et al. (2021) provides a greedy algorithm. Whenever all entries of the weight matrix are nonzero, Dai shows that it is possible to convert the additive error to multiplicative (Dai, 2023).
How to bypass the barrier of Razenshteyn et al. (2016) while still getting provable guarantees? Li et al. (2016) draws inspirations from matrix completion literature and views the problem as a low rank matrix recovery problem: suppose the matrix M E R'** is a noisy, full-rank observation that can be decomposed into M M + N where M the rank-k ground truth and N is the rank-(n - k) noise matrix. They then analyze the performance of alternating minimization when 1). the ground truth is incoherent, 2). weight matrix has a spectral gap to all-I's matrix, and 3). weight matrix is non-degenerate. Under these assumptions, they show that the alternating minimization algorithm provably finds a pair of matrices X, Y Baxk such that || XYXXXTI I| S O(k) ||W o N| +e, where || is the spectral norm of a matrix. This provides a solid theoretical ground on why alternating minimization works for weighted low rank approximation.
While the Li et al. (2016) analysis provides a polynomial time algorithm for weighted low rank approximation under certain assumptions, the algorithm itself is still far from efficient. In particular, the alternating minimization framework requires one to solve (n) different linear regressions exactly per iteration. The overall runtime of their algorithm is O([[]W]l k22 + nk log(1/e)) where [Wllo llo denote the number of nonzero entries in W, making it inefficient for practical deployment. Moreover, their analysis is non-robust, meaning that it cannot account for any error at each step. This is in drastic contrast with practice, where floating point errors and inexact solvers are used everywhere. In fact, there are good reasons for them to mandate exact regression solvers, as their algorithm only requires log(1/e iterations to converge and any fast but approximate regression solver might break the nice convergence behavior of the algorithm. Hence, we ask the following question:
Is il possible to obtain faster and more robust alternating minimization- algorithm with a similar convergence rate?
In this paper, we provide positive answer to this question. Specifically, we show that the alternating updates can be computed in nearly linear time each iteration and polynomially large errors can be tolerated. Both of these results rely fast, randomized and high-accuracy regression solver that uses sketching to compute a preconditioner. We summarize our main result in the following theorem:
Theorem 1.1 (Informal version of Theorem 4.6). There is an algorithm (see Algorithm 1) that runs in O( k + n log(L/e)) time and outputs a rank-k matrix M such that
|| M M S O(KT) ||W NI +
where T is the condition number of M and Ö(-) suppresses polylogriiii factors in n and k.
2

---

Published as a conference paper at ICLR 2025
Algorithm 1 Main Algorithm. The CLIP procedure zeros out rows whose & norm are large, and the QR procedure computes the QR decomposition of the matrix and outputs the orthonormal factor Q-
WENTTRRWEEGGIILLWRRRKKNN R=x", W Rax", e,k) 2: T + O(log(1/e)
1: procedure
Theorem 1.1
3: dak 1/ poly(n,
4:
> T is an estimate of the condition number of M.. Yo RANDOMINIT(H, k) Initialize Yo random Rademacher variables, scaled by vx
< 1/ (oly(n,
5:
6:
for f = 1 to T
7:
X, t- PASTMULLPLLREERRSIIONNNAA Y,-1, W, Eake Oak) Solve ((n) regressions
using sparsity of W and Algorithm 2.
8:
X, f- CLIP(X,)
Clip rows with large & norms.
9:
X, + QR(X,)
Y, < FASTMOUTTHPLRREGGESSSO(((((T X,, W Cski 6sk)
11:
V, < CLIP(Y,)
12:
Y₁ < QR(Y))
13:
end for
return M + &TYTL1
14:
15: end procedure
Remark 1.2. The general structure of our main algorithm (Algorithm 1) is based on the traditional alternating minimization method described in Li et al. (2016). We replace the exact update with an approximate update (lines and 10) based Algorithm 2, which makes the overall algorithm both faster and more robust. The remainder of the paper is dedicated to presenting a theoretical guarantee for its efficiency and robustness.
Roadmap. Section 2, we introduce several basic notations and definitions which we will use throughout this paper. In Section 3, we give brief overview of our techniques. In Section 4, present our main result. In Section 5, we give conclusion for this paper.
PRELIMINARY
In Section 2.1, introduce the basic notation used in this paper. In Section 2.2, present the background of the sketching technique, including the SRHT matrix and oblivious subspace embedding. In Section 2.3, we present the mathematical background and assumptions related to the weighted low rank approximation problem.
2.1 NOTATION
Let n, be arbitrary positive integers. We define a set [n] as (1,2,-..,,) We use R, RT R% and [ax to denote the sets containing all the real numbers, m-dimensionnll vectors with real entries, m-dimensional vectors with non-negative real entries, and n x m matrices with real entries.
Let E I and w T Let E [m). Let I;; E R represent the i-th entry of v. We use vx E R
to represent a vector satisfying (V²), = /zz. We define || = (Ei 1/2
Let A, W be two arbitrary matrices in gx Let i [n] and j [m]]. We use A,, E to represent column vector that is equal to the i-th row of A and A E R" represent a column vector the i-th row
that equal to the j-th column of A. Ar, E represents a entry of A, located and j-th column. diag(z) E Rux* represents the matrix satisfying diag(x);, = Z; if i j and diag(:)., = 0 if i + muz(A) represents the number of nonzero entries of A.
Suppose that n 2 We denote the spectral norm of A as |||] = suPzeRm || 12222/1/// /|////2 denote the Frobenius norm of A as || which is equal to (E EØ= 2/2, and denote ||o as max max(maxe (max || maxjejm) [[
3

---

Published as a conference paper at ICLR 2025
Further, let UEVT be the singular value decomposition (SVD) of A. Then, we have U R and Σ, V E Rm where U. V have orthonormal columns and Σ is a non-negative diagonal matrix. The Moore-Penrose pseudoinverse of matrix A is At VE -UT If E is a sorted diagonal matrix and
T1 ... Jm represent the diagonal entries of Σ, then we use aj to represent the i-th singular value of A, namely o; (A). We define miin (A) == min, ||Ax]|2// 1 / and m (A) max, | |ar||2
Now, suppose that m nt, namely A, W Raxx. We define ||A]| == V V Σラ= 2j=1 W A
Wo o A is matrix whose entries are defined as (WoA)ij Ai.. We define Dw: diag (W.)) If A is invertible, then the true inverse of A is denoted as A- and ||.A|| = m (Ain(AA))) If A is symmetric, then we define as UAUT T. the eigenvalue decomposition of A, where A diagonal matrix. Let ',,-.. An represent the entries on diagonal of A E R=*". X is called the i-th eigenvalue, namely (B). Furthermore, the eigenvalue and the singular value satisfy (A (AT A). Given two n. x real symmetric matrices A and B, we use A N B to denote the matrix B A is positive semidefinite, i.e., for any x E TR.", x ( A)x 2
2.2 SKETCHING
An important algorithmic subroutine the Subsampled Randomized Hadamard Transform SRHT:
Definition 2.1 (SRHT (Lu et al., 2013)). The SRHT matrix of size m x n is the following matrix: PHD,, where D E R= is diagonal matrix with diagonal being Rademacher random
S = -. 7mm
variables, H E RRxx is the Hadamard matrix and P E Rmxx is a row sampling matrix that samples m rows with replacement.
The key property we would like leverage from SRHT is the subspace embedding property:
Definition 2.2 (Oblivious subspace embedding (Sarlos, 2006)). Let d be positive integers and ê E (0, 1) be parameters, we say a distribution II over m x n real matrices satisfy (e, ó, n, d)- oblivious subspace embedding (OSE) if for any fixed matrix A E Boxd and S ~ II, with probability at least 1 - 6, we have for any x E Rdd
(1 - ) Ax|l2 < ||SA S (1 Ax]12
Via standard matrix concentration inequalities such as matrix Chernoff bound (see e.g. Rudelson (1999); Ahlswede Winter (2002)), one can show SRHT with m = O(e ld (log((n/)))) satisfying (e, 6, n, d)-OSE. Moreover, since H is Hadamard matrix, applying S to an n-dimensional vector can be done in On log n) using FFT. Thus, computing SA takes O(nd log n.) time.
2.3 BACKGROUND ON WEIGHTED Low RANK APPROXIMATION
The weighted low rank approximation can be treated as a generalization the noisy matrix completion problem, where the goal is recover target matrix M Rnxn from few observations (sublinear in n²) where the weight is chosen as a Boolean matrix Pn E RR**n It is hence natural to impose and generalize assumptions from matrix completion if we would like to obtain any provable guarantees. Following Li et al. (2016), we make three assumptions and we will justify them one by one.
Assumption 2.3. Given a noisy, possibly higher-rank observation M E Rax* such that M M N, where M the rank-k ground truth we want to recover and N is the noise matrix. We assume:
1. M is y-incoherent: Let M = UVVT T be its SVD, we assume
max(]|us. IV. 1. S jk
We use to denote the condition number of M Tmax (MI* 3/ m (M*).
2. Weight W has y-spectral gap all-1's matrix:
I - 1, || $ yn.
4

---

Published as a conference paper at ICLR 2025
3. Weight W is (a, (-)-boundedd Let M UEVT be its SVD, assume for any i [n]] and < a $ $ S.
al 1 U Dw: U 31, al A V Dw: V = BI.
Assumption I states that the largest row norms of the left and right singular factors should not be too far away from the average. Such matrix incoherence assumption has been very standard in the context of matrix completion (Candès Recht, 2012) as it effectively eliminates the degenerate case where the ground truth M* has very weak signals. Consider the extreme case where M eiee in such a scenario, the weight W rather uniform over all entries and N is a dense noise matrix with its first entry has small magnitude compared to other entries, then recovering M will be next to impossible. The incoherence assumption makes sure that the row and column space of 1 are spread over coordinates. Incoherence is also commonly observed in practice (Mohri & Talwalkar, 2011).
Assumption 2 is a natural generalization of the random sampling assumption for matrix comple- tion (Jain et al., 2013; Hardt, 2014). In particular, if W is a Boolean matrix where each has Mlogg n) entries chosen uniformly at random, then 7 O Vlogn). Generalize to a non-negative weight setting, it also bounds the largest possible magnitude of any entry in W to avoid degeneracy.
Assumption 3 is also best understood when W is Boolean matrix, so that Dw, selects subset of rows of U and V, and the condition essentially reduces to Assumption A2 of Bhojanapalli & Jain (2014). is a strengthening and weighted generalization of the strong incoherence property as directly implies the assumption in Candès & Tao (2010), which is necessary for matrix completion.
Having justified the assumptions we impose on the ground truth and the weight, we are in the position to state the weighted low rank approximation problem.
Problem 2.4. Let M E RRxx be a noisy, higher-rank matrix with M = M + N where M is the rank-k ground truth and N is higher-rank noise matrix. Let W E RRxx be non-negative weight matrix. Suppose both M and W satisfy Assumption 2.3. The goal is find rank-k matrix M E gaxa such that
IM 11 [oNN
by observing the matrix W a M.
When W is Boolean matrix, Problem 2.4 reduces to the noisy matrix completion problem where one needs to recover the rank-k ground truth by observing few entries of a higher-rank noisy matrix.
3 TECHNIQUE OVERVIEW
In this section, we provide a preliminary overview of the techniques we use in this paper. Before diving into our algorithm and analysis, let us first review the algorithm of Li et al. (2016). At each iteration, the algorithm alternates by solving two weighted multiple response regressions: starting with an initial matrix Y, it tries to find a matrix then they zero out the rows of X with large le norms and use the QR factor of X to proceed. Then, they alternate and solve miny out large rows and QR, the algorithm proceeds to the next iteration. The main runtime bottleneck is to solve the weighted multiple response regression per iteration.
E Ra that minimizes IIW o (M - (M-XY XY
| (M XY T I|? given the new X. After properly zeroing
Following the trend of low rank approximation (Clarkson & Woodruff, 2013) and fixed parameter tractable algorithm for weighted low rank approximation (Razenshteyn et al., 2016), is natural to consider using sketching to speed up the multiple response regression solves. Let consider
||W o (M XY )) YER
min
Let D D/w, denote the diagonal matrix that puts / /W, on the diagonal, where W, is the i-th column of W. It is not hard to verify that (1) can be cast into n linear regressions (see details in Claim C.I), each of which is in the form of
min I| yERE
VW. M - D VW, D.w,Yyl3
5

---

Ö(k) Published as a conference paper at ICLR 2025
To solve these regressions faster, one can pick random sketching matrix S R where ( 2k) and instead solve
min |SD H.D/WW,.. N., S----SSww... SD.w.Kyl@ yER
By picking a sparse sketching matrix S (Nelson & Nguyên, 2013), the above regression can be solved in O nnz(X (X) + Eo * time with high probability, and the output solution has cost at most (1+ 60) OPT where OPT is the optimal regression cost. Aggregate over n regressions, this gives an Ö( n nnz(X) + co 2 3) time per iteration (see Lemma C.6).
This approach, however, has several drawbacks that make it infeasible for application. The first is the error guarantee of such approximates regression solves. Essentially, we compute a matrix Y E Rank such that
| (M - XXTII S (1 66) min YEREXA
o o (M-XY
in other words, the approximate solution Y provides relative forward error. Unfortunately, the forward error is much less helpful when want analyze how close Y is the optimal solution Y. i.e,, the backward error. It is possible to convert forward error to backward error at the expense of dependence on other terms such as the cost of the regression and the spectral norm of X the pseudo- inverse of X. To cancel out the effect of these extra terms, we will have to set the error parameter eo to be very small, thus, a polynomial dependence on 6 in the running time unacceptable.
This motivates us to design a fast and high precision regression solver whose eo dependence is (og(L/co) (see Lemma C.10). Given an algorithm that produces an (1 + Eu) relative forward error of regression in log(1/ /eo) iterations, we can set 6g to inverse proportionally to OPT [[ Xy*I- As the spectral norm of (W x)1 polynomially bounded, this incurs an extra log n term in the runtime. It remains to devise regression solver with such runtime behavior. Our approach is to preconditioner: we pick dense sketching matrix S E R=*= with
use the sketch rows such that for any k-dimensional vector xr ||Sti||2 = (1 + O(1) || We then apply $ to D D.w, X to form short and fat matrix and compute the QR decomposition of this matrix. It turns out that the right QR factor of SD - a constant approximation the regression problem and utilize as starting point. The algorithm then iteratively performs gradient descent to optimize towards an an-approximate solution. Overall, such an algorithm takes log(1/co) iterations to converge, and each iteration can be implemented in (nkk) time. Plus the extra Onk + I.²) time to compute the initial solution, this yields an algorithm that runs in Ö(Ank + X-3) log(1/co) /co)) time to compute an cu forward error solution. Note here we sacrifice the input sparsity time in exchange of sketching matrix that works with high probability. This also accounts for the fact that both X and Y are quantities changed across iterations and the sparsity cannot be controlled.
X is a good preconditioner to D X. We then use S to find
The runtime can be further improved by leveraging the sparsity of the weight matrix W. Again, consider the regression min |DD.W, M., - D D.w, Xylž if W; only has a few nonzero entries, then the diagonal matrix D Vw: will effectively zero out most rows of X and entries of M.,. This means that are solving regression of size OW (luk) instead of (nk)). As iterate through all n regressions, the total instance size is then O (W.llok) = ([okk))kk) and we can effectively solve these regressions in an overall O((IIW + nkk²) log(1/eo)) time. We note that in matrix completion, [ l/o is oftentimes Ö(n pkly(k)). making much smaller than O(?²) and an algorithm that exploits its sparsity is therefore much more valuable.
We want remark that our high precision and dense regression solver not only works for weighted low rank approximation, but for any alternating minimization frameworks that require one to solve (11) multiple response regressions per iteration. Due to the good error dependence. the overall log(1/e) /e) convergence is well-preserved even though each iteration is only solved approximately. We believe this high precision solver will also find its use problems like (low rank) matrix sensing and tasks in which backward error for multiple response regression is required.
In addition to our high-accuracy, high probability solver, we also devise a robust analytical framework for alternating minimization, which is the core to enable us with fast approximate solvers. In particular, show that if we only output matrix y that is close the exact regression solution Y in the

---

Published as a conference paper at ICLR 2025
spectral norm, then the alternating minimization still converges the fixed point I| MII with good speed. Our analysis uses different strategy from Li et al. (2016) where they heavily rely on the closed-form of the regression solution. In contrast, errors, one can accumulate the error caused by approximate solves to the additive e term and thus polylogarithmically
show that by clever decomposition of
rounds of the iterative solve suffices give good guarantees. When adapting our analysis to noiseless matrix completion (Gu et al., 2024b), we recover their result in both runtime and sample complexity, while offering a much simpler proof.
OUR RESULTS
In Section 4.1, we analyze the weighted multiple response regression. In Section 4.2, show that the alternating minimization framework is robust, namely this alternating minimization framework can tolerate the error induced by the approximate solver and error conversion. In Section 4.3, present the formal version of our main result. Finally, in Section 4.4, we compare our results and contribution with those of prior works.
4.1 WEIGHTED MULTIPLE RESPONSE REGRESSION
One of our cornerstone results is novel adaptation of a high-accuracy regression solver based on sketching. Its root can be perhaps traced back to Rokhlin & Tygert (2008), and our two new insights are: 1). This type of high-accuracy regression solvers also be generalized to weighted case, where the design matrix and target vector are scaled by some non-negative weights. 2). We can convert the the cost of the regression the error on the solution. This step is crucial, as bridge the gap between our fast, approximate solves and the exact solutions used in Li et al. (2016), it is essentially quantifyy the difference between solutions.
Lemma 4.1. Let A Raxd b E R" and 18 Let 0.1) be accuracy parameter and 5 E (0, 1) be the failure probability. Suppose T(n, ed,,),,))))).)) ) is the runtime of black-box regression solver that produces a vector x' E Id such that
|| bl22 (1 min
|Ax
bl[2
ZERd
with probability at least Then, there exists an algorithm that runs in time ((nz((A) + T(n, d, 6)
and outputs a vector x Bd such that with probability at least 1-60
|| - bllw S (1 + +e) min | - bllu-
DERd
The proof relies on simple observation: the weights could be applied by scaling rows of A and entries of b, which in turn could be implemented in nearly linear time. This simple reduction allows us deploy fast off-the-shelf regression solver for weighted regression. To facilitate the analysis, we also require a conversion from the regression cost to how close our approximate solution is to the optimal solution.
Lemma 4.2. Let A E Baxd with n 2 d and full rank, b E Rh and let xoPT be the exact solution to the regression problem minzers ||Ax - bl[2. Suppose there exists a vector x' Rdd with
|| bll2 I (1 + AxOPT - b112.
then we have
1
| ZOPT|12 $ <((Ee)
. ||AxOPT - bl[2-
Tmin (A)
The conversion from forward to backward error is standard (Price et al., 2017; Gu et al., 2024b), and means that we will have solution. We combat this issue by employing a high-accuracy regression solver.
be polynomially small in amin ( and the cost of the optimal
The rough idea behind Algorithm is compute a quick preconditioner using sketching. Let S E mm be SRHT matrix with 2 (n/))))))) rows, it is (0.01, 8, n, OS-OSE,
7

---

Published as a conference paper at ICLR 2025
Algorithm 2 High precision solver.
1: procedure HICHIPRECISIONREGCA E Ruxd be DERR IR" (0, (0, e (0, > Lemma C.10
2:
c1 < 0.01
3:
m + - d
4:
Let S E RRxxa an SRHT matrix
5:
Compute QR decomposition of 5 QR-
6:
To f arg eRd ||S |SAR Sbll2
7:
T C (L) for sufficiently large constant C
8:
for = 0 do
9:
T++1 21 R (T (R+:)
ARx,)
10: end for
return Rer
12: end procedure
therefore with high probability, the singular values of S'A are close to A. The QR decomposition of SA provides an orthonormal basis Q and non-singular upper triangular matrix R-1 which serves as a good preconditioner for A. We can then proceed with preconditioned gradient descent using R. This procedure is particularly fast because the most tonecconumingg step to compute the QR decomposition, but is performed on an d matrix. Further, SA can be carried out in nearly linear time, and all subsequent steps in gradient descent can performed in manner that takes nearly linear time. The property of SRHT also ensures our initial point is constant approximation of the optimal point, therefore the algorithm converges in ((log(I/e) iterations, as desired. For more details, refer readers to Appendix C.
Lemma 4.3. Given a matrix A E gaxd and a vector b R", let e (0,0.1) and å E (0, 0.1,, there exists an algorithm that takes time
O(ndlogn + d log (n/0)) log(L/e))
and outputs x' E ERd such that
|| - -bl2 S (1 c)
min AA bll2 DERd
holds with probability 1
4.2 ROBUSTNESS ANALYSIS FOR APPROXIMATE UPDATE
Now that we have the regression solvers that can compute an approximate update in nearly linear time, we need show that the alternating minimization framework robust enough tolerate the large error induced by the approximate solver. We introduce generalized incoherence notion. Definition 4.4. Let max 11:
Ra*k, we define the generalized incoherence of A as p(A) 2
As our analysis crucially exploits the interplay between exact and approximate updates, we summarize the notations in the following table simplify the discussion.
Table 1: Summarization of notations regarding exact and approximate regression solves. By "clipped", we mean zeroing out rows with large l norms.
Notation
Meaning
X
Matrix for exact regression solve
X
Clipped matrix of X
X
QR factor of X = XR
X Matrix for approximate regression solve
X Clipped matrix of X
Lemma 4.5. Let Y BRnk be matrix with orthonormal columns and & and 6st be parameters and D2 be a parameter depends on &, Ea Let X. X, X, and X be defined as in Table with the clipping threshold being 46. Moreover, have I| f. XI. Este/ Finally, let M U.VT.
8

---

Published as a conference paper at ICLR 2025
Then, we have
j M S A:
Au S 0. (M*), then (i, I 8 (M));
. 1fAu S 0. 0.1 (M*), then p(X) S 8p//mmn(MM")
where dist(U, minge0x*x IUQ - X|| where Okxk the set of all k orthogonal matrices.
Let us interpret the above lemma and explain why it's crucial to our final convergence analysis. For simplicity, suppose the noise N and M X-YT. since Y has orthonormal columns, M = XX* and the first part states that if we solve the regression approximately and clip rows with large norms, then the approximate clipped matrix X is close to X The next two parts state that long as X and X close enough, then two crucial properties between the space spanned by left singular vectors and X, the QR factor the clipped matrix X, is small and 2). the generalized incoherence of X is small. These guarantees lead to natural inductive argument: suppose A, small enough, then by our algorithm, we know that X and M are close and consequently dist XX, U) and x(X) are small. These two conditions serve as a basis to prove that for the next iteration, still have X and M is small enough and the induction can proceed.
guaranteed: 1). the distance
We want to highlight the major challenges in proving these assertions. Note that the induction argument effectively provides bounds on both subspace distance and generalized incoherence, and both notions heavily rely on the conditioning of intermediate matrices. The original analysis of Li al. (2016) gives quantitative bounds on condition numbers assuming the updates are computed exactly, but the picture becomes much less clear when the updates are only computed approximately. Nevertheless, we prove that when the approximate updates are close enough to the optimality, then these bounds still hold. To compute these updates to high-precision, we utilize the high-accuracy. weighted multiple response solver being developed. One could view our proof as mixture of algorithm and analysis: our analysis mandates the algorithm to provide strong guarantees, and we in turn design algorithms to achieve these goals. For more details, we refer readers to Appendix D.
4.3 MAIN RESULT
Our main theorem is as follows:
Theorem 4.6 (Formal version of Theorem 1.1). Given a noisy, possibly higher-rank observation M e Raxa where M is the rank-k ground truth and N is the noise matrix that satisfies Assumption 2.3. There is an algorithm (Algorithm 1) uses random initialization, runs in (log((/c)) iterations and matrix M such that
generates
I
M 1 O N[|+G
The total running time is
O(II k n2.3) log(1/e)).
Due space limitation, delay the proof of Theorem 4.6 Appendix J. We briefly remark that our algorithm can be easily extended to cases where both W and M are rectangular matrices of size X n as none of our analyses rely on the matrix being square. One could replace the factor n in our runtime by max(m (m, n when dealing with rectangular weighted low rank approximation.
4.4 COMPARISONS WITH RECENT WORKS
In this section, we provide brief overview and comparison with other recent works, which could be classified into 3 categories: 1). slower, exact alternating minimization for weighted low rank approximation (Li et al., 2016); 2). faster, approximate alternating minimization for noiseless low rank matrix completion, strictly simpler problem (Gu et al., 2024b) and 3). new metries for measuring the effectiveness of low rank matrix factorizations (Yalcin et al., 2022; Zhang et al., 2024).
Compared to the result of Li et al. (2016), we significantly improve the running time from O([[ nn-²) log()) (1)) to O([I]W]lkk + log)) log(1/e) For moderately large k (say

---

Published as a conference paper at ICLR 2025
k OVVT)) and dense weight matrix (say || = Dn,)), the Li et al. (2016) algorithm would take Onn log(1/e)) time, while ours only takes Ō(n25 log(L/e)) time. In the noisy ma- trix completion setting Kelner et al. (2023), W is a Boolean matrix with || = ((n)2-(11) applying our algorithm leads to an overall runtime of ((nk3+ol1) log(1/c)), nearly matches the tfte-of-thee-rt (Kelner et al., 2023) In contrast, the Li et al. (2016) algorithm has runtime of O(nk++ol)) log(1/e)). Moreover, our analysis accounts for the approximated computation at each step, thus it opens up the gate for further speedup. This also better depicts the picture of practical alternating minimization algorithms, where updates are computed approximately both due to floating point errors and efficiency concerns. We believe this also lays foundation for theoretically verifying why alternating minimization with approximate updates has great empirical success.
Compared to the result of Gu et al. (2024b), note they show for the simpler problem of noiseless matrix completion, the alternating minimization procedure could be sped up and run in Ö(1 W llo k (og((/()) time. Our result, even in the matrix completion setting, is a strict generalization of theirs, as they assume access to the entries of the ground truth M* In contrast, our model can only access noisy entries M = M + NN thus our recovery result suffers an error in the form of O(T)) IIW NNI. which is O if = xn We also provide spectral norm error guarantee rather than Frobenius norm error, which is the objective Gu et al. (2024b) obtains. The spectral norm is oftentimes considered more robust than the Frobenius norm. In terms of analysis, the proof approach of Gu et al. (2024b) is particularly geared towards noiseless matrix completion, while analysis is much more general, as it can account for noisy matrix completion and weighted low rank approximation. We believe the generality and simplicity of our framework could be further extended to analyze alternating minimization for other problems, such as robust PCA and multi-view learning.
Compared the results of Yalein al. (2022): Zhang al. (2024), we focus on providing theoretical guarantees on the algorithm's 's performance, whereas Yalcin et al. (2022); Zhang et al. (2024) focus more on analyzing the optimization landscape and proposing complexity metries for low rank matrix problems than developing specific algorithms for weighted low rank approximation. Specifically, the main contribution of Zhang et al. (2024) is developing new complexity metric to characterize the difficulty of the nonconvex landscape arising from the Burer-Monteiro factorization. This metric aims to quantify when local search methods can successfully solve the factorized problem. The main contribution of Yalcin et al. (2022) is constructing class of low-complexity matrix completion problem instances that can be solved in polynomial time, but for which the popular Burer-Monteiro factorization approach fails. Yalcin et al. (2022) also shows the existence of problem instances in this class that have exponentially many spurious local minima when using the factorization, even though the original problem has unique global solution. It would be interesting to study whether alternating minimization could also provide provable guarantees against these metrics and problems and in turn be accelerated.
5 CONCLUSION
In this paper, we study the weighted low rank approximation problem and efficient algorithm to solve it under mild recovery assumptions. Alternating minimization has been shown to be a powerful algorithmic prototype for this problem (Li et al., 2016), and we provide a fast, approximate implementation together with a robust error analysis for the framework. To this end, we improve the
running time of Li et al. (2016) from [[[[ lluk² + n)=²) log(1/c) 1c)) to W[]kk + log()) (1/c)). Our error analysis also serves as a theoretical explanation of why alternating minimization works well in practice especially when these updates are computed approximately for better efficiency.
We would also like to point out that the runtime of our algorithm is nearly linear in terms of solution verification. Given the weight matrix W and a pair of low rank factors X and Y, it takes ((k) to verify a single entry of W (XYT) and we would need to verify a total of [|W]| entries. However, it is also worth noting that such runtime can only be achieved when random initialization is used as if one resorts to SVD initialization, the initialization time becomes O which would dominate the overall runtime. It will be an interesting open problem whether we further speed up the initialization using procedures such as random SVD and obtain a nearly linear time algorithm for alternating minimization with SVD initialization.
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGEMENT
This work was mostly done when Junze Yin was at Boston University. Junze Yin is supported by the Rice University graduate fellowship. Lichen Zhang is supported in part by NSF CCF-1955217 and NSF DMS-2022448.
REFERENCES
Rudolf Ahlswede and Andreas Winter. Strong converse for identification via quantum channels. IEEE Transactions on Information Theory, 2002.
Josh Alman and Zhao Song. Fast attention requires bounded entries, Advances in Neural Information Processing Systems, 36, 2023.
Josh Alman and Zhao Song. The fine-graing grained complexity of gradient computation for training large language models. In The Thirtty-ighuth Annual Conference on Neural Information Processing Systems, 2024. URL https::/p//o wpwnreveevveveeeeeee ew net/ Forum? tWwwol..........
Josh Alman, Jiehao Liang, Zhao Song, Ruizhe Zhang, and Danyang Zhuo Bypass exponential time preprocessing: Fast neural network training via weight-data correlation preprocessing. In Advances in neural information processing systems, NeurIPS'22 23, 2023.
Sanjeev Arora, Yuanzhi Li, Yingyu Liang, Tengyu Ma, and Andrej Risteski. A latent variable model approach to pmie word embeddings. Transactions of the Association for Computational Linguistics, 4:385-399, 3999,,9,.20 2016.
Frank Ban, David Woodruff, and Richard Zhang. Regularized weighted low rank approximation. In Advances in neural information processing systems (NeurIPS), 2019.
Aditya Bhaskara, Aravinda Kanchana Ruwanpathirana, and Maheshakya Wijewardena. Additive error guarantees for weighted low rank approximation. In Marina Meila and Tong Zhang (eds.), Pro- ceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pp.
PMLR, 2021.
Srinadh Bhojanapalli and Prateek Jain. Universal matrix completion. In Eric Xing and Tony Jebara (eds.), Proceedings of the 31 International Conference on Machine Learning, volume 32 of Proceedings of Machine Learning Research, pp. 188 1881-, Bejing, China, 22-24 Jun 2014. PMLR.
Song Bian, Zhao Song, and Junze Yin. Federated empirical risk minimization via second-orde method. arXiv preprint arXiv:2305.17480, 2023.
Christos Boutsidis and David P Woodruff. Optimal cur matrix decompositions. In Proceedings of the forty-sixth annual ACM symposium Theory of computing (STOC), pp. 353-362, 2014.
Christos Boutsidis, David Woodruff, and Peilin Zhong. Optimal principal component analysis in of the 48th Annual ACM SIGACT
distributed and streaming models. In Symposium Theory of Computing, 2016.
Emmanuel Candès and Benjamin Recht. Exact matrix completion via convex optimization. Commun. ACM, 2012.
Emmanuel J. Candès and Terence Tao. The power of convex relaxation: Near-optimal matrix completion. IEEE Trans. Inf. Theor, 2010.
Moses Charikar, Kevin Chen, and Martin Farach-Colton. Finding frequent items in data streams. In International Colloquium on Automata, Languages, and Programming, pp. 693-703. Springer, 2002.
Chao Chen, Dongsheng Li, Yingying Zhao, Qin Lv, and Li Shang. Wemarec: Accurate and scalable recommendation through weighted and ensemble matrix approximation. In Proceedings of the 38th international ACM SIGIR conference on research and development in information retrieval, pp. 303-312, 2015.
11

---

Published as a conference paper at ICLR 2025
Yifang Chen, Jiayan Huo, Xiaoyu Li, Yingyu Liang, Zhenmei Shi, and Zhao Song. Fast gradient computation for rope attention in almost linear time. arXiv preprint drXiv:2412.1731,, 2024.
Kenneth L Clarkson and David Woodruff. Low-rank approximation and regression in input sparsity time. In STOC, 2013.
Michael B Cohen. Nearly tight oblivious subspace embeddings by trace incqualities. In Proceedings of the twenty-seventh annual ACM-SIAM symposium on Discrete algorithms, pp. 278-287. SIAM, 2016.
Michael B Cohen, Yin Tat Lee, and Zhao Song. Solving linear programs in the current matrix multiplication time. Journal of the ACM (JACM), 68(1):1-39, 2021.
Yucheng Dai. On algorithms for weighted low rank approximation. Master's thesis, Carnegie Mellon University Pittsburgh, PA, 2023.
Yichuan Deng, Zhao Song, and Omri Weinstein. Discrepancy minimization in input-sparsity time. arXiv preprint arXiv:2210. 1246,, 2022.
Yichuan Deng. Zhihang Li, and Zhao Song. Attention scheme inspired softmax regression. arXiv preprint ar 11111:2304.1011,, 2023a.
Yichuan Deng, Zhao Song, and Junze Yin. Faster robust tensor power method for arbitrary order. arXiv preprint arXiv:2306.0044 2023b.
Huaian Diao, Zhao Song, Wen Sun, and David Woodruff. Sketching for kronecker product regression and p-splines. In International Conference Artificial Intelligence and Statistics, pp. 1299-1308. PMLR, 2018.
Huaian Diao, Rajesh Jayaram, Zhao Song, Wen Sun, and David Woodruff. Optimal sketching for kronecker product regression and low rank approximation. Advances in neural information processing systems, 32, 2019.
Yeqi Gao, Zhao Song, Weixin Wang, and Junze Yin. fast optimization view: Reformulating single layer attention in Ilm based tensor and svm trick, and solving it in matrix multiplication time. arXiv preprint htt.:2309900411.. 2023a.
Yeqi Gao, Zhao Song, and Junze Yin. Gradientcoin: A peer-to-peer decentralized large language models. arXiv preprint arXiv:2308.1 10502, 2023b.
Yeqi Gao, Zhao Song, and Junze Yin. An iterative algorithm for rescaled hyperbolic functions regression. The 28th International Conference on Artificial Intelligence and Statistics, 2025. URL http::ttp/// /ttp::://ppeeee wpwprevieeeveeeeeeeeeeeeee net. htte:/ooommmmmxx0000e00
Nicolas Gillis and François Glineur. Low-rank matrix approximation with weights missing data is np-hard. SIAM Journal on Matrix Analysis and Applications, 32(4):1149-116 2011.
Yuzhou Gu and Zhao Song. A faster small treewidth sdp solver. arXiv preprint arXiv:2211.06033, 2022.
Yuzhou Gu, Nikki Lijing Kuang, Yian Ma, Zhao Song, and Lichen Zhang. Log-concave sampling from a convex body with barrier: a robust and unified dikin walk. In The Thirty-eighth Annual Conference Neural Information Processing Systems, 2024a. URL httpt: /tttps://openrevieee net/ feet/foze forum? xd=xXr
Yuzhou Gu, Zhao Song. Junze Yin, and Lichen Zhang. Low rank matrix completion via robust alternating minimization in nearly linear time. In Proceedings of the 12th International Conference on Learning Representations, ICLR 24, 2024b.
Yuzhou Gu, Zhao Song, and Lichen Zhang. Faster algorithms for structured linear and kernel support vector machines. In The Thirteenth International Conference Learning Representations, 2025. URL https: / /ttp:///peeeeeeee iew net/ forum? Ld=DDNFTaVQdU.
Moritz Hardt. Understanding alternating minimization for matrix completion. In 2014 IEEE 55th Annual Symposium Foundations of Computer Science, pp. 651-660. IEFF, 2014.
12

---

Published as a conference paper at ICLR 2025
Yen-Chang Hsu, Ting Hua, Sungen Chang, Qian Lou, Yilin Shen, and Hongxia Jin. Language model compression with weighted low-rank factorization. arXiv preprint arXiv:2207.00112, 2022.
Jerry Yao-Chieh Hu, Erzhi Liu, Han Liu, Zhao Song, and Lichen Zhang. On differentially private string distances. arXiv preprint arXiv:2411.05750 2024.
Prateek Jain, Pranecth Netrapalli, and Sujay Sanghavi. Low-rank matrix completion using alternating minimization. In Proceedings of the Forty-Fifth Annual ACM Symposium Theory of Computing, STOC '13 New York, NY, USA, 2013.
Shunhua Jiang, Zhao Song, Omri Weinstein, and Hengjie Zhang. A faster algorithm for solving general lps. In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing. pp. 823-832, 2021.
Jonathan Kelner, Jerry Li, Allen Liu, Aaron Sidford, and Kevin Tian. Matrix completion in almost- verification time. In 2023 IEEE 64th Annual Symposium on Foundations of Computer Science, FOCS' 23, 2023.
Lukasz Kidzinski, Francis KC Hui, David Warton, and Trevor J Hastie. Generalized matrix factorization: efficient algorithms for fitting generalized linear latent variable models to large data arrays. Journal of Machine Learning Research, 23(291):1-29, 2022.
Yehuda Koren, Robert Bell, and Chris Volinsky. Matrix factorization techniques for recommender systems. Computer, 42(3):30-37, 2009.
Joonseok Lee, Seungyeon Kim, Guy Lebanon, and Yoram Singer. Local low-rank matrix approxima- tion. In International conference on machine learning, pp. 82-90. PMLR, 2013.
Yin Tat Lee, Zhao Song, and Qiuyi Zhang. Solving empirical risk minimization in the current matrix multiplication time. Conference on Learning Theory, pp. 2140-2157. PMLR, 2019.
Chenyang Li, Zhao Song, Zhaoxing Xu, and Junze Yin. Inverting the leverage score gradient: An efficient approximate newton method. arXiv preprint arXiv:2408.11200 2024.
Yuanzhi Li, Yingyu Liang, and Andrej Risteski. Recovery guarantee of weighted low-rank approx- imation via alternating minimization. In International Conference Machine Learning, pp. PMLR, 2016
Zhihang Li, Zhao Song, Zifan Wang, and Junze Yin. Local convergence of approximate newton method for two layer nonlinear regression. arXiv preprint arXiv:2311.1539) 2023a.
Zhihang Li, Zhao Song, and Tianyi Zhou. Solving regularized exp, cosh and sinh regression problems. arXiv preprint, 2303.15725. 2023b.
Jiehao Liang, Somdeb Sarkhel, Zhao Song, Chenbo Yin. Junze Yin, and Danyang Zhuo A faster k-means++ algorithm. arXiv preprint arXiv:2211.55118 2022.
Yingyu Liang, Heshan Liu, Zhenmei Shi, Zhao Song. Zhuoyan Xu, and Junze Yin. Conv-basis: A new paradigm for efficient attention inference and gradient computation in transformers. arXiv preprint arXiv:2405. 05219, 19, 2024.
W-S Lu, S-C Pei, and P-H Wang. Weighted low-rank approximation of general complex matrices and its application in the design of 2-d digital filters. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, 44(7))650 655, 1997.
Wu-Sheng Lu and Andreas Antoniou. New method for weighted low-rank approximation of complex- valued matrices and its application for the design of 2-d digital filters. In Proceedings of the 2003 International Symposium Circuits and Systems, 2003. ISCAS'03., volume 3, pp. III-III. IEEE, 2003.
Yichao Lu, Paramveer Dhillon, Dean P Foster, and Lyle Ungar. Faster ridge regression via the subsampled randomized hadamard transform. In Advances in neural information processing systems (NIPS), pp. 369-377, 2013.
13

---

Published as a conference paper at ICLR 2025
Ivan Markovsky, Maria Luisa Rastello, Amedeo Premoli, Alexander Kukush, and Sabine Van Huffel. The element-wise weighted total least-squares problem. Computational Statistics and Data Analysis, 50(1):181-205 2006. 2nd Special issue on Matrix Computations and Statistics.
Mehryar Mohri and Ameet Talwalkar. Can matrix coherence be efficiently and accurately estimated? In Proceedings of the Fourleenth International Conference Artificial Intelligence and Statistics, Proceedings Machine Learning Research, 2011.
Jelani Nelson and Huy L Nguyên. Osnap: Faster numerical linear algebra algorithms via sparser subspace embeddings. In 2013 ieee 54th annual symposium on foundations of computer science, pp. 117-126. IEEE, 2013.
Jeffrey Pennington, Richard Socher, and Christopher Manning. Glove: Global vectors for word representation. In Proceedings of the 2014 conference on empirical methods in natural language processing (EMNLP), pp. 1532-1543, 2014.
Eric Price, Zhao Song, and David P Woodruff. Fast regression with an foo guarantee. ICALP. 2017.
Lianke Qin, Zhao Song, Lichen Zhang, and Danyang Zhuo An online and unified algorithm for projection matrix vector multiplication with application to empirical risk minimization. In AISTATS, 2023.
Ilya Razenshteyn, Zhao Song, and David P Woodruff. Weighted low rank approximations with provable guarantees. In Proceedings of the forty-eighth annual ACM symposium on Theory of Computing, pp. 250-263, 2016.
Aravind Reddy, Zhao Song, and Lichen Zhang. Dynamic tensor product regression. In NeurIPS, 2022.
Geneviève Robin, Julie Josse, Éric Moulines, and Sylvain Sardy, Low-rank model with covariates for count data with missing values. Journal of Multivariate Analysis, 173:416-434, 2019.
Vladimir Rokhlin and Mark Tygert. A fast randomized algorithm for overdetermined linear least- squares regression. Proceedings of the National Academy of Sciences, 105(36): 132122))::322200 13217, 2008.
M. Rudelson. Random vectors in the isotropic position. Journal of Functional Analysis, 164(1): 60-72, 1999.
Tamas Sarlos. Improved approximation algorithms for large matrices via random projections. In 2006 47th annual IEEE symposium on foundations of computer science (FOCS'06), pp. 143-152. IEEE, 2006.
Dale J Shpak. Beighted-least-qquaar matrix decomposition method with application the design of two-dimensional digital filters. In Proceedings of the 33rd Midwest Symposium on Circuits and Systems, pp. 1070-1073. IEEE, 1990.
Anshumali Shrivastava, Zhao Song, and Zhaozhuo Xu. Sublinear least-squares value iteration via locality sensitive hashing. In AISTATS. arXiv preprint arXiv:2105.08285, 2023.
Ritwik Sinha, Zhao Song, and Tianyi Zhou. A mathematical abstraction for balancing the trade-off between creativity and reality in large language models. arXiv preprint arXhv:2366.02295, 2023.
Zhao Song and Zheng Yu. Oblivious sketching In International Conference on Machine Learning, pp. 9835-9847. PMLR, 2021.
central path method for linear programming.
Zhao Song, David P Woodruff, and Peilin Zhong. Low rank approximation with entrywise 11-norm еттor. In Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing, pp. 688-701, 2017.
Zhao Song, David Woodruff, and Peilin Zhong. Towards a zero-one law for column subset selection. Advances in Neural Information Processing Systems, 32, 2019a.
14

---

Published as a conference paper at ICLR 2025
Zhao Song, David Woodruff, and Peilin Zhong. Average column subset selection for entrywise no-nm loss. Advances Neural Information Processing Systems, 32, 2019b.
Zhao Song, David P Woodruff, and Peilin Zhong. Relative enror tensor low rank approximation. In Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pp. 2772-2789. SIAM, 2019c.
Zhao Song, David Woodruff, Zheng Yu, and Lichen Zhang. Fast sketching of polynomial kernels of polynomial degree. International Conference on Machine Learning, pp. 9812-9823, PMLR.
2021a.
Zhao Song, Shuo Yang, and Ruizhe Zhang. Does preprocessing help training aver-parameterizer neural networks? Advances in Neural Information Processing Systems, 14.22890-22906 2021b.
Zhao Song, Lichen Zhang, and Ruizhe Zhang. Training multi-layer over-parametrized neural network in subquadratic time. arXiv preprint 12X:2:21 2021c
Zhao Song, Zhaozhuo Xu, Yuanyuan Yang, and Lichen Zhang. Accelerating frank-wolds algorithm using cow-dimensional and adaptive data structures. arXiv preprint vtXiv:2207.09002, 2022a.
Zhao Song, Zhaozhuo Xu, and Lichen Zhang. Speeding up sparsification using inner product search data structures. arXiv preprint arXiv:2204.03209 2022b.
Zhao Song, Xin Yang, Yuanyuan Yang, and Tianyi Zhou. Faster algorithm for structured john ellipsoid computation. arXiv preprint arXiv:2211.14410 2022c.
Zhao Song, Weixin Wang, and Junze Yin. A unified scheme of resnet and softmax. arXiv preprint MiXiv:2309.13482. 2023a.
Zhao Song, Yitan Wang, Zheng Yu, and Lichen Zhang. Sketching for first order method: Efficient algorithm for low-bandwidth channel and vulnerability. ICML, 2023b.
Zhao Song. Xin Yang. Yuanyuan Yang, and Lichen Zhang. Sketching meets differential privacy: fast algorithm for dynamic kronecker projection maintenance. In International Conference on Machine Learning, pp. 32418-32462. PMLR, 2023c.
Zhao Song, Mingquan Ye, Junze Yin, and Lichen Zhang. A nearly-optimal bound for fast regression with foo guarantee. In International Conference Machine Learning, pp. 32463-32482. PMLR, 2023d.
Zhao Song, Mingquan Ye, and Lichen Zhang. Streaming semidefinite programs: o(yn) passes, small space and fast runtime. arXiv preprint arXiv:2309.05115, 2023e.
Zhao Song, Junze Yin, and Ruizhe Zhang. Revisiting quantum algorithms for linear regressions: Quadratic speedups without data-dependent parameters. arXiv preprint 1rXiv:2311. 148222, 2023f.
Zhao Song, Junze Yin, and Lichen Zhang. Solving attention kernel regression problem via pre- conditioner. In International Conference on Artificial Intelligence and Statistics, pp. PMLR, 2024.
Zhao Song, Weixin Wang, Chenbo Yin, and Junze Yin. Fast and efficient matching algorithm with deadline instances. In The Second Conference on Parsimony and Leurning (Proceedings Track), 2025. URL http:///pppeeeieeeeeeeee net. fest//form? d== Ine.Xrre..
Nathan Srebro and Tommi Jaakkola. Weighted low-rank approximations. In Proceedings of the 20th international conference machine learning (ICML), pp. 720-727, 2003.
Elena Tuzhilina, Trevor Hastie, and Mark R Segal. Principal curve approaches for inferring 3d chromatin architecture. Biostatistics, 23(2):626-66 642, 2022.
Peter D. Wentzell, Darren T. Andrews, David C. Hamilton, Klaas Faber, and Bruce R. Kowalski. Maximum likelihood principal component analysis. Journal of Chemometrics, 11, 1997.
15

---

Published as a conference paper at ICLR 2025
Zhaozhuo Xu, Zhao Song, and Anshumali Shrivastava. Breaking the linear iteration cost barrier for some well-known conditional gradient methods using maxip data-structures. Advances Neural Information Processing Systems, 34:5576-5589 2021.
Baturalp Yalcin, Haixiang Zhang, Javad Lavaei, and Somayeh Sojoudi. Factorization approach for low-complexity matrix completion problems: Exponential number of spurious solutions and failure of gradient methods. In International Conference on Artificial Intelligence and Statisties, pp. 319-341. PMLR, 2022.
Gale Young. Maximum likelihood estimation and factor analysis. Psychometrika, 6(1):49-53, 1941.
Haixiang Zhang, Baturalp Yalcin, Javad Lavaei, and Somayeh Sojoudi. A complexity metric for nonconvex rank-one generalized matrix completion. Mathematical Programming, 207(1):227-268, 2024.
Lichen Zhang. Speeding up optimizations via data structures: Faster search, sample and maintenance. Master's s thesis, Carnegie Mellon University, 2022.
16

---

Published as a conference paper at ICLR 2025
APPENDIX
Roadmap. In Section A, we provide several basic definitions and tools. In Section B, we discuss more related work. In Section C, we deseribe the fast multiple response regression solver used to speed up the alternating minimization step. In Section D, we provide our key lemmas for the update step. In Section we prove our induction lemma. In Section we state several tools from previous work. In Section G, we analyze the SVD initialization and present our main result. In Section H, we present the random initialization algorithm and analyze its properties. In Section I, we show how to prove the final guarantee of our main Theorem. In Section J, we present the complete proof of our main theorem.
A BASIC DEFINITIONS AND ALGEBRA TOOLS
In Section A.1, we present the properties of the weight matrix. Moreover, explain the algebra tools which are used for later proofs. In Section A.2, we present some basic algebraic inequalities. In Section A.3, we state few simple facts about the norm properties.
A.1 PROPERTIES WEIGHT MATRIX
Here, we present the properties of weighted matrices.
Definition A.1. For weight matrix W, we define
W |o,1 == max (max [W.,lh, max [[.hh
(e(n)
jein]
and
| == max max W., 1, max [|W5l
iE/n]
jeIn
Lemma A.2. Let 7 > 0, if |[ - 1n1 S ym, then we have
Part |W -W-1,11] S 1n
Further |||| F I n 1, + n
Part 2. |W - 1 1o0,1 S n 1,7
Further [W|lo,, S n 1. +
Proof. Proof of Part 1. We have
IW - 1, 1, S n||W - 11
< n (m²)
< n 322
where the first step follows from Part 4 of Fact A.7, the second step follows from the assumption from the lemma statement, and the last step follows from simple algebra.
Moreover, by the triangle inequality, have
|[]] IW 1 + 1
< W 1, 1,111 || + 1
< n 1, 1T
1.5g n,
where the first step follows from simple algebra, the second step follows from the triangle inequality, the third step follows from Eq. (2), and the last step follows from the definition of the Frobenius norm.
17

---

Published as a conference paper at ICLR 2025
Proof of Part 2. Note that
I[ 1n
S vn- |W - -1,1, 20.2
S Vū || -1,1,
$ n 'n
By the triangle inequality, we have
W S I 1, boo, +
1.5 S n + n.
Lemma A.3. Bounds on ry lead to bounds on | Specifically,
Part 1. If 1/(10m)/)), then we have
7. [[[ (10
Part 2. If < 1/1 then we have
7 (loo.) 1/2
Remark A.4. In previous work (Li et al., 2016), they wrote the final bound as 7 [/([[W]l...//n) where f are factors not depending on For example, f poly(a Tt)) Their bound technically not complete, because | also function yf So, in our work, our Lemma A.3 further calculates the actual condition required by Li et al. (2016) and hence completes their correctness proof.
Proof. Proof of Part 1. We need that
7 ((]W
suffices to show that
7- 7-(nn 15 7 n//n)1//3
The above equation is equivalent to is sufficient to show that
(n 7 + 1)1122 <1
1.5 0.25
<
Thus, as long as the promised bound is held.
7 < 1/(10n2/6)
Proof of Part We need that
7. (- <1
suffices to show that is sufficient to show that
7. (n + 1/2
u.75 + 75 <1
Thus, as long as have the desired result.
< 1/(10n2/2)
18

---

Published as a conference paper at ICLR 2025
A.2 BASIC ALGEBRAIC INEQUALITIES
In this section, we introduce some basic inequalities. Fact A.5. For any y and 6 (0, 1), we have
(xr + y)² 2 (1 -1- ()=6)a² t
Proof. It suffices show
2²2 + 2ry + y² 2 (1 c)x² - -22
Re-organizing the above terms, we have
ex² + 2xy (1 + E 1) 20
Thus suffices to show that
ex
2 + 2xy + E = 20. 0.
is obvious that
xx² + 11 3² 2 2/zy|-
Thus, complete the proof.
Fact A.6. Let n be an arbitrary positive integer. Let 0g 2 0 and b; 2 0 for all E [n]. Then, the following two inequalities hold
min {a,} sE(n)
min(a,)
b,
a&b I max(a;
Σ
b
kE(n)
iE(n]
iE(n]
iE(n)
min {b,} iE[n]
Σ
a
a&b S max
Σ a+.
vE(n
iE|n]
iE|n]
(E
A.3 PROPERTIES OF NORMS
We state some standard facts about norms without providing proofs.
Fact A.7. We have the following facts about norms:
Part 1. For any matrix A E Rax", let A, denote the j-th column of A. Then we have
Σ 1, = ||||
Part 2. For any psd matrix A, for any vector x, T Ax 2 mmin
Part 3. Let U E R denote an orthonormal basis. Then for any lo x & matrix B, we have || |B]]
Part 4. For any matrix Bonk we have ||||| N |||||/ I A
Part 5. For any matrix and B, (A) 2 min ( - || - B||
Part For any matrix A g and any orthonormal basis Q E R Omin (A)
muin ( (AQ).
Part 7. For any vector E I and for any orthornomal basis Q E Rkxk we have | |2 = ||Qx|22-
A.4 GENERALIZED MATRIX INCOMERENCE
In this section, provide generalized notion of matrix incoherence, denoted by p. Definition A.8. Let gaxk, The generalized incoherence of A is denoted p(A), i.c..
p(4) = 7 max ( k iE|n]
(3)
Claim A.9. When A E x has orthonormal columns, 1 S p(A) <
Proof. Since A is an orthogonal matrix, ||, < 1 for all i E [n], and thus p(A) < x. In addition,
E |t. 12 kt we have max (1l A,, [1]) 2 and then p(A) 2
19

---

Published as a conference paper at ICLR 2025
A.5 ANGLES AND DISTANCES BETWEEN SUBSPACES
An important metric we use in this paper to quantify the progress of our algorithm is the distance between subspaces. We illustrate these definitions below.
Definition A.10. Let X, be n x k matrices with orthonormal columns, i.e., X X 4x and y Ig-
We define tan (Y, X) to be equal to
I X(Y (x)xxI
-I
We define cos @(Y, X) to be equal to
uin XYT)):
we define sin (Y, X) be equal
(T
Let Og be a set containing all k x k orthogonal matrices. We define dist(Y, X) to equal to
||&Q - X||- QEO,
min
Note that by their definitions, we can get
cos (YY(((Y)))))) 1/11(Y
cos8(Y,X) $1,,
sin sin ||LY] Y X||, and
sin (i,, 1
Lemma A.11 (Structural lemma for orthonormal columns, Lemma A.5 of Gu et al. (2024b)). We ler X and Y be arbitrary matrices in B and both are orthogonal. Then, we can get
Lemma A.12 (Lemma A.7 of Gu et al. (2024b)). We let X and Y be two matrices in gaxk and both have orthonormal columns. Then, we have
sin 0(Y,
tan ((Y COS 8(Y,X)
Lemma A.13 (Lemma A.8 of Gu et al. (2024b)). Let X, Y B be orthogonal matrices. Then, we can get
sin 0(Y, + cos² 0(Y 1.
Lemma A.14 (Lemma A.9 of Gu et al. (2024b)). Let X and V be two matrices in Bank with orthonormal columns, then, we can get
tan XX, 2 sin O(Y, X)
( ta > cos (Y.X)
dist(Y, 2 sin 0(Y, X)
sin + 1-cos O(Y.)) cos
> dist(Y, X)
tan (Y,, 2 dist(YY X)
20

---

Published as a conference paper at ICLR 2025
B MORE RELATED WORK
Sketching To achieve the crucial speedup, we utilize sketching- preconditioners and we therefore provide overview of the sketching literature. Roughly speaking, given a tall dense matrix A, the goal of sketching is to design family of random matrices II such that, if we randomly sample II, we have
S
S has much smaller number of rows than A (thus the matrix SAA is close to square matrix rather than rectangular);
A preserves singular values of A with high probability;
S can be quickly applied to A.
Given such a family 11, it is natural to apply an to A then solve the smaller problem directly. This is the so-called sketch-and-solve paradigm. Sketch-and-solve has led to the development of fast algorithms for many problems, such as linear regression (Clarkson & Woodrufff 2013; Nelson & Nguyên, 2013; Song et al., 2023f,d), linear and kernel SVMs (Gu et al., 2025). low rank approximation with Frobenious norm (Clarkson & Woodruff, 2013; Nelson & Nguyên, 2013), matrix CUR decomposition (Boutsidis & Woodruff, 2014; Song al., 2017: 2019c), weighted low rank approximation (Razenshteyn et al., 2016), entrywise & norm low rank approximation (Song et al., 2017: 2019b), tensor regression (Song et al., 2021a; Reddy et al., 2022; Diao et al., 2018: 2019), tensor low rank approximation (Song et al., 2019c), tensor power method (Deng et al., 2023b), and general norm column subset selection (Song al., 2019a).
As modern machine learning centers around algorithms that are iterative in nature. Sketching can also be adapted to an iterative process reduce the cost of iteration. This the so-called Iterate- and-sketch approach and it has led to fast algorithms for many fundamental problems, such as linear programming (Cohen et al., 2021; Song Yu, 2021; Jiang et al., 2021), empirical risk minimization (Lee et al., 2019; Qin et al., 2023), semi-definite programming (Gu & Song, 2022; Song et al., 2023e), John Ellipsoid computation (Song et al., 2022c), Frank- Wolfe algorithm (Xu et al., 2021; 1; Song et al., 2022a), hamming estimation (Hu et al., 2024), reinforcement learning (Shrivastava al., 2023), k means clustering (Liang et al., 2022), online weighted matching problem (Song et al., 2025), barrier functions (Gu et al., 2024a), noftmax-inssirres regression (Deng et al., 2023a; Gao et al., 2025; Li et al., 2023b; Sinha et al., 2023; Li et al., 2023a; Song et al., 2023a; 2024), leverage score inspired regression (Li et al., 2024), federated learning (Song al., 2023b; Bian al., 2023), discrepancy problem (Deng et al., 2022; Song et al., 2022b), non-convex optimization (Song et al., 2021b;e; Alman et al., 2023; Zhang, 2022), and attention approximation (Gao al., 2023a:b).
(Weighted) low rank approximation Low rank approximation has emerged as crucial technique in machine learning and numerical linear algebra, enabling the extraction of essential structures from high-dimensional data while reducing computational costs. The goal is to find X, Y E B which minimizes IM XY IF. has been applied to numerous fields, including training deep neural networks (Song et al., 2021c), approximating attention mechanisms (Alman & Song, 2023; 2024; Chen et al., 2024), maintaining dynamic Kronecker products (Song et al., 2023c), and tensor product regression (Reddy al., 2022). In many practical scenarios, certain entries of M hold greater significance than others, giving rise to weighted low-rank approximation, where the objective is to minimize || (M XYT)Il for some weight matrix W R (Li et al., 2016; Razenshteyn et al., 2016; Gu al., 2024b; Liang al., 2024).
C WEIGHTED MULTIPLE RESPONSE REGRESSION SOLVERS
In this section, we show how to solve weighted multiple response regression by solving standard linear regressions. We present randomized and fast regression solvers based on sketching and preconditioning.
C.1 GENERIC REDUCTION AND ERROR CONVERSION
In this section, we present a generic framework to reduce the weighted multiple response regression problem to solving (nn) ordinary least-square regressions. This simple and efficient reduction enables
21

---

Published as a conference paper at ICLR 2025
us to deploy fast regression solvers to handle the approximate updates in alternating minimization. We also present a tool that converts the relative error on the regression cost to the quality of approximate solution.
The first lemma states that the cost of a weighted multiple response regression can be decomposed into summation of n weighted linear regressions.
Claim C.1. Given matrices M, W _ and X, Y E gu we have
XY min |M COM-XYTI0 li = Σ min | DDwwYYe, - D Dww, M..112 XER**
X2 Elct
and
min YERuXi
|M - X
Σ
min ER
| 1D XY,: - Dw: M.
Proof. Since the two equations can be proved a similar way, we only prove the first one.
min
||
XY T l
min XERoxk
WXr -M))
XERext
min XER*XA
Σ I 1Dvw, (VX;,; Mi.))
min XER
Σ ID WDW,,Y YX D M...
Σ
ID, D Dww, M..Il XR,
min
where the 1st step due ||Al definition, the 2nd step is by rewriting each row as an independent regression problem, the 3rd step follows from simple algebra, and the last step follows from the fact that there no X in minxera** Σ.= D| DD/.WY X,: Mr. l12 but only X Thus, we complete the proof.
The next lemma provides a simple conversion of weighted linear regression to ordinary least-squares, via a scaling trick.
Lemma C.2 (Lemma B.6 of Gu al. (2024b)). Let A be real n d matrix with 2 d, b be an n n-dimensional real vector and w be a non-negative n-dimensional vector (weight). Let co E (0. 0.1) be accuracy parameter and 8, E (0, 0.1) controls failure probability. Suppose that T T(n, d, co. 86) the running time of regression solver, and x' E Rdd is the outpul of the regression solver safisfying
A bll2 S (1 + Eo) min || - bll2
ZERd
with probability at least 1 1-80 80.
Then, there exists an algorithm whose running time is O(niz(A)) + [m, d, co, 6g)
and outputs a vector x ERd,, which satisfy
A| bllw S co)
min zERd
||
bll..
with probability at least 1 - 50.
One of the main reasons Li et al. (2016) resorts to exact weighted multiple response regression is that most approximate solvers provide backward error guarantees on the cost of regression. On the other hand, we would like the approximate solution of the regression to be close to the exact solution. The following lemma converts the backward error the cost, the forward emror on the solution.
22

---

Published as a conference paper at ICLR 2025
Lemma C.3 (Backward error, Lemma B.5 in Gu et al. (2024b)). Let A be real n X d matrix with n 2 d, be an n--dimensional real vector. Let ZOPT be the exact solution to the regression problem
min | bl12-
Suppose that there exists a vector x' E RRa, satisfying
|| bl[2 S (1 ++)
min wER
| b112-
Then, we have
1
| - ZOPT|12 S ((ve)
-|AAOPT - bbll-
mmin
J (A)
Before wrapping up this section, we present meta algorithm for solving weighted multiple response regression.
Algorithm 3 Fast, high precision solver for weighted multiple response regression
1: procedure MULTIPLERIGRSSSIONN. E R"x*, B E Ba*k W E R"x")
2:
A; is the i-th column of A
3:
b W, is the i-th column of W
4:
Dw. is a diagonal matrix where we put W; on diagonal, other locations are zero X; <- min [Dw, Br Dw, Aill2
5:
6: return X
D
X E R
7:
end procedure
8:
9: procedure FASTMUULTPPRREEGRSSSSONAA. Raxn B Runk W (Rxx)
10:
>> A; is the i-th column of A
11:
b W; is the i-th column of
Dw. is diagonal matrix where put W, on diagonal, other locations 13: X, < HICHPPRECISOONREODWW.... 6,)
12:
è Algorithm
14: return X
X E R
15:
end procedure
C.2 Low ACCURACY SOLVER
We provide an algorithm that uses sparse sketching matrix to obtain a low accuracy solution (inverse polynomial dependence on accuracy parameter e).
Definition C.4 (OSNAP matrix, (Nelson & Nguyên, 2013)), For every sparsity parameter s, target dimension m, and positive integer d, the OSNAP matrix with sparsity is defined as
Srj
8rg Ordi
V8
for all [m].j E [d]. where and independent Rademacher random variables and 8,3 Bernoulli random variables with
For every [d]. Ere(m) 8 = s, which means each column of contains exactly nonzero entries.
For all r E [m] and E [d], E(o,,l = s/m.
vT [m] [d]. EI (r.te 5,,1 S Irrer E| (s/m)iTI, i.e. 8y, are negatively correlated.
Crucially, the OSNAP matrix produces subspace embedding with nearly linear in d row count.
23

---

Published as a conference paper at ICLR 2025
Lemma C.5 (CCohen, 2016)). Let S Rmx"" be an OSNAP matrix as in Def. C.4.
Let e, (0,1) be parameters.
For any integer d n,
= Oog'dd/g(d//0))
=s=
==0de=log(d//)) log(d/8)).
then an s-sparse OSNAP matrix S is (c oblivious subspace embedding, i.e., for any fixed orthonormal basis U E Igad with probability least 6, 6, and the singular values of SÚ lie in [1 +c
To distinguish with c, 5 for our final algorithm, here we use 6g for the subroutine (approximate linear regression).
Lemma C.6 (Input sparsity and low accuracy regression). Given matrix A E Raxd and vector b E RR", let Eo E 0.10 and óu E (0, 0.1). there exists an algorithm that takes time
O((E7) (nz((A) + eu 2 log(d/66))
and outputs x' E IRd such that
|Ax bll2 S (1 cD)
min zERd
|| - bl12
holds with probability 1 6g-
Proof. To obtain desired accuracy and probability guarantee, we pick S to be OSNAP (Defini- tion C.4) with
m (dlo((d))()
and
O() log(d/66)).
We simply apply S to then solve the sketched regression min ||SA Sb|l2.
As is a matrix where each column only has nonzero entries, the time to compute SA is
( (n)( O(55' nnz log(d/6))).
The regression can then be solved via normal equation, i.e.,
(AT ST A STb.
The time to form the Gram matrix is
Ommd))
computing the d x d inversion takes ((d²) time, and forming the final solution takes another
(mdd)
time. Overall, this gives runtime of
O 2d log(d/6o)).
Thus, the overall runtime is
Ol(5' mnz(A) 202) log(d/50)).
n
24

---

Published as a conference paper at ICLR 2025
C.3 HIGH ACCURACY SOLVER
Our key algorithmic ingredient a high accuracy, iterative and sketching- solver for regression. The sketching matrix will be using is the dense subsampled randomized Hadamard transform SRHT due to loss of structure in the iterative process.
Definition C.7 (Subsampled randomized Hadamard transform (SRHT), (Lu al., 2013)). Let P be a random sampling matrix in {0, 1 1 and for each row of P, there exists a 1 at uniformly random position.
Let H {-1,1)mxm be the Hadamard matrix.
Let E IR"*" be a diagonal matrix, whose diagonal entries are all {-1, +1} with the same probability.
We define the SRHT matrix S R x x as
S = PHD.
mm
Remark C.8. For real matrix A it takes Ond n) time to apply S to A.
Lemma C.9. Let S E Rmxx be an SRHT matrix (see Definition C.7), c, 8 (0,1) be parameters. Let d be an arbitrary integer, which is less than or equal to Suppose U E Raxd be a fixed orthonormal basis.
dd
(n/o)). (on=)))) Let
We say that S is an (c. -oblivioos subspace embedding if the singular values of the matrix SU are in the interval [1 c, + e), with probability least 1 8.
Lemma C.10 (Dense and high accuracy regression, Lemma B.1 in Gu et al. (2024b)). Let A be a real n x matrix, be a real n-dimensional vector, e 0.1) be an accuracy parameter and 5 E (0, 0.1) be the failure probability. Then, there exists an algorithm that takes time
O(Idlogn + de logg(n//0)) log(1 /c))
and outputs a vector x' IRd satisfying
|| - b122 S (1 E) min | - bl|2
wERd
with probability at least 1 -
D KEY PROPERTY FOR ROBUST UPDATE
In this section, we prove crucial properties of the algorithm that enable the approximate updates. In Section D.I, formally define several necessary notations needed to analyze our robust updated step. In Section D.2, we analyze the key properties of our robust update step.
D.1 DEFINITIONS FOR UPDATE STEP
We present closed-form solution for linear regression via normal equation.
Fact D.1. Define := arg min, || - bl]2, then we have
x = (A A) ATb.
Similarly, for weighted regression, we define
= arg min Dw, Ax Dw.bll2-
Then, have
= (AT Dw, A) AT Dw:
Definition D.2. We define ç as
{ captures the maximum incoherence of the ground truth.
25

---

Published as a conference paper at ICLR 2025
Definition D.3. We define to be parameter that distinguish random and SVD initialization.
For random initialization, we set := pkk
For SVD initialization, set 1j
We next define the value choice of 7, which controls how far away the weight matrix W be from the all-1's matrix.
Definition D.4. Let
1
7 S 100
poly(k, T, n0
where eg fixed constant between (0, 1/2]-
For the convenience of analysis, we define the following threshold parameters: Definition D.5. Let ž 10 denote sufficiently large constant. We define
Aa Co 1.5 " X2I (o.,nn) + Ca nukk? ,
A Ca Ink:
The choice of Ad is decided in Eq. (9) and Eq. (15). The choice of Ay is decided in Eq. (16). The next two definitions capture the error gap of our algorithm. Definition D.6. Let Aa and A be defined as Definition D.5. We define
A == - dist V) + A |[ NJI.
and
Ag = 0.01 dist (ooisttY,, V) + 0.01A, || N|| + 2
By properly controlling the error 6sk, we can show that Ag is constant factor smaller than A.- Claim D.7. If the following condition holds
Esk S 10 A03 |IW NII?
then have
Ag S.11A.
Proof. We have
A, 0.01 dist(Y, V) + 0.01A, ||W o N +2
< 2.00 (YistYY,V) V) 0.01 ||W NI| 0.02A, WIW N
< 0.1A
where the second step follows from condition on 6sk, and the last step follows from the definition of A...
By setting the error and failure probability appropriately, show the extra blowups in our algorithm are of the order poly log n.
Claim D.8. By the choice of este we have
log(n/csk) (log(n))
By choice of failure probability (6o) of sketch,
log(n/60 /0)) log(n log (1/e))
26

---

Published as a conference paper at ICLR 2025
Proof. By assumption W, we can see that
|| S poly(n).
Since W N is noisy part, it is also natural to consider
|[ Mi poly(n),
otherwise, it is not interesting.
Since the noise cannot be 0, thus it is natural to assume that
| 2 1/ poly(n).
Thus, we know
1/ 1/poly(m) S IIW o N S poly(n).
We also know that k S n.
Now, we can compute
log(n/cak) - (1og((n/4))l1) o N))
5 (log((m/lll [/ll NI)?))
- (log(())
where the first step follows from choose 6s O IW N1 the second step follows from A 2 1, the third step follows from | N 2 1/ poly(n) (see Eq, (4)).
Sum over all the T = O(log(1/e)) iterations, so
log(n/66) = O(log (og(1(()))
D.2 KEY LEMMA FOR ROBUST UPDATE STEP
Table 2: For convenience, we provide table to summarize the notations in Lemma D.9.
Notation
Meaning
X
Optimal matrix for exact regression
X
Clipped matrix of X
X
Optimal matrix for sketched regression
X
Clipped matrix of X
Algorithm 4 Clipping rows whose norms are larger than a constant factor of Ç.
1: procedure CLIP(X R=**)
2: nk
3: for 1 to n do X1. 1 4g then
4: if
5:
X,, X
6: else
7:
X, +
8:
end if end for
9:
return X
11:
end procedure
Here, analyze the key properties of the robust update step.
27

---

Published as a conference paper at ICLR 2025
Lemma D.9 (Key lemma for update step). Let Y E R be (column) orthogonal matrix. Let be defined as Definition D.2. Let Au be defined as Definition D.6.
We define matrix X E B as follows:
X
=
min XeRuxk
|M - SYTlw,
We define matrix X E gunn as follows
X
if 4
X
otherwise
We define X E Raxk R Rk to the QR decomposition of X, Le. X = XR.
We define X E Buxk to be the sketch solution such that for all [m]
Xx.I1² & Cak/n.
We define X Raxk denote the clip of the sketched solution. Recall that M UEVT. Then,
Part 1.
- - < A
Part If Du < 0.10mia(MM) then
dist(U, S 8mn/0min (M*)
Part 3. D,, < 0.10min (M*), then
p(X) S 8p/ m (M*)²
Proof. Proof of Part 1. Recall the weighted multiple response regression
min XER
|| | XY llv
The above problem can be written different regression problems. The i-th linear regression has the formulation
|I YX,; D Dvw, X..
min
We have
1. M.. Dw,Y Y Dw: V)-
-(MM), + NNT) Dw, Dw,Y YY Dw:Y)-
= =(M),, Dw, YY Dw,Y)-++N Dw..Y(YDD..))
Dw,Y YY Dw
(5)
where the first step follows from the Fact D.1, and the second step follows from M, = (M)), +.,,,,+N,,. (because M = M + N), and the third step follows from simple algebra.
Given M U.VT, the first term in Eq. (5) can be rewritten follows:
(M*),. Dw,Y YY Dw,Y)-1
= -,,,- ΣV Dw, Y Dw,Y)-1
-U,,. ΣVT(YYT
Dw
U,T. EVTYY Dw. Dw, + yy Dw,Y Dw,Y)-1
U.J ΣVTY -Y+UU -UEVVVIL Dw,Y Dw.Y)
(6)
28

---

Published as a conference paper at ICLR 2025
where the first step follows from the fact that (M"), U,,VV, the second step follows from YYT + VLYL. the third step follows from simple algebra, and the last step follows from AA-1=I.
Combining Eq. (5) and Eq. (6) we have
X__ - -U, = U T D. Dw, + Nw: Dw;
(7)
We define set C [n] as follows
T= {i [n] Tmin (Y Dw, Y) S 0.250/n)
(8)
We upper bound (T) in different ways for SVD initialization and random initialization.
SVD case. We have 7 = 1. Using Lemma F.1 and choose c = (1), we have
| I 10 a -3 "p2k²>² || VV Y
10 a -3 p22k0y pkk (IIW]l,,,1m)) IIV V11
$ 0.143 (V,Y)//E
- /
where the second step follows from & = pk/n, the third step follows from Definition of Aa (Definition D.5).
(In particular, the third step requires
A 2 A(x 215 [[ (n)
(9)
Random case. We have 1) = pu. Using Lemma H.I, with high probability we know that [T'] = 0 S Ag/E.
In the next analysis, we unify the SVD and random proofs into same way.
For each i [T.] have
IXT U. YI2
VUTTETYY Dw, Y NT, DECIYYDDDWYY.
Dw, Dw;Y)) -)-1
II Y.Y Dw,Y Y + Dw: 11 Dw, 12-
s 20a 20a 312 -I.,,VVV YıYī Dw: +W, Y)
< 20a -,, 21I Y Dw: Y112 + ||N Dw: YI
I 40a -1p² (I1TT..113))) V Dw, Y + IIN.T. Dw. Y
< 40a 2122 pk
Σ PIVV
Dw. Y1 VI Dw; Y (1a)
(10)
where the first step follows from Eq. (7), the second step follows from | < ||||| ||>||2 the
third step follows from min (Y Dw: Y 2 0.25a/n / (for all E [n]][T see Eq. (8)). the fourth step
follows from (a + 6)2 < 2a + 262. the fifth step follows from |_Az||a S |||||| || the sixth step
follows from | S pk/n.
Taking the summation over i E [n]]T coordinates (for Eq. (10)), we have
Nkk n
Σ I v
U <
Σ sE|n]yr
40a
1072
Σ
Dw. Y INN Dw: F1
(11)
For the first term in Eq. (11) (ignore coefficients 40a 272 and uek | we have
Σ
IVTY Dw: 211² < k (V) dist v)
VT
29

---

Published as a conference paper at ICLR 2025
S ny² poypuo (E) dist(Y, v)2
(12)
where the first step follows from Lemma F.2, the last step follows from p( S / (E.) For the second term in Eq. (11) (ignore coefficients 40a=2 have
Σ /T
IN Pw: Y12 s
Σ IN iE(m]
II( NYY
S (W N)
(13)
where the last step follows from Fact A.7. Loading Eq. (12) and Eq. (13) into Eq. (11), we have
[ U. r I 40a 312 pk | ny uo min (E) dist(Y, k|l o N))
Σ
40a 312 (x?2) dist ( + 40a k// o 21
(14)
where the last step follows from |ل = 1 and 5 = max O2i)/mmin(2)
mmin
Thus, we have
Σ [ U. y 40a 317 (1222T97)) (Ys 40a k o NII2
0.01 dist(Y, + 0.014} |[ o N||2
0.11
where the first step follows from Eq. (14), the second step follows from Definition D.5, and the last step follows from Definition D.6.
(In particular, the second step above requires
A 2 N( m 7.55))
(15)
and
D5 2 (nk))
(16)
)
By Definition D.6 and choosing 6sk to be sufficiently small Claim D.7, we know that
A S 0.014
Then, show
Σ T
|[T -
U.
Σ
|| - Σ -
*E|M](
iE[M](T
< A yz
Note that
| U.E T 1 < pk/n= E.
If | 4,, then
X - U..VVTY > 2 - V V
which implies that
(17)
30

---

Published as a conference paper at ICLR 2025
We define set S C [n] as follows
fi= | 2 48).
Then we have
|S| = lfd E [n] | 2 4EH1
< E [n]] | - -U,.VV Y 2 E)
/
where the first step follows from the definition of S, the second step follows from Eq. (17), the third step follows from EiemoNV [ 1,: We can show
U. Y| - g'
|
- USV
Σ
U B
Σ IETUS
[ UTT.V
Σ P - U..VV Y
iTTUS
Σ iETUS
|| U. EV 1
iES
Σ
2T1 + |I. + Σ
[ - -T
IETUS
iETUS
IT TUS) S| 2- (45 +5) + +
(TSS 100 + A
< 50Ag
A
where the first step follows from the definition of X LVV T ||2, the second step follows from
S S [n], the third step follows from XT. = XT when é & S, the fourth step follows from the triangle inequality, the fifth step follows from || < 46 and |U]]EVTY]] < E, and the sixth step follows from [T' U SI < /. the last step follows from Claim D.7.
Proof of Part 2. We let B ΣV Y and have
sin 0U,
|IUXXI
-IUI((X-UB]E -UBBR
I1(XX -UBIII ||
1(X UB)Il
minu (X)
D
(18)
miiu (X)
where the first step follows from the definition of
sin XU,
(see Definition A.) the second step follows from
X == (X-
the 3rd step is due to the Cauchy- Schwarz inequality, and the 4th step is because of
|| R ||-1|||
o miin (X)
31

---

Published as a conference paper at ICLR 2025
and the 5th step follows from
I1X Y)|| - A
which proved in part 1, and it infers 11(X - -UB]] - Au for B -V Y.
Using ||xX UB) < Aus we have
Mmin (X)
mmin (UB) Au
min mmin (UVVT)) Au
minn (EVT)) - Au
miin (M*) Cos cos 0(Y, V) A..
min (M (M2 Au
min (M //4,
(19)
where the second step follows from how defined B, the third step follows from U has orthonormal columns, the third step follows from
min (ΣVTY)= 2 min (M" cos O(Y, V), where the forth step follows from cos 2 1/2, and A. S minn (M*) 10. Then, by combining Eq. (18) and Eq. (19), we have
sin (U, mmiu X) Au
s 4Du/ Tninn (MT*)
(20)
Therefore,
dist XU, 2in OU,
S <8 /mn/ (M*).
where the first step follows from Part of Lemma 1.14, and the last step follows from Eq. (20).
Proof Part 3. Given X = XRR we have X X,T. and and
|| 1 S
X 1 11 111
&
min (X)2
(mii (M*) 2
S 8E/ mmiu (M*)²
(21)
where the first step follows from X X R and Cauchy-Schward inequality, the second step follows from | < & and || Tuin (X)' the third step follows from (X)= S
(Tmin N 2 and the last step follows from A. I 0.10mmin (M*). (M
To see this, we have
R 1 Amax (R (R-
Amax (X (X-x)11
Amin (X X
1
2 (X)*
32

---

Published as a conference paper at ICLR 2025
where the first step follows from ||A]] Amax (AAT), )2 the second step follows from X X RTR. the third step follows from Amax ( (A-1) = Amia (A) and the last step follows from X(A A) (A).
Then,
P(X) = max x,:,,:112
Y1
re(n) k
12 max
11,X,,:112
4: iE[n]
88
mmin (M+)2
8tu
Fmiu (M-
where the first step follows from the definition of p(X). the second step follows from simple algebra. the third step follows from Eq. (21), and the last step follows from μ.
E THE ANALYSIS OF THE INDUCTION LEMMA
The goal of this section is to prove induction.
Lemma E.1. Suppose
Au S A dist v)+ + As || |[Wo]]] NiI.
and
Aa S 100
Umin (M*)
For any 2
1
dist(X,,U) < 2
+ 100 (M*)-1 D, |IW o N]],
dist (Y,, V) < 2 + 100 (M0"" A, ||W NN],
(22)
Proof. By using induction, we show that Eq. (22) holds.
Base case: By Lemma F.3, Yı satisfies
dist(Y,, V) s _ + " D, IIW o NJ
Inductive case: Suppose that it holds for the first cases. By definition A., we have
A < Aa disti(Y, V) + + |W o Nl
(23)
We have
Mist(X,-1,.4 U)
A.
m (M*)
8
(Ad dist(Yi, V) + Ay |[ a NI)
mmin (M
$ dist(Y,, V) + 80 (M*) A, || N
where the first step follows from Part 2 of Lemma D.9, the second step follows from Eq. (23) the last step follows from A (M*) 1/100.
33

---

Published as a conference paper at ICLR 2025
F TOOLS FROM PREVIOUS WORK
In this section, we state several tools from previous work. In Section F.I, we introduce varied version of Lemma from Li al. (2016) bound the eigenvalues. In Section F.2, we bound (IT)) ( Dw,Y Y112. In Section F.3, we summarize the base case lemma.
F.1 BOUNDING EIGENVALUES
Now, in this section, we start to bound the eigenvalues.
Lemma F.1 (A variation of Lemma 10 in Li et al. (2016)). Let Y be a (columm) orthogonal matrix in Baxk Let c E (0, 1). We have
Ki [n] min (Y
10a
IV Y
da
Proof. Let j be an arbitrary integer in [n]. Let be greater than 0. j is called "good" if
|| Villž s g².
We define S C [n]] as follows
Sg bj== v3113 s g²}}
For convenience, we define Sg C [n] as follows
S :- Sg.
We choose y to satisfy the following condition
2a
g² 20 Wl
(24)
Let be an arbitrary unit vector in RR., Thus, we have
a TT Dw,Ya= =
Σ (Dw.)
sE(m]
Σ
(Dw.) (a,
JES,
Σ (Dw,) (ow))(((,,, Vz) + V v2))2
JES,
2 (1 e/4)
Σ (ww:), (a, v2)² - -2 Σ
(Dw.)slayY, -V2)
jES,
jES,
2 c/4) Σ
(Dw.); (a, V2)² -4--19 (Dw.);
jES,
jES,
2 (1 e/4) Σ
(Dw,)/(a v2)² 11g Σ (Dw.),
jES,
jE(n]
2 (1 e/4) Σ
(Dw.), (V,)) Σ
(Dw.), vy)² - Σ (Dw.)s
jE(n)
1E3,
jE(n)
uk
2 c/4) Σ
(aw.), v2)²
Σ (Dw.), - Σ (Dw.)ss (25) JE(n)
jC(m)
JEB,
where the first step follows from simple algebra, the second step follows from S [n]. the third step follows from the property of the inner product, the fourth step follows from Fact A.5, the fifth step follows from the definition of Sz. the sixth step follows from (Dw:); 0, the
34

---

Published as a conference paper at ICLR 2025
seventh step follows from 1- c/4 I 1, and the last step follows from the property of V (e.g.
(a, Vj)² < lalß 1v x IV 1 & pk/n).
We show that
Σ (DWw))(a, =1) =V"VT(Ww:) jE
T (V (Dw.j,V
a.
(26)
where the second step follows from Fact A.7 and the third step follows from definition of a (see Definition 3) and minn (A) & miin (B) if A 3 B.
Moreover, recall
|
max
Σ
1(w.)sb
(27)
vE
je(n]
We can show that
4E 4
(Dw.), S $: 46 9
jE(n)
6²a
46
| 10
20 W
o0,1
Ea
(28)
where the first step follows from the definition of || W from the definition of g2 (see Eq. (24)), and the last step follows from simple algebra.
20.1 (see Eq. (27)), the second step follows
We define
T : fie e [n] mmin (WDWYY (1 - e)a}
Let us consider
(Dw.)..
JES,
We define
S = (i E [n]
pk Σ
(Dwi)i
2 6a
4 }.
(29)
JEB,
If i $ S, then we have
1k (Dw:); -19² Σ
a Ty Dw: Ya (1 (1-4/4) Σ (Dw,) 12)²
(Dw.)a
jE(n)
JEB,
JE(m]
uk
2 (1 e/4)a
(Dw.); - -1 9
(Dw:)s
jEZ,
jC(m)
2 (1 e/4)a ea/4 4E 111 122 (Dw,),
jc(m]
2 (1 c/4)a ca/4 /4- - ca/4
2 (1-c)a.
where the first step follows from Eq. (25), the second step follows Eq. (26), the third step follows from the Definition of S (see Eq. (29)), the fourth step follows from Eq. (28), and the last step follows from simple algebra.
35

---

Published as a conference paper at ICLR 2025
In summary, know that if S, then E T. By taking its contraposition, have that if T, then E S.
Thus, can show that
TT| S |S]
In the next few paragraphs, will explain how upper bound [S]. Using Fact A.7
IV - V YI JE
By simple counting argument (you have positive values, their summation is JV Yl you can't /l1//9" of them that are bigger than g²), we have
have more than |V
5 S I 9l
(30)
Let us I be the indicator vector of S, i.e.
i S; otherwise & S.
Vi [n]. (us)
0
Let Ugg E IR" be the indicator vector of 3g i.e.
{1
if E ggi otherwise $ Sg
Vi [n], (ug)s
Then, know that
us Wug = Σ Σ (Dw.)s
iES JEB,
|S|
min ES
Σ (Dw.)s
JE
can S
2
(31)
4pk
where the first step follows from simple algebra, and the second step follows from simple algebra, the third step follows from Definition (29).
On the other hand,
Ug W Ug = ug 1, 1 Ug + ug (W 1,, )m,
||us lla || + u (W 1, 1, )mg
||ug 1 ||0/]] 1 + ||uslle I|W 1, 17 |
< ||u
[||]
1|,,]11 +m |lus||2 -1/u,112
< [5| + ym |||||
(32)
where the first step follows from simple algebra, the second step follows from simple algebra, the third step follows from x Ay |><||| |||llll the fourth step follows from Definition 2, and the last step follows from ug and Mg are indicator vectors.
By combining Eq. (31) and Eq. (32), have
B + nn
B3,//S||
>
60m
(33)
4pk
Note that if A + B 2 C. Then if A C/2, then B C/2.
36

---

Published as a conference paper at ICLR 2025
For the terms in Eq. (33), define A, B and C' follows
A 1
B = 7
C ccm
4pk
If
15 < can
8yck
then, have
mm V /3/1/1| 2
60m 8yk
The above equation implies
|S| S 500
[5g
e²a
p²
< 500
(IIV - //²
2a2
1 k²
20 1
500
||V - YI
(34)
2a²
62cx
where the second step follows from Eq. (30) and last step follows from Eq. (24).
We have
lfi E [n]
Omin
(TTDw, D) < (1 (aatlsls
1222?²
I 10
11 - Y ||
$ 10
||
k. VVV Y|I
e'a
where the first step follows from the definition of and the second step follows from combining Eq. (34) and Eq. (24), and the last step follows from Fact A.7. This completes the proof.
F.2 BOUNDING (VLY) Dw,
In this section, we bound [VTY Dw: Y by multiplicative factor times |Y VII Lemma F.2 (A variation of Lemma 11 in Li et al. (2016)). Let Y be (column) orthogonal matrix in Raxk Let i E [n]. Then we have
Σ II s 2 VI2
iE(n)
For the completeness, we still provide the proof.
Proof. Let j'.j be two positive integers in [k]. Y, represents the matrix Y's j-th column. V represents the matrix Y YLYTV's j-th column. We define wid' Rd as
= (V))
We need to show that the spectral norm of
VIY.T.D..Y.
37

---

Published as a conference paper at ICLR 2025
is bounded. Note that (V,, Yy) 0, which implies that
Σ <E(n)
x 0.
For V Dw,Ys,
JJYYXIDw: Dw;Yy
E (V),)s(()))(())) mE(m)
(Dw, s
(35)
(E(n]
where the first step follows from the definition of
V,TY,Y, Dw.Ys
and the second step follows from (V)))(Yy)) = x/.". It implies that
Σ
Dww. (ww:)add') |1 113
*E[n] =E[m]
II( 1, 1 Jart
S I|W 1,, 1,, 1|
S $2
where the first step follows from the definition of ||Wa:'''' 113. the second step follows from 1, 20 0, the third step follows from ||z||2 ||A||||= 2, and the last step follows from ||W 1, S yn (see Definition 2).
Observe that
|| 112 =
(aridd
sE|m]
Σ (V)))
mE/m
p(Y)
(VA)
(E(M)
(YY) 11V,12
p(Y)k
| VII
p(Y)
(Y - -V)1 [
yk
||Y ||,,YTI -1|-111 VII
JY -I2
where the first step follows from the definition of |للا5' 113. the second step follows from (a'hd)² (YG)) the third step follows from definition of p. the fourth step follows from the definition of || 112, the fifth step follows from the fact that V, is defined to be the j-th column of VLYLV. V, the sixth step follows from YT y = 0, the seventh step follows from ||A||| < ||A]| [[B]] and the last step follows from |Y 1.
38

---

Published as a conference paper at ICLR 2025
It implies that
Σ VE(N) (E(n]
Dww)saddy S p()))ll VII?.
(36)
Now we are ready to bound V Y Dw, Y. Note that
Y. Dw, YII YV Dw: YI
Σ Y
Y Y,)²
Σ
(37)
J
sE(n)
where the first step follows from ||A|] < || for all matrix A, the second step follows from definition of || [|r, and the third step follows from Eq. (35).
This implies that
Σ iE|m
II Dw:
Σ
Σ
Σ
(Dw;) (ww))agd''
iE[n] j.k' sEln
Σ momYlll VII
EEE
E
E 72 m VII
SE(R
kky? p( (nki)llly - V112
7 ((kk) V
where the first step follows from Eq. (37), the second step follows from Eq. (36). the third step follows from simple algebra, the fourth step follows from the property of Σ (e.g. E a = an), and the last step follows from simple algebra. This completes the proof.
F.3 SUMMARY OF BASE CASE LEMMA
We state general base case lemma that covers both random initialization and SVD initialization.
Lemma F.3 (General base case lemma). Let A 20a 20a-ank For the base case, we have
dist(Y V) < 80 (M*) -1. A IW o N|]F.
The proofs are delayed into Section H and Section G in which we analyze the random and SVD initializations with different parameters.
G SVD INITIALIZATION AND MAIN RESULT
In Section G. we introduce our assumption on In Section G.2, we bound II(W - 1, o HJJ In Section G.3, we analyze the property of the rank-k SVD. Section G.4, analyze the properties of dist and In Section G.5, we present our main result.
G. ASSUMPTION
Here, we set the parameter 8. Definition G.I. We assume that
= 0.001
|| N| < ax MaminnM? /M*
39

---

Published as a conference paper at ICLR 2025
G.2 BOUNDING I(W 11,1, H
In this section, we bound I(W - 1, 1 o H|] in terms 7, the rank & and the top singular value d Lemma G.2 (Spectral lemma, Lemma in Li et al. (2016)). Let K and J be (column) orthogonal matrices, whose sizes are n. x n. Let H be an arbitrary matrix in I satisfying
H = AEB
(38)
where E Rank B R and E E Rixk Note that and B might not be orthogonal, but Σ is diagonal. The matrix W E R=x is entry wise non-negative matrix, which has artificial spectral gap, satisfying
W 1,1, + +nnJEw KT,
where
Ew =1
Let 01 == max or (E). Then, we have
I(W - 1, H|| S Tko1 vP(A)p(B).
Proof. For each i [ni] let A,, denote the i-th row of matrix A E Bank For each [kt], let A, denote the r-th column of matrix A E Baxk Then we have
1
Σ
|| A, o x11
r=1
Σ 11, 12
# ((A)
k ((A)
(39)
n
where the first step follows from the definition of 11 the second step follows from simple algebra,
the third step follows from the definition of 13 the fourth step follows from the definition of (Definition A.8), and the last step follows from Σラ 2 <1 1.
Let 2, y E 12" be two arbitrary unit vectors. Then, we have
xT ((W 1, o H)y x [(W 1,,1,, o (AΣBT)) )
Σ 0, T(W - 1,,1,) A,B,B)
(, x) JEW KJ (B, o y)
y
Σ 1, x1|2 ||EW KT || [B, o 2112
ym
S 7m
Σ
a A, 0 x/ f||2 | o w/12
< \no
Σ
A o x 2|2 [B, o 2lll
r=1
40

---

Published as a conference paper at ICLR 2025
|| o x113) x11))
Σ
|| B, o y12)
r=1
r=1
S sonrs V n (() n =
200 k. (0 (o(A)
where the first step follows from the definition of H (see Eq. (38)), the second step follows from E a diagonal matrix, the third step follows from the definition of W. the fourth step follows from ||Ax||2 S ||A]|| || the fifth step follows from [[ K || ||Ew || ||W S 1, the sixth step follows from 0 max [i] Ty, the seventh step follows from Cauchy-Schre -Schwrr- inequality, the eighth step follows from Eq. (39), and the last step follows from simple algebra.
The lemma follows from the definition the operator norm.
Lemma G.3 (Wedin's Theorem, Lemma 6 in Li al. (2016)). M a matrix, and o1... On are the singular values of M M is a matrix, and oi,..., on are the singular values of M. Suppose that X,Y and U, V are the first singular vectors (left and right) of M, M respectively. If there exists a which is greater than 0 and satisfies
max rE(K+1,-
o S min
m}
(e(1, .}]
then
M M
max ((, 0(, X))
a
G.3 PROPERTY OF RANK- SVD
Now, we first define the parameter A and then analyze the properties of rank-k SVD. Definition G.4. We define A as follows
A 10(yuk + +5)
muin (M*)
Lemma G.5 (Lemma in Li et al. (2016)). Assume that W and M satisfy every assumption. We define (X, Σ, Y) == rank-k SVD(W o M). Let A1 be defined as Definition G.4 and assume that A 0.01. Then, we have
max(tan U), tan O(Y, V)} 5...A,
Proof. We know that
I M M
W (M + +N - MI
W o M -M"I N
IW o M +
S yukom (M") + 5
ypuk + 8
(40)
where the first step follows from the definition of M, the second step follows from triangle inequality and the third step follows from Definition G.I, the fourth step follows from Lemma G.2, the last step follows from o mmax (M*) = 1.
Therefore.
max (r o MM re[k+1,n]
max (W M - -M") + max (r(M)) rE(k+1,n]
rE(k+11.n
< max (W o M - -M*+++ +0
re[k+1,m)
< I[ o oM- M*
yuk +8
41

---

Published as a conference paper at ICLR 2025
0 min (M' 0
0 min ( [M
2
where the first step follows from triangle inequality, the second step follows from the fact that M has rank-k, the third step follows from Eq. (40), the fourth step follows from yuk < 0.1 (M*).
and the last step follows from ô - 0.10 (M
Now, by Wedin's theorem (see Lemma G.3) with
1
2 min (M*),
for
(X, rank- M),
we have
2 + 8) min (M*
(ax(sin@(U,X)) sin YV) v
By our choice of parameters
sin 0 1/2.
Using Lemma A.14, we have
tan 0 S sin 0.
Then the lemma follows.
G.4 INITIAL PROPERTIES OF DISTANCE AND P
We analyze the properties of distance and during initialization. We show that as long as A is chosen properly, the distance and p can be bounded.
Lemma G.6 ((SVD initialization, variation of Lemma 8 in Li et al. (2016)). Assume that W and M satisfy every assumption. Let A be defined as Defintion G.4. Assume that A S 0.01/k. kk. Then, have
.
Part 1. dist(V.Yi) < 1/2.
Part 2. Let p(-) be defined as Definition A.8. We have p(Yi) I 4u,
Proof. Proof of Part 1. First, we consider Vi E BRnk By Lemma G.5, we get that dist (Yst(yi, V) A
which means that there exists Q Olxck such that
|| VII s Az.
(41)
Hence,
Y.Q - V S Vk. || - VIl
< Vi
A
< 10
(42)
where the first step follows from Fact A.7, the second step follows from Eq. (41), and the last step
follows from A < 0.01/k.
Next, consider Vi R In the clipping step, there cases for all i [m).
42

---

Published as a conference paper at ICLR 2025
Case 1. ||1 2 &. We know
|| = ||V,,i112
_
24k
(43)
where the first step follows from Fact A.7, the second step follows from our Case assumption, and the last step follows from the definition of & (see Definition D.2).
have
| - Villa 2 Q12 1
2ph
2
12
2pt
pk
n
pk
(44)
where the first step follows from the triangle inequality, the second step follows from Eq. (43), the third step follows from the property of V,, and the last step follows from simple algebra.
By the definition of clipping, in this time, Y1, O.
11, Vil2 = 100 Vill2
||Vll2
puk
S I1Q Y1.
YY. Vill2
where the first step follows from Y1,; = 0 and the second step follows from simple algebra, the third step follows from the property of V., and the last step follows from Eq. (44).
Case 2. | E. In this case, know
Y1,= = Y
Thus,
|Y i Villa 1| Vill2-
Combining Case and Case 2, we know that for all i E [n].
|| Vill2 < ||Ý Vill2-
Taking the summation of squares, we have
|| - VI < Ý.Q VI
(45)
100
the last step follows from Eq. (42).
Eventually, we would like to show V is close to Y,. Suppose
Y1 -
where is an upper--rianguun matrix.
Then,
sin (V) = |
||
43

---

Published as a conference paper at ICLR 2025
I({ VQ- ')R"
1 VII ||
1 Min (Y 1)
1410 - VII
- V - V uuin (Y1)
where the first step follows from definition of sin, the second step follows from Fact A.7, the third step follows from definition of Y1, and the fourth step follows from ||A]]] S ||||| |B]] the fifth step follows from singular values of R and those of Y, are identical, and the last step follows from
||- S || ||
Note that
Omin (1)
mmin (1QQ)
IV IV I 7 min (V) |/ VII
min (V) || -VILF
min (V) 1/10
1
(46)
where the first step follows from Fact A.7, the second step follows from Fact A.7, the third step follows from ||- || S || I|F, the fourth step follows from Eq. (45), and the last step follows from
amin (V) ===1 1.
Thus.
sin (VV 2 || - V|IF
1
<
10
(47)
2'
where the first step follows from Eq. (46), and the second step follows from Eq. (45), the last step follows from simple algebra.
Therefore,
dist(V,Y)) - 2 tan ((V, Y1) 4in (V, Vi) < 1/2,
where the first step follows from Lemma A.14, the second step follows from Lemma 14, the third step follows from Eq. (47).
Proof of Part 2. For p(Y)), we observe that
Y1, Y
We have
12, Y 112
R '1
{ || '11
E amin Y 1) -1
$ E .2
(48)
where the first step follows from |Ax|| 2 N |A|| ||a||2 the second step follows from || V,,ill 122 S E.
and the third step follows from | = amin Y and last step follows from amin (Y)) - <
(see Eq. (46)).
44

---

Published as a conference paper at ICLR 2025
Note that
((Yi,)) E. ||,,1|2
max
ie(n
2
k 4p
where the first step follows from Definition A.8, the second step follows from Eq. (48), last step follows from Ç = =uk/n. n. This leads to the bound, which completes the proof.
G.5
MAIN RESULT
Finally, in this section, present main results.
Table 3: Summary of our results.
References
Init
Time
(Li et al., 2016)
Random
O([I 2 + nk3) log(1/c))
Theorem H.2
Random
[[[[] + nk²) log(1/e))
(Li et al., 2016)
SVD
O(n²) + [
nk (l)))
Theorem G.7
SVDD
O(ne) + ( W ll nk log (/c))
Theorem G.7 (Main result, SVD initialization). Let rj There is algorithm starts from SVD initialization runs tog(1/c) iterations and generates M, which is a matrix in Roxn and
| M S O(a n)k)) (aaaaaaaT)iil
N N| +e
O Ö((IIW + nk log(1 (og(1/c))
the total running time.
Proof. It follows directly from Lemma J.1 and Lemma J.2.
H RANDOM INITIALIZATION
In Section H.1, we present our random initialization algorithm (see Algorithm 5) and analyze its properties. In Section H.2, summarize main result.
H.1 INITIALIZATION
Now, we start to present Algorithm 5.
Algorithm 5 Random Initialization
procedure RANDOMINIT(P.)
2: Let Y E gank generated with Yi, + Vn by.g, where bi,j is drawn uniformly from {-1, +1}
3: return Y
4: end procedure
The following lemma shows that the minimum singular value of the matrix Y Dw; Y can be lower bounded with high probability over all E [n].
Lemma H.1 (Random initialization, Lemma 9 in Li et al. (2016)). Let Y E B be random matrix with Y., = A bi.g. where bi are independent and uniform variables from {-1, 1}. Let
1. Let k > 1. Let O. We assume that || Fnplog Then, have
Pr [mmnnn (Y
Dw, 2 Apk
Vi E [n] 2 1 11/n?
45

---

Published as a conference paper at ICLR 2025
Proof. Notice that
Y Dw,Y= [Dr))(Dw:)3
jE(m]
For every j E [r]],
(Dw:)y
independent, and
(Dw:)), ! (Dw.)s
Using linearity of expectation, we have
(Z))
(Dw))s
JE(n]
JEMN
Then, we can get that the following equation holds. We use this first and will prove it from Eq. (50)
Σ
(Dw.) 2
an
(49)
kp
JE(n]
Indeed, by the assumption weight is not degenerate, get that for all vectors in R.,
a DV, = Σ (Dw.),
jEIm]
min
(V, 2
Σ
(Dw))
sE(n]
jE(n]
pk (Dwi)s
jEIn
pu on
n kp
o,
where the second step follows from Fact A.6, the third step follows from the incoherence of V, the fourth step follows from our claim (see Eq. (49)), and the last step follows from simple algebra.
Then, by the incoherence of V, we have
uk
Σ JEIn]
(Dw.), < Σ (Dw.);
(50)
JE(n]
Hence,
Σ (Dw.) 2 an
kp
JC(m)
Combining everything together, we get
E
Σ (Y3) (Dw.),)2) kj
jE[m]
Define
B B: = :== (DW:))(DD) < k
((ww)s
n
46

---

Published as a conference paper at ICLR 2025
kullg*
where the first step follows from sampling procedure and the second step follows from the assumption that ||W] S koycoop
co
Since all the random variables
(Y,) (Dw.)s
are independent, applying Matrix Chernoff we get that
--
Pr[ jE(m]
S (Dw:)), 8)
a kips
S 11
FuB
(1 8)
e
S n
1op
(1-8)(1-8) ()(1-8)
Picking 8 3. and union bounding over all with probability at least
2/ for i
4.
muin (Y DTW,W
4ki
as needed.
H.2 MAIN RESULT
In this section, we summarize main result.
Theorem H.2 (Main result, random initialization). Let n be defined as Definition D.3. There is an algorithm starts from random initialization runs log(1 /e) iterations and generates M, which is matrix in Rax" and
|M M S O( ( " NII
l(kk n log(1/)))
is the total running time.
Proof. We use Algorithm 5 to initialize Y. Then, we can use the proof of Lemma D.9, and T is changed
T (i E [n]] Omin (YT (TT D,Y) 0.250/n).
where 77 = pukt
However, because of this change, T' = 0, with high probability. Then, the same calculation follows as in Lemma D.9. Note that in this case, Lemma F.1 is not needed because S = 0. Then, we can use Lemma J.1 and Lemma J.2, directly.
I BOUNDING THE FINAL ERROR
In Section I.1, we express M M simpler form that easier for further analysis. In Section I.2, we prove that |M - M*] is bounded. Both of these are used to support the proof of our main theorem.
47

---

Published as a conference paper at ICLR 2025
1.1 REWRITE M- M
In this section, we start to simplity/rewrite M M
Claim 1.1. Let M E Raxx (p, -iinchheeentt (see Definition 1). Let M be defined as Theorem 4.6. Then, we have
M'M +UEV +UEVTAGY A +RY
Proof. We start expanding the difference by definition:
M M
=XTT_ - MM
-UEV (VQ A,) R)(QQ A,) U
-UEV TVQ AvQ++++,)) A,) R(VQ+A))
USV
UEV VQ(VQ + A,) +UEV (VQ A, +
A,)" + +,, UEVT
-UEV VEETVQQ -UEV +UEV -UIVIVQAJ + A A,) + +RVQ + A
UEV -UEV VQ4, +UEV +UEVTA,(VQ + AA) +RVQ+ A,) T
= UEQAJ +UEV A, + RYT,
where the first step follows from M TY the second step follows from X = UVVV + R. Y VQ Ay. and M UEV, the third step follows from simple algebra, the fourth step follows from the simple algebra. the fifth step follows from that for all matrices A, B, ()+BBT A BT. the sixth step follows from V E = is an orthogonal matrix and Q E R kx is rotation matrix,
and the last step follows from = VQ + Ag.
1.2 BOUNDING ||M - M
In this section, we bound JM M*
Claim 1.2. Let M E Rexn is (p, --iincoherent (see Definition 1). Let M be defined as Theorem 4.6. Then, we have
| - (2Ar)
Proof.
M -- = [UXQA] + UUVT Y,Y + IRT
< | + |UEV |||A, + IIIIII
x 0(1) || + ||]]||
< O(dist(Y, V)) Ay. IW c NII
I 2
+ AF. W N + A
J[W o N|
1
= (AF + WW o
N
+
2
1
W NI|
2
< W N|| +e
where the second step is due to the inequalities ||A + B| < |||++BB and |AB|| |||||]] B|| the
third step is supported by ||E|| = || ||Q] 1, | |||
(See Definition 1), ||]| =
the second last step follows from A > Ar. and the last step follows from L O 1
48

---

Published as a conference paper at ICLR 2025
J PROOF OF MAIN RESULT
We dedicate this section to the proof of Theorem 4.6. There are two parts of Theorem 4.6: the correctness part and the running time part. In Section J.1, present the proof of the correctness part of Theorem 4.6. In Section J.2, we display the proof of the running time part of Theorem 4.6.
J.1 CORRECTNESS PART OF THEOREM 4.6
Now, we start proving the correctness part of Theorem 4.6.
Lemma J.1 (Correctness part of Theorem 4.6). Suppose M B=x" is p-incoherent (see As- sumption ).. Assume that W has "-spectral gap (see Assumption 2) and (or, (o.,3)-oundde (see Assumption 3). Let y satisfy condition in Definition D.4.
There an algorithm (Algorithm 1) takes M + N E IR"*" as input, uses either SVD initialization or random initialization and runs in ((log(1/1) iterations and generates M, which is a matrix in Raxx and
I - M s O(a |[ N
Proof. By Lemma 1.1, we can prove, for any > 1
1
(i, -
+ 100a -1 0 miin (M*) | N)),
2 dist (2, V) v
+ 100g o (M*)-"; | o Nj.
(51)
2
Note that minn (M*) ST (see Definition 1), then the above statement becomes
1
dist(X,,,)) S
AF
|[ o NII,
2
1
dY, V) - 2
+ AF
|W NN].
where Ap == 5.,. By Lemma D.9 and Claim D.7,
X - UEV Y|lFF
A dist((y)) V) + A, IIW NII.
(52)
where Aa and Ay are defied as Definition D.5.
To promise the first term in A is less than 0.1 and using Lemma A.3, we need to choose (note that cy is defined Definition D.4)
7 S 20 poly (H, k) n
To promise the second term in A is less than 0.1, have choose
a k))
? S 20
For dYst(Y, V), let
P = arg min QEORXE
IYQ - VIJ
We define
V
and Y VPT - APT
Let Q = PT E Okxk and A, = T then
Y VQ A
49

---

Published as a conference paper at ICLR 2025
with ||Ay|| dY, YY, V).
We define
R= X-X--VVVY
Then Eq. (52) implies that
|R]| & dist(Y, + A, N
Let X = Xx+1 and Y = Yr, then
=
Using Claim L.1, have
M M = + +UEV A, + +RY..
(53)
Using Claim 1.2, we have
IM M (2AF)I N +
(54)
J.2 RUNNING TIME PART OF THEOREM 4.6
Now, we start proving the running time part of Theorem 4.6.
Lemma J.2 (Running Time Part of Theorem 4.6). The running time of Algorithm is O([I nk²) (gg(()))) with random inilialization.
Proof. Now we analyze the running time. We first compute the initialization time. The entry Yi,; of the matrix Y equal / bi,, where ba.y are independent uniform from f-1, Hence, the time complexity of random initialization is O(nk). There are T iterations. For each iteration, there are three major steps, solving regression, Clip and QR. The dominating step is to solve regression. We choose Esk Claim D.7. Using Lemma C.2 and Lemma C.3, we know that we should choose 6o Eake/ (nly(n) and 6, 1/ poly(n, log(1/ (1/e), this step takes O(I k + nk*) time. The CLIP and QR algorithms take time ((nk) and Onk?²) respectively. Hence, the T iterations take time O([[WWloo + n²*²) log(1/e)).
K EXPERIMENTAL RESULTS
We conducted two experiments showing the performance of our main algorithm (Algorithm 1) and experiment, particularly for novel nigh-precision regression algorithm (Algorithm 2). We first present the first two experiments for our main algorithm. In both experiments, we set
M = M N T
where M is the rank-k ground truth and N is higher-rank noise matrix. We 100. We generate the noise matrix N as an n x n random matrix with i...d. Gaussian entries of zero mean and variance J. We apply the sketching matrix S E Rmxa with m 150 (we choose the CountSketch matrix (Charikar et al., 2002) when solving for the regression problems (see Lines 7 and 10 from our Algorithm 1). We iterate the alternating minimization steps (see Line 6) for T 20 times. To show the performance of our algorithm, we compare the running time and error between our Algorithm and the exact solver from Li al. (2016).
800 and
Experiment The first experiment is the matrix completion problem. For each row of the weight matrix Ruxn we randomly select 400 entries to be equal to and the remaining 400 entries be O. The second experiment the general weighted low rank approximation where the weight matrix W is constructed via 1, G for G being a random matrix with standard Gaussian entries.
Below, we present our experimental results of the matrix completion problem. We generate the ground truth M XYT by pair X, Y E B with random i.i.d. entries scaled by and the distributions Laplace, Gaussian, and uniform, respectively. Our time is measured in seconds.
50

---

Published as a conference paper at ICLR 2025
Table 4: Experimental results of Algorithm for the matrix completion problem.
Distribution
Time of the ex- act solver
Time of the % speedup approximate solver
Exact enro
solver
Approximate solver error
Laplace 234
205
12.39%
1.54
10 10-4
3.56
10
Gaussian
247
216
12.55%
3.85
10-4
1.33 x 10
Uniform 248
207
16.53%
1.34 x 10-5
4.99 x 10-4
Table 5: Experimental results of Algorithm for the weighted low-rank approximation.
Distribution
Time of the act solver
Time of the % speedup approximate solver
Exact solver error
Approximate solver error
Laplace 249
222
10.84%
2.11 x 10
4.71 x 4
Gaussian
223
204
8.52%
6.29
10
2.15 10-4
Uniform
221
215
2.71%
3.00
10-5
7.13 x 10-3
Experiment We present our second experimental results as follows (recall the weight matrix is an all-1's matrix plus a noise matrix with i.i.d. standard Gaussian entries):
We note that in these two settings, our algorithm achieves speedup compared to the algorithm of Li et al. (2016). For the matrix completion setting, our speedup is in the range 12%-16%, while for the dense weights regime, speedup is roughly 10%. In 20 iterations, approximate solver obtains errors similar to those of the exact solver.
Experiment 3
We expect the speedup will be more significant once the discrepancy between and l larger, as sketching is known to work well in the regime where n > k, as evidenced by the following experiment using sketching solve the regression:
min Az b12 for A E paxk and E IR".
We test the performance of regression solvers for 10 k 500 with sketch size 5500, and run our solver for 5 iterations. The results are as follows:
Table 6: Experimental results of Algorithm 2.
Distribution
Time for exact solve
Time for approx solve
Percentage speedup
Error
Gaussian
5.149
3.695
28.24%
8.24 x 10-3
Laplace
5.480
3.866
29.45%
8.41 10-3
Power Law
5.505
4.444
19.27%
8.77 10-3
Our data matrices A and response vectors b are generated according to standard Gaussian, Laplace, and power law distribution with p = 5, and we measure the Éoo error the solution, i.c., let 2 denote the vector outputted by the solver, measure
|
11,0
where
x= = arg min A Z
bll2-
The speedup obtained here ranges from 20% to 30%, so we have strong grounds to believe that the acceleration will be even evident when is large. However, performing weighted low-rank approximation on 10 x 10th size matrices is currently out the scope of our computational power, we leave this future direction.
51


Reviews:
[{"review_id": "K6XFOXQZg7", "reviewer": "Reviewer_dS4q", "summary": "This paper studies the problem of weighted low-rank matrix approximation. In this problem one is given a weight matrix W\n and observes a matrix M', and the goal is to recover a matrix M of rank at most k that minimizes  $\\Vert W \\circ (M'- M)\\Vert$..\n\nThis paper provides new theoretical guarantees within the alternating minimization framework for weighted low-rank matrix approximation. The authors propose a new alternating minimization algorithm that uses high-accuracy regression solver based on sketching instead of exact regression solvers in the optimization step of the algorithm compared to the work of (Li et al, 2016). This allows authors to improve runtime by a factor of $k$. \n.", "questions": "Are you aware of any practical application where the improvement in the running time by a factor of k obtained in this paper may be significant to applications, and where one might expect assumptions to hold?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The problem of weighted low-rank matrix approximation is an important practical problem and an important algorithmic primitive that arises in numerical linear algebra, and various machine learning applications. The problem naturally generalizes classical low-rank matrix approximation problems and a matrix completion problem. Alternating minimization is an extremely successful practical approach to matrix completion problem, so establishing theoretical guarantees for this approach to the weighted low-rank matrix approximation problem is of high interest.\n\nThis paper builds on top algorithm analyzed in (Li et al., 2016) and the key difference between the algorithm in (Li et al., 2016) and this paper is that exact solvers for linear regression is replaced with an approximate solver that can be traced back to (Rokhlin & Tygert, 2008). This allows authors to improve runtime by a factor of $k$ and to get a nearly optimal runtime.", "weaknesses": "The assumption 2.4.2 asks that $||W - J||<\\gamma n$, where J is an all-one matrix with $\\gamma = O(n^{-c})$ (compared to O(1) in  (Li et al., 2016)). This assumption sounds quite strong to me. This means that in the matrix completion scenario, this allows one to have only o(1)-fraction of entries not observed in every row. Considering that one of the key regimes for low-rank problems assumes that the rank satisfies k = O(1), I am not sure if improvement in runtime completely justifies slightly stronger assumption on W.\n\nAnother slight weakness of the paper is that most of the key techniques used in this paper were introduced in prior work, so the paper lacks a bit of novelty in this sense.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Rwk97Ig6Y5", "reviewer": "Reviewer_MHu4", "summary": "The authors propose an efficient framework for alternating minimization, which allows the subproblems to be solved approximately. They apply the framework to the weighted low rank approximation problem, theoretically improving the runtime.", "questions": "1. Could the authors explain the connection between the problem they considered (Problem 2.4) and the original problem? Are they equivalent under some conditions?\n2. Could the authors explain the meaning of $\\epsilon_0$ in Line 351?", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The authors propose a high-precision regression solver and extend it to the weighted case. They also demonstrate the robustness of the alternating minimization framework. The theoretical analysis is a significant contribution.\n2. The motivation behind their method is detailed and well-presented in the paper.", "weaknesses": "The authors did not conduct experiments to verify the performance of their algorithm. I am interested in whether the algorithm is feasible in practice.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "13rvBCPRer", "reviewer": "Reviewer_j2kz", "summary": "This paper considers the weighted low-rank approximation problem, which can be solved using alternating minimization. In previous methods, the alternating minimization algorithm requires solving n linear regression problems each iteration. In this work, the authors develop a framework based sketching, and they prove that alternating minimization is robust enough such that the approximate solution return each iteration also converges towards the optimal solution. As a result, they obtain an improved running time compared to Li et al. 2016.", "questions": "Can the authors comment on the performance of their algorithm on some real or even synthetic problems?", "limitations": "", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "I think the theoretical results in this paper is solid: the proof that the alternating minimization framework is robust enough to tolerate large errors induced in the approximate solution is interesting, and can possibly be generalized to other similar problems in matrix optimization.", "weaknesses": "I see two main weaknesses. First, I think the writing of this paper needs to be improved, especially in section 3, where the algorithm is first introduced. Algorithm 1 is outlined in a dense way that makes it hard for the reader to parse exactly what is going on. Also, lines 270-280 contain numerous sentences that are a bit awkward or hard to understand. I suggest that the authors use a separate section to slowly present algorithm 1 step by step, and describe which parts are new. \n\nAnother weakness i see is the lack of any numerical results. In the introduction (lines 80) the authors say that the results of Li et al. 2016 are hard to deploy in practice because of the high time complexity. While this paper reduces the time complexity in theory, there are no numerical section to demonstrate that the reduction in time complexity is useful in practice. I suggest that the authors at least run some synthetic experiments to demonstrate a reduction in real runtime.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "OnqY4kSXJi", "reviewer": "Reviewer_5iC9", "summary": "In this work, the authors proposed a new alternating minimization method for the weighted low-rank matrix approximation problem. They utilized an efficient multiple-response regression solver to reduce the computational time from ${O}((\\|W\\|_0k^2 + nk^3)\\log(1/\\epsilon))$ to $\\tilde{O}((\\|W\\|_0k + nk^3)\\log(1/\\epsilon))$.", "questions": "I have a few other minor comments for the authors to consider:\n\n- Theorem 4.6: I wonder if the time complexity also depends on the scale of $M^*$ and $N$. For example, let $c>0$ be a constant. If we multiply $M^*$ and $N$ by $c$, then the solution $\\tilde{M}$ will also be multiplied by $c$. However, the value $\\epsilon$ is not changed in the error bound. In addition, Assumption 2.3 does not bound the scale of $M^*$ and $N$ either. If I did not miss anything, I feel that the number of iterations or the error $\\epsilon$ should depend on the scale of $M^*$ and $N$.\n\n\n- Line 499: it would be helpful if the authors could provide supporting evidence for the claim that the spectral norm is more robust than the Frobenius norm.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper provides a novel algorithm for the low-rank matrix approximation problem. The results in this work may also be applied to other optimization and machine learning problems. After checking the paper, I think the theoretical results in the paper should be sound, although I did not check all proofs due to the time limit.", "weaknesses": "Although this paper is mostly a theoretical paper, I think it is still beneficial to include numerical illustrations to verify the theoretical findings and show the usefulness of the proposed method. For example, the authors may generate artificial datasets that satisfy the assumptions in this work and compare their method with the methods in (Li et al., 2016). It would be even better if the authors can verify that Assumption 2.3 holds in some practical examples and exhibit the empirical performance on these real-world examples.\n\nMoreover, in Line 223, the authors claimed that Assumption 3 is a generalization to the RIP condition. Although I agree with the authors that the assumption is more general in the sense that it contains the weight matrix $W$, I think it is more restrictive than the RIP condition. Suppose that $W$ is a boolean matrix and each column of $W$ contains at most $k$ 1's. In this case, the assumption considers at most $n$ combinations of rows of $U$ and $V$. In contrast, the RIP condition requires the error bound to hold for all $k$-dimensional subspace. In addition, I feel that Assumption 3 is not directly comparable with the classical RIP condition. The definition of the classic $(r, \\delta)$-RIP condition would imply that $(1-\\delta) \\lVert K \\rVert_F^{2} \\leq \\lVert W \\circ K \\rVert_F^{2} \\leq (1+\\delta) \\lVert K \\rVert_{F}^{2}$ holds for any matrix $K$ of rank at most $k$, which is independent of $M^*$. I would suggest the authors include more detailed explanations to this claim.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)