Near-Optimal Nonconvex Matrix Completion Jian-Feng Cai∗
Xiliang Lu†
Juntao You‡
arXiv:2609.17048v1 [math.NA] 15 Sep 2026
Abstract We study nonconvex methods for matrix completion, the problem of recovering a lowrank matrix from a subset of its entries. Convex methods achieve sample complexity linear in the matrix dimension and the rank, up to logarithmic factors, whereas global guarantees for commonly used nonconvex methods require a higher polynomial dependence on the rank. We close this gap by analyzing Riemannian gradient descent (RGD) and Riemannian Gauss–Newton (RGN) methods. For an n × n matrix of rank r with incoherence parameter µ and condition number κ, the two methods achieve exact recovery with high probability from O(µnr log n log(nκ)) and O(µnr log n log(2µrκ)) observations, respectively. The methods use a multiscale residual initialization, while the analysis simultaneously controls the spectral error and incoherence. The resulting RGD iterates converge linearly, whereas RGN eventually converges Q-quadratically. Keywords. nonconvex matrix completion, sample complexity, Riemannian gradient descent, Riemannian Gauss–Newton method, multiscale initialization, leave-one-out analysis
1
Introduction
Matrix completion seeks to recover a low-rank matrix from a subset of its entries and arises in a broad range of problems in machine learning and data analysis, including collaborative filtering [18], dimensionality reduction and clustering [9], model reduction and system identification [20], and sensor network localization [9]. Given an unknown matrix X⋆ ∈ Rn×n of rank r and an observed index set Ω, a natural formulation is min rank(X)
X∈Rn×n
subject to
PΩ (X) = PΩ (X⋆ ),
(1)
where PΩ retains the entries indexed by Ω. Since rank minimization is computationally intractable in general, the seminal work of Candès and Recht [6] studied the convex relaxation obtained by replacing the rank with the nuclear norm. Candès and Tao [7] proved exact recovery from O(µ2s nr log6 n) randomly observed entries under the strong incoherence condition. Chen [8] later showed that O(µnr log2 n) observations suffice under standard incoherence. This bound is optimal in its dependence on µ, n, and r, up to logarithmic factors. Nevertheless, nuclear norm minimization can be computationally expensive. To reduce this computational cost, a number of efficient nonconvex methods have been developed that exploit the rank constraint directly. Their global recovery guarantees, however, generally require a higher polynomial dependence on r; see Table 1. A natural question is whether efficient nonconvex methods for matrix completion can attain a sample complexity linear in n and r, up to logarithmic factors. ∗
Department of Mathematics, Hong Kong University of Science and Technology. E-mail: [email protected]. School of Mathematics and Statistics, Hubei Center for Applied Mathematics, and Hubei Key Laboratory of Computational Science, Wuhan University, Wuhan 430072, China. E-mail: [email protected]. ‡ (corresponding author.) School of Artificial Intelligence, Hubei Center for Applied Mathematics, and Hubei Key Laboratory of Computational Science, Wuhan University, Wuhan 430072, China. E-mail: [email protected]. †
1
We address this question by establishing near-optimal sample complexity bounds for RGD and RGN equipped with a multiscale residual initialization. Under standard incoherence condition, we establish global recovery guarantees for RGD and RGN, with RGD converging linearly and RGN eventually converging Q-quadratically. With high probability, the two methods recover X⋆ with sample complexities O(µnr log n log(nκ))
and
O(µnr log n log(2µrκ)),
respectively, where κ = σ1 (X⋆ )/σr (X⋆ ) is the condition number. Related work. A common nonconvex approach to matrix completion is to factorize X = LR⊤ and optimize over the factors, as in OptSpace [18], alternating minimization [13, 17], and gradient-based methods [3, 23, 29]. For incoherent positive semidefinite matrices, gradient descent [21] converges linearly without explicit regularization while maintaining incoherence along the iterates, and scaled projected gradient descent [4, 22, 25, 32] further removes the dependence of the convergence rate on κ. Projected-gradient and hard-thresholding methods [2, 15, 16, 24] instead work directly with the matrix variable, but require a rank-r approximation after each gradient step. Riemannian methods [28, 31] avoid this large-scale truncation by restricting the search direction to the tangent space of the fixed-rank manifold. Second-order methods have also been studied, including MatrixIRLS [19] and Gauss–Newton methods [27, 33, 34]. Despite these developments, the known global recovery guarantees generally retain a higher polynomial dependence on the rank; for example, Riemannian gradient descent [31] requires O(µκ6 nr 2 log2 n) observations under certain conditions. In contrast, for Gaussian matrix sensing, Riemannian gradient descent [5] attains the optimal O(nr) sample complexity. For matrix completion, however, the sampling operator does not satisfy a uniform restricted isometry over low-rank matrices, and the iterates must be shown to remain incoherent. This makes the analysis more challenging. Another issue is the dependence on the condition number in initialization. The ordinary spectral estimator, widely used in nonconvex methods, can incur an additional factor κ2 in the sampling requirement (see Theorem 3.4 in this work). Stagewise methods such as SoftDeflate [14] and the projected-gradient method [16] reduce or remove this dependence, but the recovery guarantees still have a higher polynomial dependence on r. Our multiscale initialization instead controls the spectral, row, column, and entrywise errors simultaneously, which avoids an additional polynomial dependence on r in the sampling requirement. Our contributions and proof strategy. Our main contribution is to establish near-optimal global sample complexity bounds for efficient nonconvex RGD and RGN methods in matrix completion. As discussed above, two sources of additional sample complexity arise in the usual analysis: the conversion from spectral to Frobenius error can introduce an additional factor r, while ordinary spectral initialization can incur a factor κ2 . To control the loss in rank r, we keep track of the spectral, row, column, and entrywise errors simultaneously. For Z ∈ Rn×n , define kZk2,∞ := max eTi Z 2 , i
ZT
2,∞
:= max kZej k2 , j
and introduce the sharp-norm as r r 1 n 1 n kZk2,∞ , ZT kZk♯ := max kZkop , 2 µr 2 µr
kZk∞ := max |(Z)ij |, i,j
n kZk∞ . , 2,∞ 4µr
(2)
The multiscale initialization maintains this sharp-norm control and reaches the required local regions without introducing an additional polynomial dependence on r in the sampling requirement. For RGD, the eventual conversion to the Frobenius norm affects only the number of 2
Table 1: Complexity comparison for matrix completion under the stated recovery guarantees. All results assume incoherence, while Factorized GD [21] also assumes PSD and κ = O(1); Riemannian GD [31] additionally assumes spikiness. Method
Iteration complexity
Nuclear norm Polynomial minimization time; solver [7, 8, 10] dependent
Cost per iteration
Sample complexity
Total computational complexity
Solver dependent
O(µnr log(2µr) log n)
Solver dependent
O(|Ω|r + nr)
O(µ3 nr 3 log3 n)
O µ3 κ2 nr 4 log3 n log 1ε O n3 log 1ε
Factorized GD[21]
O κ2 log 1ε
SVP/PGD [10, 15, 30]
O log 1ε
O(n3 )
O(µ2 κ4 nr 2 log n)
ScaledPGD [25]
O log 1ε
O(|Ω|r + nr 2 )
O µκ2 nr 2 (µκ2 ∨ log n)
Riemannian O log 1ε GD[31]
O(|Ω|r + nr 2 )
RGD O log 1ε (this paper)
O(|Ω|r + nr 2 )
RGN 1 b + nr 2 ) Ω|r (this paper) O log log ε O Jk (|
O µκ2 nr 3 (µκ2 ∨ log n) log 1ε
O max{µ0 , µ21 }κ6 nr 2 log2 n O max{µ0 , µ21 }κ6 nr 3 log 2 n log 1ε
O(µnr log n log(nκ))
O µnr 2 log2 (nκ) log n + log 1ε
O(µnr log n log(2µrκ))
O µJnr 2 log2 n log(2µrκ) log 1ε
initialization steps; for RGN, the sharp-norm control is continued through the initial iterations before the analysis enters the Frobenius regime. To control the loss in κ2 , we use a multiscale residual initialization in place of the ordinary spectral estimator. Successive residual reconstructions reduce the sharp-norm error by a fixed factor at the sampling level p & µr log n/n, and we establish both the sampling and computational complexities of this procedure. We further show that the µκ2 r/n sampling scale of the ordinary spectral estimator is necessary over an explicit family of incoherent matrices. The proof combines these initialization estimates with the local convergence analysis. For RGD, the initialization reaches the required Frobenius neighborhood, where a uniform tangentspace sampling estimate yields linear convergence. For RGN, a finite leave-one-out argument propagates the sharp-norm control through the initial iterations; once the iterates enter a sufficiently small Frobenius neighborhood, a local deterministic argument yields quadratic convergence. Organization and notation. Section 2 introduces the observation model, the Riemannian algorithms, and the multiscale initialization. Section 3 states the global recovery guarantees and the lower bound for ordinary spectral initialization. Section 4 presents the proof framework, including the local convergence and initialization results. Their proofs, together with the proofs of the global theorems, are given in Sections 5–7. Sections 8 and 9 present the numerical experiments and concluding remarks. The appendices collect the supporting probabilistic and geometric estimates, the proof of the spectral lower bound, the sharp-norm estimates for spectral reconstruction, the leave-one-out analysis, and the implementation and complexity analysis. Throughout the paper, bold lowercase letters denote vectors and bold uppercase letters denote matrices, while scalars are written in ordinary type. The vector ei denotes the ith standard basis vector. For a vector x, kxk2 denotes the Euclidean norm. For a matrix Z, kZkop and kZkF denote the operator norm and Frobenius norm respectively. We write σi (Z) for the i-th largest singular value of Z, and rank(Z) and range(Z) for its rank and column space. The symbols I and I denote the identity matrix and identity operator, respectively. We use 3
O(·) for bounds up to an absolute numerical constant independent of the problem parameters.
2
Problem Formulation and Riemannian Algorithms
We first formulate the matrix completion problem and then describe the Riemannian algorithms considered in this paper.
2.1
Problem setup
Suppose that X⋆ ∈ Rn×n , n ≥ 2 is an unknown matrix of rank r, where 1 ≤ r < n. Assume we observe its entries independently according to the Bernoulli sampling model: ( (X⋆ )ij , with probability p, 1 ≤ i, j ≤ n, Yij = ∗, with probability 1 − p, where 0 < p ≤ 1. Let Ω := {(i, j) : Yij 6= ∗} denote the set of observed indices, also denoted by Ω ∼ Bernoulli(p). The associated sampling operator PΩ is defined by ( Zij , (i, j) ∈ Ω, (PΩ (Z))ij = 0, (i, j) ∈ / Ω. The matrix completion problem is to recover X⋆ from the observed entries PΩ (X⋆ ). Let X⋆ = U⋆ Σ⋆ V⋆T ,
Σ⋆ = diag(σ1 , . . . , σr ),
σ1 ≥ · · · ≥ σr > 0,
be a compact singular value decomposition of X⋆ . We assume that X⋆ satisfies the following standard incoherence condition. Assumption 2.1 (Incoherence [6–8]). For some 1 ≤ µ ≤ n/r, it holds µr 2 2 T T max max U⋆ ei 2 , max V⋆ ej 2 ≤ . i j n The standard incoherence condition was introduced by Candès and Recht [6] for low-rank matrix completion. It requires the left and right singular spaces of X⋆ to be weakly correlated with the canonical basis, preventing the matrix from being concentrated on only a few entries. For the completion problem, we consider the following nonconvex formulation: min fΩ (X),
X∈Mr
fΩ (X) :=
1 kPΩ (X − X⋆ )k2F , 2p
(3)
where Mr := {X ∈ Rn×n : rank(X) = r} is the manifold of rank-r matrices.
2.2
Riemannian gradient descent
We first present the RGD method for (3). The algorithm updates the current iterate along the negative gradient direction in the tangent space of the fixed-rank manifold and then retracts the tangent update back onto the manifold. We first recall the geometry of the fixed low-rank manifold. Let X = U ΣV T ∈ Mr be a compact singular value decomposition. The tangent space of Mr at X is [28] TX Mr = U Z1T + Z2 V T : Z1 , Z2 ∈ Rn×r , 4
and the orthogonal projection onto TX Mr is PTX (Z) = U U T Z + ZV V T − U U T ZV V T .
(4)
Since ∇fΩ (X) = p−1 PΩ (X − X⋆ ), the Riemannian gradient of fΩ at X is grad fΩ (X) = PTX ∇fΩ (X) = p−1 PTX PΩ (X − X⋆ ). A tangent update does not in general belong to Mr . We use the orthographic retraction [1] to map it back onto the fixed-rank manifold. For ξ ∈ TX Mr , define −1 T RetrX (ξ) := (X + ξ)V U T (X + ξ)V U (X + ξ), (5)
whenever U T (X + ξ)V is nonsingular. Let ξk = − grad fΩ (Xk ). Since h∇fΩ (Xk ), ξk i = − kξk k2F , exact line search along this tangent direction [31] gives αk = arg min fΩ (Xk + αξk ) = α∈R
kξk k2F
p−1 kPΩ (ξk )k2F
,
(6)
whenever ξk 6= 0. The RGD update is then Xk+1 = RetrXk (αk ξk ), as summarized in Algorithm 1. If ξk = 0, the algorithm terminates. Algorithm 1 Riemannian gradient descent (RGD) Input: p, PΩ (X⋆ ), and X0 ∈ Mr . 1: for k = 0, 1, 2, . . . do 2: ξk ← −p−1 PTXk PΩ (Xk ) − PΩ (X⋆ ) . 3: αk ← kξk k2F / p−1 kPΩ (ξk )k2F if ξk 6= 0; otherwise, stop. 4: Xk+1 ← RetrXk (αk ξk ). 5: end for The tangent gradient can be evaluated using sparse matrix–factor products in O(|Ω|r + nr 2 ) operations. The exact line-search stepsize can be evaluated in the same order by computing kξk kF and the entries of ξk on Ω. The orthographic retraction can be computed from low-rank factors using two thin QR factorizations and an r×r singular value decomposition in O(nr 2 +r 3 ) operations [1]. Thus, since r < n, each RGD iteration costs O(|Ω|r + nr 2 ) operations. The implementation details are given in Section F.
2.3
Riemannian Gauss–Newton
Riemannian Gauss–Newton first computes a search direction in the tangent space and then b ⊂ Ω denote the subset of observations used retracts the tangent update back onto Mr . Let Ω for the RGN iterations. At X ∈ Mr , since D RetrX (0)[ξ] = ξ,
ξ ∈ TX Mr ,
linearizing the sampled residual along the retraction gives the Gauss–Newton subproblem 1 2 PΩb (X − X⋆ + ξ) F . ξ∈TX Mr 2 min
(7)
Its first-order optimality condition is the tangent normal equation PTX PΩb PTX ξ = PTX PΩb (X⋆ − X). 5
(8)
When the sampled tangent normal operator is positive definite on TX Mr , (7) has a unique minimizer. We compute this direction by applying the conjugate gradient method to (8) on TX Mr . At the k-th nonterminal RGN iteration, let Jk ≥ 1 denote the number of CG iterations used to solve the normal equation. CG is started from zero and run to the exact solution. If the right-hand side is zero, the algorithm terminates; the sequence is then continued by its final iterate for the convergence statements. In exact case, the finite-termination property of CG [11, Theorems 2.3.2 and 3.1.1] gives Jk ≤ dim(TXk Mr ) = r(2n − r). The resulting iteration is summarized in Algorithm 2. Algorithm 2 Riemannian Gauss–Newton (RGN) Input: PΩb (X⋆ ) and X0 ∈ Mr . 1: for k = 0, 1, 2, . . . do 2: Apply Jk CG iterations, starting from zero, to solve the following normal equation for ξk : PTXk PΩb PTXk ξ = PTXk PΩb (X⋆ − Xk ). 3:
Xk+1 ← RetrXk (ξk ).
4: end for
We also consider a regularized variant of RGN. At the k-th iteration, define λk =
1 PTXk PΩb (X⋆ − Xk ) σr (Xk ) F
(9)
and add λk ξ to the left-hand side of the normal equation (8). At every nonterminal iteration, λk > 0, so the regularized normal equation is positive definite on TXk Mr . Using the same CG and termination conventions as above gives Algorithm 3. Algorithm 3 Regularized Riemannian Gauss–Newton Input: PΩb (X⋆ ) and X0 ∈ Mr . 1: for k = 0, 1, 2, . . . do 1 2: λk ← PTXk PΩb (X⋆ − Xk ) . σr (Xk ) F 3: Apply Jk CG iterations, starting from zero, to solve the following normal equation for ξk : PTXk PΩb PTXk ξ + λk ξ = PTXk PΩb (X⋆ − Xk ). Xk+1 ← RetrXk (ξk ). 5: end for 4:
The normal equations in Algorithms 2 and 3 can be solved matrix-free without forming the b + nr 2 ) operasampled tangent normal matrix. Each normal-operator application costs O(|Ω|r tions, and the regularization term in Algorithm 3 does not change this order. The orthographic retraction costs O(nr 2 + r 3 ) operations. Hence, since r < n and Jk ≥ 1, the k-th iteration of either method costs b + nr 2 ) O Jk (|Ω|r
operations. The matrix-free implementation and detailed complexity analysis are given in Section F.
6
2.4
Multiscale initialization
The local convergence results for both RGD and RGN require an initial point sufficiently close to X⋆ . By Theorem 3.4, the standard spectral estimator can require a sampling probability of order µκ2 r/n to reach a constant sharp-norm neighborhood. We instead use a multiscale residual spectral initialization, related to residual spectral updates in singular value projection and iterative hard thresholding [2, 15, 24] and to stagewise constructions for matrix completion [14, 16]. At each scale, the current approximation is corrected by the observed residual and then truncated at a decreasing spectral level. Let Ω(ℓ) ⊂ Ω denote the observations used at scale ℓ, with sampling probability q. Starting from Z0 = 0, set τ0 := 2q −1/2 kPΩ(0) (X⋆ )kF and τℓ := 4−ℓ τ0 . At iteration ℓ, the residual correction satisfies Zℓ + q −1 PΩ(ℓ+1) (X⋆ − Zℓ ) = X⋆ + q −1 PΩ(ℓ+1) − I (X⋆ − Zℓ ). Thus, the sampling perturbation acts on the current error X⋆ − Zℓ . Accordingly, we use the following update: 0 ≤ ℓ < K, (10) Zℓ+1 = Tτℓ Zℓ + q −1 PΩ(ℓ+1) (X⋆ − Zℓ ) ,
where Tτℓ denotes P the spectral truncation operator defined below. For a singular value decomposition Y = j σj uj vjT , first define the hard spectral thresholding operator H≥λ (Y ) :=
X
σj uj vjT .
σj ≥λ
For τ > 0, Tτ (Y ) is computed as follows: Starting from Q0 = qf(G), where G ∈ Rn×r has independent standard Gaussian entries, compute τ2 0 ≤ t < ⌈12 log n⌉ , (11) Qt , Qt+1 = qf Y (Y T Qt ) + 4096 2
τ where qf denotes the orthogonal factor in a thin QR factorization. The term 4096 Qt preserves T the eigenspaces of Y Y and keeps the block iteration well defined. Writing Q for the final factor, define Tτ (Y ) := QH≥τ /8 (QT Y ). (12)
Thus each reconstruction uses matrix–factor products and a compressed singular value decomposition. Its sharp-norm approximation property is established in Theorem 4.3, and the computational cost of one reconstruction is as follows. Proposition 2.2. Suppose that rank(Z) ≤ r and |Λ| = m. Then Tτ Z + q −1 PΛ (X⋆ − Z) can be computed in matrix-free form using O (mr + nr 2 ) log n operations. Evaluating the observed residual costs O(mr) additional operations. Proof. The proof is given in Subsection F.1. As for the topping rule in the above initialization stage, we adopt the following residual test. For the rank-r iterates, we use the observed residual Rℓ := q −1/2 kPΩ(ℓ+1) (X⋆ − Zℓ )kF , and stop the initialization at Rℓ ’s first increase or after K reconstructions. The output is denoted by ZKb . The procedure is summarized in Algorithm 4. 7
Algorithm 4 Multiscale initialization Input: rank r, sampling probability q, maximum number of reconstructions K, and {PΩ(ℓ) (X⋆ )}K ℓ=0 . 1: Z0 ← 0, τ0 ← 2q −1/2 kPΩ(0) (X⋆ )kF , j ← ∅. 2: for ℓ = 0, . . . , K − 1 do 3: Zℓ+1 ← Tτℓ Zℓ + q −1 PΩ(ℓ+1) (X⋆ − Zℓ ) . 4: τℓ+1 ← τℓ /4. 5: If ℓ + 1 < K, set Rℓ+1 ← q −1/2 kPΩ(ℓ+2) (X⋆ − Zℓ+1 )kF . 6: Stop and return Zj if ℓ + 1 < K, j 6= ∅, rank(Zℓ+1 ) = r, and Rℓ+1 > Rj . 7: Set j ← ℓ + 1 if ℓ + 1 < K and rank(Zℓ+1 ) = r. 8: end for 9: return ZK . We now specify the observation sets used in the analysis. Set ( K + 1, for RGD, B := q := 1 − (1 − p)1/B . K + 2, for RGN, The two choices of the upper bound K are specified in Section 3. For each (i, j) ∈ Ω, independently draw bij ∈ {0, 1}B \ {0} from the product Bernoulli(q) distribution conditioned on being nonzero, and set Ω(ℓ) := {(i, j) ∈ Ω : (bij )ℓ = 1} for 0 ≤ ℓ ≤ B − 1. Under the unconditional S (ℓ) law, Ω(0) , . . . , Ω(B−1) are mutually independent Bernoulli(q) subsets with Ω = B−1 ℓ=0 Ω , where the subsets may not be disjoint. The first K + 1 components are used for initialization. For b := Ω(K+1) is reserved for the subsequent iterations. RGN, the last component Ω
Remark 2.3. The above decomposition is used only in the analysis. In numerical implementation, the same observation set Ω is used throughout the algorithms, with q replaced by p.
3
Main Results
In this section, we state the global recovery guarantees for RGD and RGN, and then give a lower bound for ordinary spectral initialization.
3.1
Global convergence of RGD
For RGD, let X0 be the output of Algorithm 4 with the upper bound √ K = 5 + log4 (κ nr) .
Starting from X0 , let {Xk }k≥0 be the iterates of Algorithm 1 on Ω. The following theorem establishes linear convergence from this initialization.
Theorem 3.1 (Global convergence of RGD). Suppose that Assumption 2.1 holds and Ω ∼ Bernoulli(p). Let X0 be the output of Algorithm 4 with K specified above, and let {Xk }k≥0 be generated by Algorithm 1. If µr log n log(nκ) p ≥ C1 , n where C1 > 0 is a sufficiently large absolute constant, then, with probability at least 1 − n−10 , the initialization and all subsequent iterates are well defined, and k 3 kX0 − X⋆ kF , kXk − X⋆ kF ≤ 8 8
k ≥ 0.
Proof. The proof is deferred to Subsection 6.5. Thus O(µnr log n log(nκ)) observations suffice for exact recovery with high probability. After Algorithm 4, relative Frobenius accuracy ε is attained within O(log(1/ε)) RGD iterations. The complete initialization and computational costs are given in Theorem F.1.
3.2
Global convergence of RGN
For RGN, set l m K = 6 + log4 (µr 3/2 κ) ,
q = 1 − (1 − p)1/(K+2) ,
K̄ := ⌈log2 log2 (4n)⌉ .
We first run Algorithm 4 with this upper bound K using the RGN subsets specified in b The following theorem establishes global Subsection 2.4, and then run Algorithm 2 on Ω. convergence and eventual quadratic convergence of the resulting RGN iterates.
Theorem 3.2 (Global convergence of RGN). Suppose that Assumption 2.1 holds and Ω ∼ Bernoulli(p). Let X0 denote the output of Algorithm 4 with K specified above, and let {Xk }k≥0 b If be generated by Algorithm 2 on Ω. p ≥ C2
µr log n log(2µrκ) , n
where C2 > 0 is a sufficiently large absolute constant, then, with probability at least 1 − n−10 , the initialization and all subsequent iterates are well defined. Moreover, σr (X⋆ ) kXk − X⋆ k♯ ≤ 280µr and
7 25
2k
,
kXk − X⋆ k2F kXk+1 − X⋆ kF ≤ 4 √ , q σr (X⋆ )
0 ≤ k ≤ K̄,
k ≥ K̄.
Proof. The proof is deferred to Subsection 7.5. The RGN can also be regularized and admits a similar global recovery guarantee with the eventual quadratic rate. Corollary 3.3. Under the conditions of Theorem 3.2, let {Xk }0≤k≤K̄ be the RGN iterates b Then, with in Theorem 3.2. Starting from XK̄ , replace Algorithm 2 by Algorithm 3 on Ω. probability at least 1 − n−10 , all subsequent iterates are well defined and satisfy kXk − X⋆ k2F , kXk+1 − X⋆ kF ≤ C̄2 √ q σr (X⋆ )
k ≥ K̄,
where C̄2 > 0 is an absolute constant. Proof. The proof is deferred to Subsection 7.6. Thus O(µnr log n log(2µrκ)) observations suffice for both RGN and its regularized variant to converge to X⋆ with high probability, with quadratic convergence after finitely many iterations. Relative Frobenius accuracy ε is attained within O(log log(1/ε)) subsequent RGN iterations. The computational cost additionally depends on the CG iteration counts Jk ; see Theorem F.1.
9
3.3
Lower bound for ordinary spectral initialization
We next explain why the multiscale initialization is needed in place of the ordinary spectral estimator. The latter may require an additional factor κ2 in the sampling probability to reach the sharp-norm neighborhood required by the local convergence results. For Y ∈ Rn×n , let Hr (Y ) ∈ arg
min
rank(Z)≤r
kY − ZkF
denote a fixed best rank-r approximation. The following theorem shows that a sampling probability of order µκ2 r/n is necessary over an explicit family of incoherent matrices, with the initialization error measured in the sharp norm (2). √ Theorem 3.4. Let r, s ≥ 2 be integers with rs ≤ n, set µ = n/(rs), and let 1 ≤ κ ≤ s. There exist an absolute constant C3 > 0 and a rank-r matrix X⋆ satisfying Assumption 2.1 with coherence parameter µ and condition number κ such that, for every 0 < q < C3
µκ2 r , n
the spectral estimate
satisfies
X0 = Hr q −1 PΛ (X⋆ ) ,
Λ ∼ Bernoulli(q),
1 1 Pr kX0 − X⋆ k♯ ≥ σr (X⋆ ) ≥ . 2 2
Proof. The proof is deferred to Section A. The theorem identifies the µκ2 r/n sampling scale for the ordinary spectral estimator over this family. In Algorithm 4, reconstruction is applied to successive residuals whose sharp error decreases geometrically.
4
Proof Framework: Local Convergence and Initialization
The global results are obtained by combining local convergence of the fixed-sample iterations with the multiscale initialization. We first state the local results for RGD and RGN and then show that Algorithm 4 reaches their respective hypotheses.
4.1
Local convergence
The two local results use different error controls: 1 √ p σr (X⋆ ), 128 σr (X⋆ ) , kX − X⋆ k♯ ≤ 1000µr
kX − X⋆ kF ≤
for RGD,
(13a)
for RGN.
(13b)
Thus RGD is controlled in a Frobenius neighborhood whose radius depends on the sampling probability, whereas RGN requires the stronger sharp-norm control. A common geometric identity underlying both analyses is the exactness of the graph retraction for the population tangent correction. If kX − X⋆ k♯ ≤ 41 σr (X⋆ ), then a later established Lemma 5.3 shows that the graph core is invertible and RetrX (PTX (X⋆ − X)) = X⋆ . 10
(14)
Hence the local analysis reduces to controlling the error introduced by sampling. For a nonterminal RGD step, let ξ = PTX p−1 PΩ (X⋆ − X) and let α be the stepsize in (6). Then αξ − PTX (X⋆ − X) = (α − 1)PTX (X⋆ − X) + αPTX (p−1 PΩ − I)(X⋆ − X). The tangent sampling isometry controls both terms, which are first order in the error. For RGN, if ξ denotes the tangent correction, then PTX q −1 PΩb PTX (ξ − PTX (X⋆ − X)) = PTX q −1 PΩb PT ⊥ (X⋆ − X). X
(15)
The right-hand side is driven by the normal component, which is quadratic in the error. Indeed, under the local sharp-norm condition, a later established Lemma 5.2 gives 2
16 kX − X⋆ k♯ . PT ⊥ (X⋆ − X) ≤ X 3 σr (X⋆ ) ♯
(16)
The two mechanisms are summarized in Figure 1. True-tangent sampling isometry Lemma 5.6
Quadratic normal component PT ⊥ (X⋆ − X) X
♯
≤
(16/3) kX − X⋆ k2♯ /σr (X⋆ ) Lemma 5.2
Uniform tangent transfer and exact line search αξ − PTX (X⋆ − X) Lemma 5.7
Sampled tangent correction and local conditioning PTX q −1 PΩb PT ⊥ (X⋆ − X) X Lemma 7.4
Linear RGD contraction
Quadratic RGN convergence
Figure 1: Local mechanisms for RGD and RGN. For RGD, the tangent sampling isometry controls the stepsize and the first-order tangent perturbation. For RGN, the sampled correction is driven by the quadratic normal component. Finite leave-one-out control brings the iterates into a Frobenius neighborhood where quadratic convergence follows deterministically. For RGD, the high-probability event is uniform over the entire Frobenius neighborhood in (13a). Consequently, the initial point need not be independent of the observation set in the analysis. Theorem 4.1 (Local convergence of RGD). Suppose that Assumption 2.1 holds and Ω ∼ n , where C4 > 0 is a sufficiently large absolute constant, then, Bernoulli(p). If p ≥ C4 µr log n with probability at least 1 − n−10 /48, the following holds simultaneously for every X0 ∈ Mr satisfying 1 √ kX0 − X⋆ kF ≤ p σr (X⋆ ) : 128 Algorithm 1 on Ω is well defined, and its iterates satisfy k 3 kX0 − X⋆ kF , k ≥ 0. kXk − X⋆ kF ≤ 8 Proof. The proof is deferred to Subsection 5.3. b in the analysis. This independence For RGN, the initial point need be independent of Ω permits the leave-one-out argument used to control the initial RGN iterates, after which the quadratic Frobenius recursion applies. 11
b ∼ Theorem 4.2 (Local convergence of RGN). Suppose that Assumption 2.1 holds and Ω σr (X⋆ ) b Bernoulli(q). Let X0 ∈ Mr be independent of Ω and satisfy kX0 − X⋆ k♯ ≤ 1000µr . Set
n K̄ := ⌈log2 log2 (4n)⌉. If q ≥ C4 µr log , then, conditional on X0 , with probability at least n −10 b 1 − n /12 over Ω, the iterates of Algorithm 2 are well defined and satisfy, simultaneously for all 0 ≤ k ≤ K̄, k σr (X⋆ ) 7 2 kXk − X⋆ k♯ ≤ . 280µr 25
Moreover, for every k ≥ K̄, kXk − X⋆ k2F kXk+1 − X⋆ kF ≤ 4 √ . q σr (X⋆ ) Proof. The proof is deferred to Subsection 7.4.
4.2
Multiscale initialization
We now show that Algorithm 4 reaches the two local conditions above. For a fixed approximation Z and an independent observation subset Λ, the residual correction satisfies Z + q −1 PΛ (X⋆ − Z) = X⋆ + (q −1 PΛ − I)(X⋆ − Z).
(17)
Thus its sampling perturbation is determined by the current error, rather than by the largest singular value of X⋆ . The first result bounds one spectral reconstruction. Theorem 4.3. Suppose that Assumption 2.1 holds. Let Z be fixed with rank(Z) ≤ r and kZ − X⋆ k♯ ≤ τ , where τ > 0. Let Λ ∼ Bernoulli(q) be independent of the Gaussian matrix used in Tτ . If µr log n q ≥ C5 , n where C5 > 0 is a sufficiently large absolute constant, then Z + = Tτ Z + q −1 PΛ (X⋆ − Z) satisfies, with probability at least 1 − n−12 /24,
1 Z + − X⋆ ♯ ≤ τ. 4
rank(Z + ) ≤ r,
(18)
Proof. The proof is deferred to Subsection 6.2. The proof controls the sampling perturbation and its products with the true singular spaces before estimating the reconstructed matrix. The finite computation in (11) attains the required accuracy without a gap between adjacent singular values. The supporting estimates are proved in Subsection B.2 and section D. Theorem 4.4 (Multiscale initialization). Suppose that Assumption 2.1 holds and 1 ≤ K ≤ n. Let {Zℓ }K b be the output of Algorithm 4. If ℓ=0 denote the complete sequence in (10), and let ZK b ≤ K, q ≥ C5 µr log n/n, then, with probability at least 1 − n−10 /3, Algorithm 4 is well defined, K and kX⋆ kF ≤ τ0 ≤ 3 kX⋆ kF ,
kZℓ − X⋆ k♯ ≤ 4−ℓ τ0 ,
rank(Zℓ ) ≤ r,
√ In particular, K = ⌈5 + log4 (κ nr)⌉ gives a rank-r output satisfying ZKb − X⋆ F ≤
1 √ q σr (X⋆ ), 128
12
0 ≤ ℓ ≤ K.
(19)
(20)
whereas K = 6 + log4 (µr 3/2 κ) gives a rank-r output satisfying ZKb − X⋆ ♯ ≤
σr (X⋆ ) . 1000µr
(21)
Proof. The proof is deferred to Subsection 6.4. The two choices of K allow at most O(log(nκ)) and O(log(2µrκ)) reconstructions, respectively. For RGD, q ≤ p and (20) imply the hypothesis of Theorem 4.1. For RGN, (21) gives the sharp-norm hypothesis of Theorem 4.2, while the initialization, including the stopping decision, b The following proposition justifies the residual test, which is independent of the reserved set Ω. shows a residual increase cannot cause a return before the corresponding local convergence condition is satisfied. Proposition 4.5. Under the assumptions of Theorem 4.4, with probability at least 1 − n−10 /6, if Algorithm 4 returns at the residual test, then its output satisfies ZKb − X⋆ F ≤
σr (X⋆ ) . 1024n
(22)
Proof. The proof is given in Subsection 6.3.
5
Local Convergence of RGD
We first establish the geometric estimates used in the convergence and initialization arguments. We then prove a sampling estimate that holds uniformly over a Frobenius neighborhood of X⋆ and use it to prove Theorem 4.1. The initialization results are proved in Section 6.
5.1
Deterministic geometry
We begin with bounds for the singular subspaces and the normal component of the error, followed by perturbation bounds for the graph retraction. The proofs are given in Section C. For X = U ΣV T ∈ Mr , define GX := U T X⋆ V ,
LX := (I − U U T )X⋆ V ,
RX := (I − V V T )X⋆T U .
(23)
We also write NX := PT ⊥ (X⋆ − X)
tX := PTX (X⋆ − X),
X
(24)
for the tangent and normal components of the error. Here GX is the representation of X⋆ in the current singular bases, while LX and RX describe its components outside the current singular subspaces. The following lemma bounds the incoherence of these bases and the variation of the associated projectors. Its operator-norm estimates follow the argument in [31, Lemma 4.1], while the rowwise estimates use the sharp norm. Lemma 5.1. Suppose that Assumption 2.1 holds. Let X = U SV T ∈ Mr be a compact orthonormal factorization, and set e♯ := kX − X⋆ k♯ . If e♯ ≤ σr (X⋆ )/8, then r o n µr . max kU k2,∞ , kV k2,∞ ≤ 2 n
The singular-space projectors satisfy the rowwise bounds n max U U T − U⋆ U⋆T 2,∞ , V V T − V⋆ V⋆T Moreover,
PTX − PTX⋆ F→F ≤ 13
2,∞
o
32 ≤ 7
16 e♯ . 7 σr (X⋆ )
r
µr e♯ . n σr (X⋆ )
Proof. The proof is deferred to Subsection C.2.1. The next lemma gives a sharp-norm counterpart of the quadratic normal-component estimate in [31, Lemma 4.1], together with the graph-factor representation used below. Lemma 5.2. Suppose that Assumption 2.1 holds. Let X = U SV T ∈ Mr be a compact orthonormal factorization, and set e♯ := kX − X⋆ k♯ . If e♯ ≤ σr (X⋆ )/8, then GX is invertible, and the graph factors in (23) satisfy T NX = LX G−1 X RX ,
and
(25)
o n max kLX kop , kRX kop ≤ e♯ , r n o µr e♯ . max kLX k2,∞ , kRX k2,∞ ≤ 4 n Moreover, the normal component satisfies G−1 X op ≤
4 , 3σr (X⋆ )
(26)
2 16 e♯ kNX k♯ ≤ . 3 σr (X⋆ )
Its row, column, and entrywise norms satisfy n o 16 r µr e2 ♯ T , max kNX k2,∞ , NX 2,∞ ≤ 3 n σr (X⋆ )
2 64µr e♯ kNX k∞ ≤ . 3n σr (X⋆ )
Proof. The proof is deferred to Subsection C.2.2.
The following identity shows that the population tangent correction recovers X⋆ exactly. It is the local inverse formula for the orthographic retraction [1, Sec. 3.2.3]; its proof verifies the required core invertibility under the stated error bound. Lemma 5.3. Let X ∈ Mr satisfy kX − X⋆ k♯ ≤ σr (X⋆ )/4. Then RetrX (tX ) = X⋆ . Proof. The proof is deferred to Subsection C.2.3. We next bound the effect of perturbing the population tangent correction. Lemma 5.4. Suppose that Assumption 2.1 holds. Let X ∈ Mr , set e♯ := kX − X⋆ k♯ , and let η ∈ TX Mr . If e♯ ≤ σr (X⋆ )/8 and kηk♯ ≤ σr (X⋆ )/16, then RetrX (tX + η) is well defined and kRetrX (tX + η) − X⋆ − ηk♯ ≤ 13
e♯ kηk♯
σr (X⋆ )
+6
kηk2♯
σr (X⋆ )
.
Proof. The proof is deferred to Subsection C.2.4. The following Frobenius estimates will be used for local RGD convergence and for the eventual quadratic convergence of RGN. Part (i) is a variant of [31, Lemma 4.1]. Lemma 5.5. Let X ∈ Mr . Then the following statements hold. (i) Suppose that kX − X⋆ kF ≤ σr (X⋆ )/4. Then PT ⊥ (X⋆ − X) X
F
≤2
kX − X⋆ k2F , σr (X⋆ )
PTX − PTX⋆ F→F ≤ 3
kX − X⋆ kF . σr (X⋆ )
(ii) Let η ∈ TX Mr . Suppose that kX − X⋆ kF ≤ σr (X⋆ )/40 and kηkF ≤ σr (X⋆ )/320. Then the graph retraction is well defined and its nonlinear remainder satisfies kRetrX (tX + η) − X⋆ − ηkF ≤
11 kX − X⋆ kF kηkF 11 kηk2F + . 5 σr (X⋆ ) 10 σr (X⋆ ) 14
Proof. The proof is deferred to Subsection C.3. We next establish the standard sampling isometry at the tangent space of X⋆ [8, 31] and transfer it uniformly to nearby tangent spaces.
5.2
Sampling isometry and uniform transfer
Lemma 5.6. Suppose that Assumption 2.1 holds, and let Λ ∼ Bernoulli(pΛ ). If pΛ ≥ c1 µr log n/n, where c1 > 0 is a sufficiently large absolute constant, then, with probability at least 1 − n−10 /48, 1 PTX⋆ p−1 . Λ PΛ PTX⋆ − PTX⋆ F→F ≤ 16 Proof. The proof is deferred to Subsection B.1. We next transfer the preceding isometry to nearby tangent spaces. A related local tangentspace estimate appears in [31, Lemma 4.2]. The following deterministic lemma holds simultaneously throughout the stated Frobenius neighborhood, so the matrix X may depend on Λ. Lemma 5.7. Let Λ be an index set and let 0 < pΛ ≤ 1. Suppose that PTX⋆ p−1 Λ PΛ PTX⋆ − PTX⋆ F→F ≤
1 . 16
(27)
Then, simultaneously for all X ∈ Mr satisfying kX − X⋆ kF ≤
1√ pΛ σr (X⋆ ), 40
and all ζ ∈ TX Mr , we have 2 5 kζk2F ≤ ζ, p−1 kζk2F . Λ PΛ ζ ≤ 3 4 √ If, in addition, kX − X⋆ kF ≤ pΛ σr (X⋆ )/128, then
(28)
9 7 kζk2F ≤ ζ, p−1 kζk2F . (29) Λ PΛ ζ ≤ 8 8 √ Proof. Fix any X ∈ Mr satisfying kX − X⋆ kF ≤ pΛ σr (X⋆ )/40 and any ζ ∈ TX Mr . By homogeneity, it suffices to consider kζkF = 1. The assumed neighborhood is contained in kX − X⋆ kF ≤ σr (X⋆ )/40, and hence Lemma 5.5 applies. Since PTX ζ = ζ, we obtain ζ − PTX⋆ ζ F = PTX − PTX⋆ ζ F ≤ PTX − PTX⋆ F→F ≤3
kX − X⋆ kF 3√ pΛ , ≤ σr (X⋆ ) 40
PTX⋆ ζ F ≥ 1 − ζ − PTX⋆ ζ F ≥ 1 − Applying (27) to PTX⋆ ζ ∈ TX⋆ Mr gives r
15 −1/2 PTX⋆ ζ F ≤ pΛ PΛ PTX⋆ ζ F ≤ 16
15
r
3√ pΛ . 40
17 PTX⋆ ζ F . 16
(30)
By (27), (30), the triangle inequality, and kPΛ kF→F ≤ 1, we have −1/2
pΛ
−1/2
kPΛ ζkF ≥ pΛ r
−1/2
PΛ PTX⋆ ζ F − pΛ
PΛ ζ − PTX⋆ ζ
15 −1/2 PTX⋆ ζ F − pΛ ζ − PTX⋆ ζ F 16 r 2 ≥ . 3 ≥
F
F
By (27), (30), PTX⋆ ζ F ≤ 1, and the triangle inequality, we also have −1/2
pΛ
−1/2
kPΛ ζkF ≤ pΛ r
−1/2
PΛ PTX⋆ ζ F + pΛ
PΛ ζ − PTX⋆ ζ
17 −1/2 PTX⋆ ζ F + pΛ ζ − PTX⋆ ζ F 16 r 5 ≤ . 4 ≤
Finally, since PΛ is an orthogonal projector, the above inequalities give 2 5 2 −1 ≤ ζ, p−1 Λ PΛ ζ = pΛ kPΛ ζkF ≤ . 3 4 Rescaling the ζ proves (28). Since X was arbitrary in the Frobenius neighborhood specified in Lemma 5.7, (28) holds simultaneously throughout that neighborhood. √ To prove (29), suppose that kX − X⋆ kF ≤ pΛ σr (X⋆ )/128 and kζkF = 1. By Lemma 5.5 and orthogonality, we have ζ − PTX⋆ ζ F ≤
3 √ pΛ , 128
2
2
PTX⋆ ζ F = 1 − ζ − PTX⋆ ζ F ≥ 1 −
9 . 1282
The preceding comparison with the true tangent space now gives r r 7 9 −1/2 ≤ pΛ kPΛ ζkF ≤ . 8 8 Squaring and rescaling proves (29).
5.3
Proof of Theorem 4.1
Proof. For C4 ≥ c1 , Lemma 5.6 implies that, with probability at least 1 − n−10 /48, PTX⋆ p−1 PΩ PTX⋆ − PTX⋆ F→F ≤
1 . 16
On the event supplied by Lemma 5.6, fix any X ∈ Mr satisfying kX − X⋆ kF ≤
1 √ p σr (X⋆ ), 128
and set E := X⋆ − X and ξ := PTX p−1 PΩ E. Applying (29) and using the self-adjointness of PTX p−1 PΩ PTX , we obtain 1 PTX p−1 PΩ PTX − PTX F→F ≤ , 8 16
3 p−1/2 kPΩ PTX kF→F ≤ √ . 2 2
(31)
If ξ 6= 0, then (29) gives p−1 kPΩ ξk2F = ξ, p−1 PΩ ξ ≥
7 kξk2F > 0. 8
Thus the stepsize in (6) is well defined and satisfies 8 8 ≤α≤ , 9 7
|α − 1| ≤
1 . 7
(32)
For ξ = 0, set α = 1 in the following estimates, so that (32) still holds. Set η := αξ − tX . Decomposing E into its tangent and normal components, we have η = (α − 1)tX + α PTX p−1 PΩ PTX − PTX tX + αPTX p−1 PΩ PT ⊥ E. (33) X
The assumed neighborhood is contained in kEkF < σr (X⋆ )/4. Hence part (i) of Lemma 5.5 applies, and (31)–(33) give 3α α kEkF + √ kηkF ≤ |α − 1| + PT ⊥ E X 8 F 2 2p ≤
kEk2F 24 1 2 kEkF + √ √ ≤ kEkF . 7 3 7 2 p σr (X⋆ )
(34)
In particular, kEkF ≤ σr (X⋆ )/128 < σr (X⋆ )/40 and kηkF ≤ σr (X⋆ )/384 < σr (X⋆ )/320. Consequently, part (ii) of Lemma 5.5 applies. Since αξ = tX + η, we obtain 11 kηkF 11 kEkF kRetrX (αξ) − X⋆ kF ≤ kηkF 1 + + 5 σr (X⋆ ) 10 σr (X⋆ ) (35) 1 11 11 3 ≤ + 1+ kEkF ≤ kEkF . 3 640 3840 8 If ξ = 0, then RetrX (αξ) = X, so (35) implies X = X⋆ . Thus every terminal iterate in the stated neighborhood equals X⋆ . Since X was arbitrary in the Frobenius ball of Theorem 4.1, (35) holds simultaneously throughout that ball. In particular, the neighborhood is invariant under the RGD update. Therefore, if X0 satisfies the hypothesis of the theorem, induction gives k 3 kX0 − X⋆ kF , k ≥ 0. kXk − X⋆ kF ≤ 8 The same induction verifies the hypotheses of Lemma 5.5 and the positivity of every nonterminal line-search denominator, and therefore guarantees that all subsequent iterates are well defined. This proves the theorem.
6
Multiscale Initialization and Global Convergence of RGD
We now prove Theorems 4.3 and 4.4. We first bound the approximate reconstruction error in the sharp norm, then justify the residual stopping rule and verify the two local convergence conditions. Combining these estimates with Theorem 4.1 also proves the global RGD result.
6.1
Sharp spectral reconstruction
For a fixed error matrix E and an independent subset Λ ∼ Bernoulli(q), the sampling perturbation is W = (q −1 PΛ − I)E. To control both row and column errors, we use the symmetric dilation 0 W U⋆ 0 W := , Q⋆ := . (36) WT 0 0 V⋆ 17
The columns of Q⋆ are orthonormal, and Assumption 2.1 gives kQ⋆ k2,∞ ≤ ing lemma controls powers of the sampling perturbation on these columns.
p
µr/n. The follow-
Lemma 6.1. Suppose that Assumption 2.1 holds. Let E be fixed with kEk♯ ≤ τ , where τ > 0. If q ≥ c2 µr log n/n, where c2 > 0 is a sufficiently large absolute constant, then, with probability at least 1 − n−12 /48 over Λ, r µr τ j τ j , j ≥ 0. (37) , W Q⋆ 2,∞ ≤ 2 kW kop ≤ 512 n 64 Proof. The proof is deferred to Subsection B.2. We next give a deterministic reconstruction estimate. Its hypotheses concern approximate singular factors, rather than an exact truncated SVD. This distinction allows the computation in (11) to stop after a prescribed number of steps. Lemma 6.2. Suppose that Assumption 2.1 holds. Let Y = X⋆ + W , and suppose that (37) bΣ b Vb T be a compact holds for some τ > 0. Let Z + have rank at most r. If Z + 6= 0, let Z + = U singular value decomposition and suppose that b ≥ τ, σmin (Σ) 8 r µr τ bΣ b ≤ Y Vb − U , 64 n op
b = Vb Σ, b Y TU
Y − Z + op ≤
7τ . 32
(38a) (38b)
If Z + = 0, suppose instead that kY kop ≤ 7τ /32. Then kZ + − X⋆ k♯ ≤ τ /4. Proof. The proof is deferred to Subsection D.1. The iteration in (11) is a randomized subspace iteration; see, e.g., [12]. The next lemma verifies the approximation conditions for Tτ . Only the singular components above the current error scale must be accurately represented. No gap between adjacent singular values is assumed. Lemma 6.3. Let Y ∈ Rn×n satisfy σr+1 (Y ) ≤ τ /512, where τ > 0. Then Z + = Tτ (Y ) has rank at most r and, with probability at least 1−n−12 /48 over its Gaussian initial matrix, satisfies (38). If Z + = 0, the conclusion is kY kop ≤ 7τ /32. Proof. The proof is deferred to Subsection D.2.
6.2
Proof of Theorem 4.3
Proof. Set W = (q −1 PΛ − I)(X⋆ − Z),
Y = X⋆ + W .
Choose C5 ≥ c2 . By Lemma 6.1, except with probability n−12 /48, (37) holds,
σr+1 (Y ) ≤ kW kop ≤
τ . 512
Conditional on such a realization of Y , Lemma 6.3 gives, except with probability n−12 /48, rank(Tτ (Y )) ≤ r,
(38) holds.
Hence Lemma 6.2 yields kTτ (Y ) − X⋆ k♯ ≤ The result follows by a union bound.
18
τ . 4
6.3
Proof of Proposition 4.5
Proof. We continue (10) through step K, irrespective of the stopping test. Scalar Bernstein gives kX⋆ kF ≤ τ0 ≤ 3 kX⋆ kF except with probability n−12 /48. At each reconstruction, conditional on the preceding observations and Gaussian matrices, Lemma 6.1 gives (37) except with probability n−12 /48. Conditional on the corrected matrix, (86) holds with the same failure bound and implies (38) by the proof of Lemma 6.3. Induction using Lemma 6.2 therefore gives an event, with failure probability at most (2K + 1)n−12 /48, on which (19), the sampling estimates (37), and the Gaussian estimate (86) hold throughout the complete sequence. Observed residuals. For every nonzero iterate on this event, (83) and the projection onto the true singular spaces give r µr , max{kUℓ k2,∞ , kVℓ k2,∞ } ≤ 2 n where Zℓ = Uℓ Σℓ VℓT is a compact singular value decomposition. Fix an index ℓ < K for which these preceding reconstruction estimates hold and rank(Zℓ ) = r, and condition on the observations and Gaussian matrices used to construct Zℓ . Write Eℓ = Zℓ − X⋆ . The identity Eℓ = U⋆ U⋆T Eℓ + (I − U⋆ U⋆T )Eℓ Vℓ VℓT implies
r µr kEℓ kF . kEℓ k∞ ≤ 3 n
(39)
The next observation set Ω(ℓ+1) is independent of Eℓ . For the centered sum defining Rℓ2 −kEℓ k2F , the variance and summand bounds in Lemma B.1 are at most 9µr kEℓ k4F nq
9µr kEℓ k2F , nq
and
respectively. Set W = (q −1 PΩ(ℓ+1) − I)(−Eℓ ) and T⋆ = TX⋆ Mr . Since 2µr , n
2
PT⋆ (ei eTj ) F ≤
2 the √ corresponding bounds for the vectorization of PT⋆ W are at most 2µr kEℓ kF /(nq) and 3 2µr kEℓ kF /(nq). Applying Lemma B.1 with t = 24 log n gives, after increasing C5 ,
3 5 kEℓ k2F ≤ Rℓ2 ≤ kEℓ k2F , 4 4
kPT⋆ W kF ≤
1 kEℓ kF , 16
(40)
except with conditional probability n−12 /48. These estimates also hold when Eℓ = 0. A conditional union bound over ℓ < K, together with the preceding reconstruction events, gives total failure probability at most 3K + 1 −12 n−10 n < . 48 6 We work on this joint event for the remainder of the proof. Reconstruction after rank r is attained. Fix a rank-r iterate Zℓ with ℓ < K, and write σ = σr (X⋆ ) and τ = τℓ . The preceding reconstruction retains r singular values above τℓ−1 /8 = τ /2. By Weyl’s inequality and (37), τ τ ≤σ+ , 2 128
hence
τ ≤ 3σ.
For the next corrected matrix Y = X⋆ + W , it follows that kW kop ≤
σ τ ≤ , 512 16 19
σr (Y ) ≥
15σ . 16
(41)
f = Hr (Y ) and E e=X f − X⋆ . Then E e Let X ≤ 2 kW kop . In the left and right singular op A B f e= coordinates of X⋆ , write E C D . The matrix Σ⋆ + A is invertible, and rank(X) = r gives −1 D = C(Σ⋆ + A) B. Consequently, D E 2 kW kop kW kop e ≤ e . PT⋆⊥ W , E kBkF kCkF ≤ E σ − 2 kW kop 2(σ − 2 kW kop ) F
f implies The best-approximation property of X D E 2 e ≤ 2 kPT⋆ W k E e ≤ 2 W,E e E F F
F
+
kW kop
σ − 2 kW kop
Using (41) and (40), we obtain
e E
2 F
.
1 kEℓ kF . (42) 4 F It remains to account for the finite subspace iteration in Tτ . By (41) and τ ≤ 3σ, the space range(U1 ) in the proof of Lemma 6.3 is precisely the leading r-dimensional left singular space of Y . Write Σ1 for its singular values and P = QQT for the final subspace projector. Equation (87), with t = ⌈12 log n⌉ and n32 30−t ≤ n−7 , yields τ k(I − P )U1 kop ≤ n−7 , k(I − P )U1 Σ1 kop ≤ √ n−7 . 8 2 f − X⋆ X
≤ 4 kPT⋆ W kF ≤
In particular,
p
τ . 8 Thus all r singular values are retained and Zℓ+1 = P Y has rank r. The equality of the largest principal angles between two r-dimensional spaces also gives σr (QT Y ) ≥
1 − n−14 σr (Y ) >
P (I − U1 U1T ) op = k(I − P )U1 kop ≤ n−7 . f Since Y − X
op
≤ τ /512, we have f Zℓ+1 − X
F
f + P (Y − X) f ≤ (I − P )X F F √ −7 1 1 √ + ≤ rn τ ≤ n−6 τ. 8 2 512
Together with (42), this proves 1 kEℓ kF + n−6 τℓ . (43) 4 The same argument shows that every iterate after the first rank-r iterate also has rank r. The stopping test. Therefore, whenever the residual test returns a candidate on the joint event, the compared rank-r iterates are consecutive. Write them as Zℓ and Zℓ+1 , where ℓ+1 < K and Rℓ+1 > Rℓ . By (40) and (43), r r 3 5 1 −6 kEℓ kF < kEℓ kF + n τℓ . 4 4 4 kEℓ+1 kF ≤
It follows that
σr (X⋆ ) , 1024n where the sampling condition, with C5 sufficiently large, implies n ≥ 8. Since the algorithm returns Zℓ , this proves (22). Finally, (2) gives kEℓ k♯ ≤ n kEℓ kF /(µr), and the sampling condition gives q ≥ r/n. Thus the returned matrix satisfies both local entrance conditions. This proves the proposition. kEℓ kF < 2n−6 τℓ ≤ 6n−6 σr (X⋆ ) ≤
20
6.4
Proof of Theorem 4.4
Proof of Theorem 4.4. We first consider the complete sequence in (10). Initial scale. By incoherence, |(X⋆ )ij |2 kX⋆ k2F
≤
X |(X⋆ )ij |2
µ2 r 2 , n2
i,j
kX⋆ k2F
= 1.
Scalar Bernstein applied to the observed squared entries gives ) ( kPΩ(0) (X⋆ )k2F 3n2 q n−12 1 ≤ 2 exp − , ≤ > − 1 Pr 2 28µ2 r 2 48 q kX⋆ k2F after increasing C5 , since µr ≤ n. Thus, outside this event, kX⋆ kF ≤ τ0 ≤ 3 kX⋆ kF .
(44)
In particular, τ0 > 0 and kX⋆ k♯ = σ1 (X⋆ ) ≤ τ0 . Induction step. Let Fℓ contain the observations in Ω(0) , . . . , Ω(ℓ) and the Gaussian matrices used in the first ℓ reconstructions. Conditional on Fℓ , the matrix Zℓ and the scale τℓ are fixed, whereas Ω(ℓ+1) and the next Gaussian matrix are independent. Let G0 be the event in (44), and define o n Gℓ+1 := Gℓ ∩ rank(Zℓ+1 ) ≤ r, kZℓ+1 − X⋆ k♯ ≤ 4−(ℓ+1) τ0 . On Gℓ , the induction hypothesis and Theorem 4.3 imply Pr(Gℓ \ Gℓ+1 | Fℓ ) ≤
n−12 . 24
Taking expectations and proceeding by induction proves (19) through step K on GK . Probability estimate. Since K ≤ n and G0 ⊇ · · · ⊇ GK , c Pr(GK )≤
n−12 n−12 2K + 1 −12 n−10 +K = n < . 48 24 48 6
(45)
√ Verification of the local hypotheses at step K. For K = ⌈5 + log4 (κ nr)⌉, the difference √ ZK − X⋆ has rank at most 2r. Using (19), τ0 ≤ 3 r κσr (X⋆ ), and q ≥ r/n, we obtain √ r √ √ q 3 2 r −K σr (X⋆ ) ≤ σr (X⋆ ). kZK − X⋆ kF ≤ 2r 4 τ0 ≤ 1024 n 128 For K = 6 + log4 (µr 3/2 κ) , the same invariant gives kZK − X⋆ k♯ ≤ 4−K τ0 ≤
σr (X⋆ ) 3 σr (X⋆ ) ≤ . 4096µr 1000µr
In both cases, kZK − X⋆ kop < σr (X⋆ ),
rank(ZK ) ≤ r,
so Weyl’s inequality gives rank(ZK ) = r. Residual stopping. Intersect GK with the event in Proposition 4.5. The total failure probb = K and the preceding ability is at most n−10 /3. If no residual increase is detected, then K bounds apply. Otherwise, (22), (2), and q ≥ r/n give √ q σr (X⋆ ) ZKb − X⋆ F ≤ σr (X⋆ ), ZKb − X⋆ ♯ ≤ . 128 1000µr The residual test returns only an iterate of rank r. This proves (20)–(21) in both cases. We can now apply the local RGD theorem to the initialization output. The corresponding b in Subsection 7.5. RGN argument uses the independence of Ω 21
6.5
Proof of Theorem 3.1
Proof. Choose C1 ≥ 24(C4 + C5 + 1). Since K + 1 ≤ 12 log(nκ),
p = 1 − (1 − q)K+1 ≤ (K + 1)q,
log n . It also implies K + 1 ≤ n. Hence Theorem 4.4 gives, the sampling condition gives q ≥ C1 µr 12n with probability at least 1 − n−10 /3, √ √ q p σr (X⋆ ) ≤ σr (X⋆ ). kX0 − X⋆ kF ≤ 128 128
The event in Theorem 4.1 holds simultaneously throughout this neighborhood, so no independence between X0 and Ω is required. A union bound proves the convergence statement.
7
Convergence of RGN
We prove Theorem 4.2 by controlling the initial RGN iterates in the sharp norm and then establishing a quadratic Frobenius recurrence. For the initial iterates, we compare the original sequence with auxiliary sequences obtained by completing one row or column. Once the Frobenius error is sufficiently small, Lemma 5.7 gives the required estimate uniformly, without further leave-one-out comparisons. b ∼ Bernoulli(q) and R0 := q −1 P b . The initial point X0 is Throughout this section, Ω Ω b Inverses of tangent normal operators are taken on their corresponding tangent independent of Ω. b q). spaces. Whenever Lemma 7.2 is used below, it is applied with (Λ, pΛ ) = (Ω,
7.1
Dependence and the leave-one-out construction
b so concentration cannot be applied to TX Mr as though it The difficulty is that Xk = Xk (Ω), k b Following the row- and column-deletion construction were independent of the observation set Ω. in [21, Sec. 7.2 and Algorithm 5] and [10], define, for each row i and column j, Rr,i (Z) := ei eTi Z + q −1 PΩb (I − ei eTi )Z , and set
Rc,j (Z) := Zej eTj + q −1 PΩb Z(I − ej eTj ) ,
A := {0} ∪ {(r, i) : 1 ≤ i ≤ n} ∪ {(c, j) : 1 ≤ j ≤ n}. Each completed operator replaces one row or column by its population counterpart. The corresponding auxiliary sequence is independent of the Bernoulli variables in that row or column. These sequences are used only for 0 ≤ k ≤ K̄. The following lemma gives a sampling isometry that holds simultaneously for the original operator and all completed operators. Lemma 7.1. Suppose that Assumption 2.1 holds. If q ≥ c1 µr log n/n, with c1 as in Lemma 5.6, then, with probability at least 1 − n−10 /24, max PTX⋆ Rα PTX⋆ − PTX⋆ F→F ≤ α∈A
Proof. The proof is deferred to Subsection E.1.1.
22
1 . 16
For every α ∈ A, set X0α = X0 . We use the deterministic error bound σr (X⋆ ) ρk := 280µr
7 25
2k
,
ρk+1 = 280µr
ρ2k . σr (X⋆ )
(46)
To define each auxiliary sequence on every outcome, we use the following stopping rule for 0 ≤ k < K̄. Given Xkα , let tXkα be its tangent error component from (24). If the tangent least-squares problem 1
min
ξ∈TX α Mr 2 k
hXkα + ξ − X⋆ , Rα (Xkα + ξ − X⋆ )i
has a unique minimizer, let ξkα be that minimizer and set ηkα := ξkα − tXkα . Otherwise, set ηkα := 0.
ξkα := tXkα ,
When the minimizer is unique and the graph retraction is defined, denote the candidate next state by cα := RetrX α (ξ α ). X k k+1 k We accept the update only if the minimizer is unique, the graph retraction is defined, and cα − X⋆ X k+1
kXkα − X⋆ k♯ ≤ ρk ,
♯
≤ ρk+1 ,
α cα . Otherwise, leave X α unchanged and set the subsequent := X hold. In this case, set Xk+1 k k+1 states equal to X⋆ : for every k + 1 ≤ ℓ ≤ K̄,
Xℓα := X⋆ ,
ξℓα = ηℓα := 0.
After stopping, we use the fixed compact SVD X⋆ = U⋆ Σ⋆ V⋆T . Thus the auxiliary sequence is defined on every outcome and remains independent of the variables in its completed row or column. For these stopped processes, define hcorr := max ηkα − ηk0 F . k
α 0 dloo k := max Xk − Xk F ,
α∈A
α∈A
For 0 ≤ k ≤ K̄, write Xk := Xk0 . If the updates at indices 0, . . . , k − 1 are accepted, this sequence agrees with Algorithm 2 through index k. We will prove that all these updates are accepted on the event used in the convergence proof.
7.2
Simultaneous sampling events
We first bound the number of observations in each row and column. These bounds will be used to control the sampled tangent operators. Lemma 7.2. Let Λ ∼ Bernoulli(pΛ ), and set δij := 1{(i,j)∈Λ} . Then, with probability at least 1 − 2ne−pΛ n/3 , X X δij ≤ 2pΛ n. δij ≤ 2pΛ n, max max i
j
j
i
Proof. For every row or column degree d ∼ Binomial(n, pΛ ), the multiplicative Chernoff bound [26, Corollary 5.2], with deviation parameter one, gives Pr{d > 2pΛ n} ≤ e−pΛ n/3 . Therefore, a union bound over the 2n row and column degrees proves the result. 23
The next lemma gives concentration bounds that hold simultaneously for every completed row, every completed column, and every index 0 ≤ k < K̄. Lemma 7.3. Let X0 be the initial point and write δij := 1{(i,j)∈Ω} b . For the stopped row-i and column-j trajectories, define " #T " # zkr,i := PT ⊥ (X⋆ − Xkr,i ) − ηkr,i X
r,i k
and set (wkr,i)a :=
zkc,j := PT ⊥ (X⋆ − Xkc,j ) − ηkc,j ej ,
ei ,
X
δia − 1 (zkr,i )a , q
(wkc,j )a :=
c,j k
δaj − 1 (zkc,j )a . q
Let Vkr,i and Ukc,j be the right and left singular factors selected by the fixed compact-SVD conb and satisfies vention. Suppose that Assumption 2.1 holds and X0 is independent of Ω kX0 − X⋆ k♯ ≤ ρ0 =
σr (X⋆ ) . 1000µr
Then, conditional on X0 , with probability at least 1 − n−10 /48, the following two inequalities hold simultaneously for all 1 ≤ i, j ≤ n and 0 ≤ k < K̄: s r r µr 3 log n µr r,i r,i T r,i r,i r,i 2 ≤8 + (Vk ) wk + zk 4 wk zk n q n 2 ∞ 2 2 p 64 µr/n log n r,i , (47) zk + q ∞ s r r µr 3 log n µr c,j c,j T c,j c,j c,j ≤8 + (Uk ) wk + zk 4 wk zk 2 n q n 2 ∞ 2 2 p 64 µr/n log n c,j + . (48) zk q ∞ Proof. The proof is deferred to Subsection E.2.
7.3
Error bounds for the initial RGN iterates
The next lemma bounds the sampled tangent normal operators uniformly over the sharp-norm neighborhood. It also bounds the difference between the sampled and population tangent corrections. Lemma 7.4. Let α ∈ A and X ∈ Mr , set e♯ := kX − X⋆ k♯ , and let tX be defined by (24). Suppose that Assumption 2.1 and the conclusions of Lemmas 7.1 and 7.2 hold, and that e♯ ≤ σr (X⋆ ) 1000µr . Then 11 9 kζk2F ≤ hζ, Rα ζi ≤ kζk2F , ζ ∈ TX Mr . (49) 10 10 α be the resulting unique sampled tangent minimizer and set η α := ξ α − t . Then Let ξX X X X α kηX kF ≤ 9µr
e2♯ σr (X⋆ )
.
(50)
Proof. The proof is deferred to Subsection E.3. To propagate the sharp-norm error bound, we also need coordinatewise control of the corrections. The following lemma relates this control to the difference between the original and completed corrections and then bounds that difference. 24
Lemma 7.5. Let 0 ≤ k < K̄. Suppose that Assumption 2.1 and the conclusions of Lemmas 7.1 and 7.2 hold, and that r µr α loo ρk . max kXk − X⋆ k♯ ≤ ρk , dk ≤ 2 α∈A n Then every correction satisfies r
kηkα k♯ ≤ 2
ρ2k n corr hk + 19µr , µr σr (X⋆ )
α ∈ A.
(51)
If, in addition, q ≥ c4 µr log n/n for a sufficiently large absolute constant c4 > 0 and the conclusion of Lemma 7.3 holds, then r ρ2k µr corr hk ≤ 128 µr . (52) n σr (X⋆ ) Proof. The proof is deferred to Subsection E.5. We next pass from the tangent corrections to the retracted iterates. By Lemma 5.4, if kX − X⋆ k♯ ≤ ρk and kηk♯ ≤ σr (X⋆ )/16, then kRetrX (PTX (X⋆ − X) + η) − X⋆ k♯ ≤ kηk♯ + 13
ρk kηk♯
σr (X⋆ )
+6
kηk2♯
σr (X⋆ )
.
(53)
The following lemma bounds the difference between retracted iterates at two nearby matrices. For X, Y ∈ Mr , write dX,Y := kX − Y kF . Lemma 7.6. Let X, Y ∈ Mr , let ηX ∈ TX Mr and ηY ∈ TY Mr , and use tX , tY from (24). Set dη := kηX − ηY kF . Suppose that
and
o n σr (X⋆ ) max kX − X⋆ k♯ , kY − X⋆ k♯ ≤ ρ ≤ 1000µr max {kηX kF , kηY kF } ≤ 9µr
ρ2 . σr (X⋆ )
Then RetrX (tX + ηX ) and RetrY (tY + ηY ) are well defined, and 5 µrρ2 kRetrX (tX + ηX ) − RetrY (tY + ηY )kF ≤ dη + 19 dX,Y . 4 σr (X⋆ )2
(54)
Proof. The proof is deferred to Subsection C.4.2.
7.4
Proof of Theorem 4.2
Proof. We prove by induction that max kXkα − X⋆ k♯ ≤ ρk , α∈A
r
dloo k ≤2
µr ρk . n
Whenever the candidate retraction is defined, the update rule gives i h cα − X⋆ = η α + RetrX α (tX α + η α ) − X⋆ − η α . X k+1 k k k k k
For the actual and an auxiliary candidate, subtracting the two updates gives h i cα − X⋆ − η α . c0 − X⋆ − η 0 − X cα = η 0 − η α + X c0 − X X k k+1 k k k k+1 k+1 k+1 25
(55)
(56)
Step 1: bounds for the initial RGN iterates. Condition on an arbitrary realization of X0 satisb this conditioning leaves fying kX0 − X⋆ k♯ ≤ σr (X⋆ )/(1000µr). Since X0 is independent of Ω, b unchanged. Fix C4 ≥ c1 + c4 + 128. By the sampling assumption, the law of Ω q ≥ C4
µr log n µr log n ≥ (c1 + c4 ) . n n
Thus the sampling-rate hypotheses of Lemmas 7.1 and 7.5 are satisfied. Suppose that the conclusions of Lemmas 7.1 to 7.3 hold simultaneously. At k = 0, all auxiliary trajectories equal X0 and dloo 0 = 0, so (55) holds at k = 0. Suppose inductively that transitions 0, . . . , k − 1 have been accepted for every auxiliary process and that (55) holds at some k < K̄. Since ρk ≤ σr (X⋆ )/(1000µr), Lemma 7.4 implies that every sampled tangent normal operator at iteration k is invertible. Hence every tangent minimizer is unique, and its correction satisfies PTX α Rα PTX α ηkα = PTX α Rα NXkα , α ∈ A. k
k
k
Combining (51) and (52), and using ρk ≤ σr (X⋆ )/(1000µr), gives kηkα k♯ ≤ 275µr
ρ2k 11 1 ≤ ρk < σr (X⋆ ), σr (X⋆ ) 40 16
α ∈ A.
(57)
The bound (55) at index k and (57) verify the hypotheses of Lemma 5.4 with X = Xkα and η = ηkα . Hence every candidate graph core is invertible and every candidate retraction is well defined. Applying (53) gives cα − X⋆ X k+1
♯
< 280µr
ρ2k = ρk+1 . σr (X⋆ )
Thus transition k is accepted for every α ∈ A. For the two-base update, Lemma 7.4 and the induction hypothesis give max{ ηk0 F , kηkα kF } ≤ 9µrρ2k /σr (X⋆ ), which verifies the correction-size hypothesis of Lemma 7.6. Using (54), (52), and the induction hypothesis, we obtain r µrρ2k loo µr 5 corr α 0 Xk+1 − Xk+1 F ≤ hk + 19 dk < 2 ρk+1 . 2 4 σr (X⋆ ) n Taking the maximum over α proves (55) at index k + 1. Therefore, all updates before K̄ are accepted, and r µr loo α ρk , 0 ≤ k ≤ K̄. dk ≤ 2 max kXk − X⋆ k♯ ≤ ρk , α∈A n In particular, kXk − X⋆ k♯ ≤ ρk for 0 ≤ k ≤ K̄. Step 2: the Frobenius error at K̄. Since the difference of two rank-r matrices has rank at √ K̄ most 2r, 2−2 ≤ (4n)−1 , 2r ≤ 2µr, and q ≥ C4 µr log n/n ≥ n−2 , we have √ √ kXK̄ − X⋆ kF ≤ 2r kXK̄ − X⋆ k♯ ≤ 2r ρK̄ √ 2r −2K̄ 1 1√ q σr (X⋆ ). 2 σr (X⋆ ) ≤ σr (X⋆ ) ≤ ≤ 280µr 40n 40 Since all updates before K̄ are accepted, the stopped sequence agrees with Algorithm 2 through index K̄. From this point onward, Xk denotes the iterates of that algorithm without the auxiliary stopping rule. Step 3: quadratic convergence in the Frobenius norm. Since 0 ∈ A, Lemma 7.1 gives PTX⋆ q −1 PΩb PTX⋆ − PTX⋆ F→F ≤ 26
1 . 16
Fix k ≥ K̄ and suppose that ek := kXk − X⋆ kF ≤
1√ q σr (X⋆ ). 40
Whenever the tangent minimizer is unique, set ηk := ξk − tXk . The update rule gives Xk+1 − X⋆ = ηk + [RetrXk (tXk + ηk ) − X⋆ − ηk ] .
(58)
b q), the sampled tangent normal operator on TX Mr is invertBy Lemma 5.7 with (Λ, pΛ ) = (Ω, k ible and r −1 3 5 −1 −1/2 ≤ , PTXk q PΩb PTXk q . (59) PΩb PTXk ≤ 2 4 F→F F→F Hence the tangent minimizer ξk is unique, and
PTXk q −1 PΩb PTXk ηk = PTXk q −1 PΩb NXk . The self-adjointness of PΩb , (59), part (i) of Lemma 5.5, and ek ≤
√ q σr (X⋆ )/40 give
3 PT q −1 PΩb NXk 2 r Xk F 3 5 ≤ kNXk kF 2 4q √ e2 3 5 1 σr (X⋆ ). ≤ < √ k 2 q σr (X⋆ ) 320
kηk kF ≤
(60)
The same neighborhood gives ek ≤ σr (X⋆ )/40. Combining (58), (60), and part (ii) of √ Lemma 5.5, and using ek ≤ q σr (X⋆ )/40, shows that the retraction is well defined and yields 11 ek 11 kηk kF kXk+1 − X⋆ kF ≤ kηk kF 1 + + 5 σr (X⋆ ) 10 σr (X⋆ ) e2 1 1√ ≤ 4√ k q σr (X⋆ ). ≤ ek ≤ q σr (X⋆ ) 10 40 Starting from k = K̄, induction proves the Q-quadratic recurrence for every k ≥ K̄. It follows that Xk → X⋆ . Conditional on the chosen realization of X0 , a union bound for the events in Lemmas 7.1 to 7.3, together with q ≥ C4 µr log n/n with C4 ≥ 128 and n ≥ 2, gives total failure probability at most n−10 n−10 n−10 + 2ne−qn/3 + < . (61) 24 48 12 The argument applies to every realization of X0 satisfying kX0 − X⋆ k♯ ≤ σr (X⋆ )/(1000µr). Therefore, (61) gives the conditional probability asserted in Theorem 4.2 and completes the proof.
7.5
Proof of Theorem 3.2
Proof. Choose C2 ≥ 32(C4 + C5 + 1). Since K + 2 ≤ 16 log(2µrκ), the sampling condition gives q≥
p = 1 − (1 − q)K+2 ≤ (K + 2)q, C2 µr log n . 16n 27
It also implies K + 2 ≤ n. Thus Theorem 4.4 gives, with failure probability at most n−10 /3, kX0 − X⋆ k♯ ≤
σr (X⋆ ) . 1000µr
The stopping index and the returned initializer depend only on Ω(0) , . . . , Ω(K) and the indeb Conditional on a successful pendent Gaussian matrices, and are therefore independent of Ω. −10 initialization, Theorem 4.2 fails with probability at most n /12. Hence the union bound gives −10 −10 the convergence statement in Theorem 3.2 fails with probability at most n 3 + n12 < n−10 .
7.6
Proof of Theorem 3.3
Proof. The proof is based on the event in the proof of Theorem 3.2. By Step 2 in the proof of (X⋆ ) Theorem 4.2, we have kXK̄ − X⋆ kF ≤ σr40n . For C2 sufficiently large, the sampling condition √
q
also gives kXK̄ − X⋆ kF ≤ 128 σr (X⋆ ). Fix k ≥ K̄ and suppose that ek := kXk − X⋆ kF ≤
√
q σr (X⋆ ). 128
Let ξbk denote the ordinary RGN direction at Xk and set bk := ξbk − tXk . η
The argument in Step 3 of the proof of Theorem 4.2, together with (29), gives √ e2 17 3 5 , ξbk ≤ ek . kηbk kF ≤ √ k 2 q σr (X⋆ ) 16 F
Moreover, by the definition of λk , the ordinary normal equation, (29), and Weyl’s inequality, PTXk q −1 PΩb PTXk ξbk λk 5 ek F = ≤ . q σr (Xk ) 4 σr (X⋆ )
Subtracting the ordinary and regularized normal equations gives λk λk −1 PTXk q PΩb PTXk + I (ξk − ξbk ) = − ξbk . q q
Again by (29), we have
ξk − ξbk
F
≤
85 e2k . 56 σr (X⋆ )
F
< 5√
Hence, with ηk := ξk − tXk , we have
bk kF + ξk − ξbk kηk kF ≤ kη
e2k σr (X⋆ ) < . q σr (X⋆ ) 320
Part (ii) of Lemma 5.5 therefore applies and yields e2 11 kηk kF 11 ek + ≤ C̄2 √ k kXk+1 − X⋆ kF ≤ kηk kF 1 + 5 σr (X⋆ ) 10 σr (X⋆ ) q σr (X⋆ ) √ for an absolute constant C̄2 > 0. The right-hand side is at most q σr (X⋆ )/128 for C2 sufficiently large. Thus the stated neighborhood is invariant, and induction proves the result for every k ≥ K̄. 28
8
Numerical Experiments
In this section, we compare factorized gradient descent (FGD), ScaledGD [25], RGD, and RGN on randomly generated matrix completion problems. All computations were performed in MATLAB R2024b on 64-bit Windows, with an Intel Core Ultra 5 125H CPU and 16 GB memory. We measure the reconstruction error by errk :=
kXk − X⋆ kF . kX⋆ kF
Convergence and running time. We first compare the convergence and running time under different condition numbers. We take n = 1000, r = 10, and p = 0.2. Let U⋆ contain the left singular vectors of an n × r matrix with independent random signs. For κ ∈ {1, 10, 100}, we set X⋆ = U⋆ diag(σ1 , . . . , σr )U⋆T , where the nonzero singular values are linearly spaced from 1 to 1/κ. Each entry is observed independently with probability p without noise. We use the same U⋆ and Ω for all three values of κ; for each κ, the four methods receive the same observations PΩ (X⋆ ). FGD and ScaledGD use the same spectral initialization and the respective updates in [25], with stepsize 0.5 for both methods. RGD and RGN share the initialization in Algorithm 4, with q = p, Ω(ℓ) = Ω, and K = 8, returning the rank-r candidate with the smallest observed residual. In (11), we additionally stop when 2 1/2 1 − r −1 QTℓ,t−1 Qℓ,t F ≤ 10−3 holds at two consecutive iterations with t ≥ 3; the iteration limit remains ⌈12 log n⌉. RGD and b = Ω. For RGN, CG is stopped at relative RGN follow Algorithms 1 and 3, respectively, with Ω residual 0.05 or after 200 iterations. All four methods terminate when kPΩ (Xk − X⋆ )kF ≤ 10−14 , kPΩ (X⋆ )kF
in at most 1000 iterations. Figure 2 plots errk against the iteration count and elapsed time for one realization at each κ. Elapsed time includes initialization and the subsequent iterations. Empirical recovery rates. We next examine how the number of observations needed for recovery varies with the rank. We set n = 500 and conduct 30 random trials for each pair r ∈ {2, 4, . . . , 30}, m/n ∈ {4, 8, . . . , 84}. In each trial, let U⋆,r and V⋆,r consist of the first r T . We columns of two independent n × 30 Haar orthonormal frames, and form X⋆ = U⋆,r V⋆,r observe m entries uniformly without replacement. For this experiment, the search directions are b = Ω. The relative observed residual tolerances are 10−12 those in Algorithms 1 and 2, with Ω for RGD and 10−8 for RGN. Recovery is declared successful if errk ≤ 10−6 . Figure 3 reports the fraction of successful trials at each (r, m). Both panels exhibit a transition from failure to successful recovery as m increases. These results show an approximately linear dependence on nr over the tested ranks. The solid lines in Figure 3 provide simple reference relations for this observed transition.
9
Concluding Remarks
We have established global recovery guarantees for RGD and RGN under standard incoherence and Bernoulli sampling. With high probability, O(µnr log n log(nκ)) and O(µnr log n log(2µrκ)) observations suffice for the two methods, respectively. The sample requirements are linear in n and r, up to logarithmic factors. RGD converges linearly, while RGN satisfies a doubly 29
relative Frobenius error relative Frobenius error relative Frobenius error
100 FGD ScaledGD RGD RGN
10-5
FGD ScaledGD RGD RGN
10-10 10-15
0
100 200 iteration
300
0
1
2 3 elapsed time (s)
4
100 10-5 10-10 10-15
FGD ScaledGD RGD RGN
0
100 200 iteration
FGD ScaledGD RGD RGN
300
0
5 10 elapsed time (s)
100 10-5 10-10 10-15
FGD ScaledGD RGD RGN
0
100 200 iteration
FGD ScaledGD RGD RGN
300
0
2
4 6 8 elapsed time (s)
10
Figure 2: Relative reconstruction errors of FGD, ScaledGD, RGD, and RGN versus iteration count (left) and elapsed time in seconds (right), with n = 1000, r = 10, and p = 0.2. From top to bottom, κ = 1, 10, 100. exponential sharp-norm error bound during its initial iterations and a quadratic Frobenius recurrence thereafter. Both methods use the multiscale residual initialization with residual stopping and different upper bounds on the number of reconstructions for their respective local conditions. The observation set is fixed at the beginning, and the sampling requirement does not depend on the target accuracy. Nevertheless, several extensions merit further study. One is to establish the same guarantees when every initialization step reuses the entire observation set. Another is to analyze RGN with an inexact tangent solve and a computable stopping criterion, so that the convergence rate can be related to the total number of CG iterations.
30
(b) RGD
1
30 28 26 24 22 20 18 16 14 12 10 8 6 4 2
4 12 20 28 36 44 52 60 68 76 84 Ratio m/ n
0.8
0.6
0.4
Empirical success rate
30 28 26 24 22 20 18 16 14 12 10 8 6 4 2
Rank r
Rank r
(a) RGN
0.2
0 4 12 20 28 36 44 52 60 68 76 84 Ratio m/ n
Figure 3: Empirical recovery rates of RGN (left) and RGD (middle) for n = 500 and κ = 1, over 30 trials at each (r, m). The horizontal axis is m/n and the vertical axis is r. White denotes success in all trials, and black denotes failure in all trials. The solid reference lines are m = (2r + 9)n for RGN and m = (2.5r + 10)n for RGD.
Declaration on the use of artificial intelligence Generative AI tools were used during revision to assist with language editing, structural reorganization, bibliographic cross-checking, and consistency checks of notation, formulas, and proofs.
A
Lower bound for ordinary spectral initialization
We prove Theorem 3.4 by constructing a block-diagonal family of incoherent matrices. Proof of Theorem 3.4. Take C3 = 2−12 . Choose pairwise disjoint sets I1 , . . . , Ir ⊆ {1, . . . , n} with |Ia | = s, and set ua = v a = s
−1/2
X⋆ = κu1 v1T +
1Ia ,
r X
ua vaT .
a=2
Then σr (X⋆ ) = 1,
2
2
max eTi U⋆ 2 = max eTj V⋆ 2 =
κ(X⋆ ) = κ,
j
i
µr 1 = . s n
Thus Assumption 2.1 holds. Set d = qs and Y = q −1 PΛ (X⋆ ). On the active coordinates, Y = diag(Y1 , . . . , Yr ),
Y1 =
κ ∆1 , d
Ya =
1 ∆a , d
2 ≤ a ≤ r,
where the ∆a ∈ {0, 1}s×s have independent Bernoulli(q) entries. Since µ = n/(rs) and κ2 ≤ s, d = qs < C3 κ2 ,
q < C3
31
1 κ2 ≤ C3 < . s 2
Suppose first that d ≤ 14 log(rs). For each active row,
π := Pr{the row is empty on its support} = (1 − q)s ≥ e−2qs = e−2d ≥ (rs)−1/2 .
Hence √ 1 Pr{an active row is empty} = 1 − (1 − π)rs ≥ 1 − e−rsπ ≥ 1 − e− rs > . 2 If i ∈ Ia is such a row, then
eTi Y = 0
=⇒
eTi Hr (Y ) = 0,
eTi X⋆ 2 ≥ s−1/2 ,
and therefore
1 1√ s eTi (Hr (Y ) − X⋆ ) 2 ≥ . 2 2 Now suppose that d > 14 log(rs). Let D be the event that every row and column degree of every ∆a is at most 16d. The multiplicative Chernoff bound gives 1 Pr(D c ) ≤ 2rs exp −(16 log 16 − 15)d < 2(rs)1−(16 log 16−15)/4 ≤ . 30 kHr (Y ) − X⋆ k♯ ≥
With M = k∆1 k2F , we have
EM 2 ≤ s2 d2 + sd,
EM = sd,
sd >
log2 (rs) d2 > . C3 16C3
Hence the Paley–Zygmund inequality yields 1 9 sd − Pr(D c ) > . 16 sd + 1 2 On D, the maximum row sum and the maximum column sum of ∆a are both at most 16d. Hence v u ! u X X u k∆a kop ≤ t max |(∆a )ij | ≤ 16d, 1 ≤ a ≤ r, |(∆a )ij | max Pr{D ∩ {M ≥ sd/4}} ≥
i
j
j
i
and therefore
kY1 k2F ≥
κ2 s , 4d
max kYa kop ≤ 16,
2≤a≤r
Thus
X j≥2
kY1 kop ≤ 16κ, σj (Y1 )2 ≥
κ2 s − 256κ2 > 4d
1 − 256 s > 256(s − 1). 4C3
σ2 (Y1 ) > 16 ≥ max kYa kop . 2≤a≤r
If Y1 has at least r singular values larger than 16, then every best rank-at-most-r approximation of Y is supported on I1 × I1 , and hence kHr (Y ) − X⋆ k♯ ≥ P(I2 ∪···∪Ir )×(I2 ∪···∪Ir ) X⋆ op = 1. Otherwise, the Eckart–Young–Mirsky theorem and σ2 (Y1 ) > 16 imply, with I = I2 ∪ · · · ∪ Ir , rank(PI×I Hr (Y )) ≤ r − 2. Since PI×I X⋆ has r − 1 singular values equal to one, kHr (Y ) − X⋆ k♯ ≥ kPI×I (Hr (Y ) − X⋆ )kop ≥ 1. In both cases,
This proves the theorem.
1 1 > . Pr kHr (Y ) − X⋆ k♯ ≥ 2 2
32
B
Probabilistic Lemmas
We collect the concentration tools used throughout the sampling arguments. We use the following standard forms of the rectangular and self-adjoint matrix Bernstein inequalities; see [26, Theorems 1.6 and 6.1]. Lemma B.1. Let Z1 , . . . , ZN ∈ Rd1 ×d2 be independent mean-zero random matrices. Suppose that kZa kop ≤ L almost surely, and define X X . E(ZaT Za ) , v := max E(Za ZaT ) a a op
op
Then, for every t > 0,
Pr
X
Za
a
op
√ 2 > 2vt + Lt ≤ (d1 + d2 )e−t . 3
If the Za are self-adjoint operators on a d-dimensional Hilbert space, then, for every u > 0, X u2 Pr Za ≥ u ≤ 2d exp − . a 2(v + Lu/3) op
B.1
Proof of Lemma 5.6
We next prove the true-tangent sampling isometry used in the local RGD analysis. Proof of Lemma 5.6. Fix c1 = 215 . Let δij := 1{(i,j)∈Λ} and aij = PTX⋆ (ei eTj ), and define (a ⊗ a)(Z) := ha, ZiF a. Since the matrices ei eTj form an orthonormal basis of Rn×n , we have X aij ⊗ aij = ITX⋆ . i,j
Consequently, we have PTX⋆ p−1 Λ PΛ PTX⋆ − PTX⋆ =
X δij i,j
pΛ
− 1 (aij ⊗ aij ).
By (4) and the orthogonality of its two summands, we have aij = U⋆ U⋆T ei eTj + (I − U⋆ U⋆T )ei eTj V⋆ V⋆T , 2
2
2
kaij k2F = U⋆ U⋆T ei 2 + (I − U⋆ U⋆T )ei 2 V⋆ V⋆T ej 2 ≤
2µr . n
Each summand has norm at most 2µr/(npΛ ). Moreover, E(δij /pΛ − 1)2 = (1 − pΛ )/pΛ and (a ⊗ a)2 = kak2F (a ⊗ a), so the variance operator satisfies 2µr 1 − pΛ X kaij k2F (aij ⊗ aij ) IT . pΛ npΛ X⋆ i,j
Since dim(TX⋆ ) = r(2n − r) ≤ 2nr, pΛ ≥ c1 µr log n/n, r ≤ n, and n ≥ 2, Lemma B.1 at threshold 1/16 gives 1 n−10 −1 Pr PTX⋆ pΛ PΛ PTX⋆ − PTX⋆ F→F > . ≤ 16 48 This proves the lemma. 33
B.2
Proof of Lemma 6.1
We prove the sampling estimates for a fixed error matrix. Independence from the current observation component is supplied by conditioning in Subsection 6.4. Proof of Lemma 6.1. By (2), 4µr τ, kEk∞ ≤ n
max{kEk2,∞ , E
Hence, for W = (q −1 PΛ − I)E, |Wij | ≤ By Lemma B.1,
For |z| = τ /64, set
4µrτ , nq
max
max i
X
T
EWij2 , max j
j
r µr τ. }≤2 2,∞ n
X
EWij2
i
n nq τ o ≤ 4n exp −c . Pr kW kop > 512 µr
≤
4µrτ 2 . nq
(62)
(63)
F (z) = (zI − W )−1 Q⋆ .
For i ∈ [2n], let W (i) be the principal minor obtained by deleting row and column i, let wi be the deleted off-diagonal column, and define ( (zI − W (i) )−1 (Q⋆ )−i,: , W (i) op ≤ τ /512, F (i) (z) = 0, otherwise. Conditional on W (i) , (62) and Lemma B.1, applied to the real and imaginary parts, give nq τ (i) (i) T (i) . (64) F (z) W ≤ (4r + 1) exp −c Pr wi F (z) > 1024 µr 2,∞ 2 Let Γ be a uniform grid on |z| = τ /64 with √ |Γ| ≤ 52 n,
τ min |z − z0 | ≤ z0 ∈Γ 1024
r
µr . n
A union bound in (63)–(64) yields, except with probability Cn3 exp(−cnq/(µr)), kW kop ≤
τ , 512
wiT F (i) (z)
2
≤
τ F (i) (z) 1024 2,∞
(65)
simultaneously for i ∈ [2n] and z ∈ Γ. Increasing c2 makes this failure probability at most n−12 /48. Fix a realization in (65). Then (zI − W (i) )−1
op
≤
512 , 7τ
kwi k2 ≤
τ , 512
and F−i,: (z) = F (i) (z) + (zI − W (i) )−1 wi Fi,: (z), 8 F (i) (z) ≤ kF (z)k2,∞ , 7 2,∞ (Q⋆ )i,: + wiT F (i) (z) Fi,: (z) = . z − wiT (zI − W (i) )−1 wi 34
Moreover, z − wiT (zI − W (i) )−1 wi ≥
τ 55τ > . 3584 66
Using Assumption 2.1 and (65), r r τ 72 µr µr 66 + kF (z)k2,∞ ≤ , kF (z)k2,∞ ≤ τ n 896 τ n
z ∈ Γ.
For arbitrary |z| = τ /64, choose z0 ∈ Γ as above. The resolvent identity gives kF (z) − F (z0 )kop ≤
512 7τ
and hence −1
sup
(zI − W )
|z|=τ /64
Finally, W j Q⋆ =
1 2πi
I
|z|=τ /64
2
6 |z − z0 | ≤ τ
128 Q⋆ 2,∞ ≤ τ
r
r
µr , n
µr . n
z j (zI − W )−1 Q⋆ dz,
(66)
j ≥ 0.
Combining this identity with (66) proves (37).
C
Deterministic geometry of fixed-rank matrices
We prove the geometric estimates used in Sections 5 to 7. We first derive the graph-retraction identity and the sharp-norm bounds at one matrix. We then prove the Frobenius estimates and the comparison bounds for two nearby matrices.
C.1
Tangent coordinates and the graph-retraction identity
Let X = U ΣV T ∈ Mr and ξ ∈ TX Mr . Define the orthogonal tangent coordinates M := U T ξV ,
B := (I − U U T )ξV ,
C := (I − V V T )ξ T U .
Then U T B = 0, V T C = 0, and ξ = U M V T + BV T + U C T . Moreover, U T (X + ξ)V = Σ + M . Hence, whenever Σ + M is nonsingular, (5) gives T RetrX (ξ) = U + B(Σ + M )−1 (Σ + M ) V + C(Σ + M )−T .
Expanding the product yields
RetrX (ξ) = X + ξ + B(Σ + M )−1 C T .
(67)
For later two-base comparisons, we also record the invariance under orthogonal changes of basis. For a compact orthonormal factorization X = U SV T , orthogonal changes of basis U 7→ U OU and V 7→ V OV leave X, PTX , and the graph retraction unchanged when the core and tangent coordinates are transformed accordingly. Thus, when comparing two matrices, we align their left and right bases separately by orthogonal Procrustes transformations and transform their cores at the same time. The estimates below are invariant under these choices; in particular, σmin (S) = σr (X).
35
C.2
Sharp-norm estimates
We begin by controlling the current singular spaces and their projectors in the sharp-norm neighborhood. C.2.1
Proof of Lemma 5.1
Proof. Let E := X⋆ −X. Weyl’s inequality [5, Supplement, Theorem 2] gives σr (X) ≥ σr (X⋆ )− e♯ . Since (I − U⋆ U⋆T )X⋆ = 0 and (I − V⋆ V⋆T )X⋆T = 0, we have (I − U⋆ U⋆T )U = −(I − U⋆ U⋆T )EV S −1 , (I − V⋆ V⋆T )V
(68a)
= −(I − V⋆ V⋆T )E T U S −T .
(68b)
For every i and j, Assumption 2.1 gives max
n
eTi U⋆ U⋆T U 2 ,
eTj V⋆ V⋆T V
2
o
≤
r
µr . n
Equations (68a) and (68b), together with e♯ ≤ σr (X⋆ )/8, give
p 3 µr/n e♯ ≤ , σr (X⋆ ) − e♯ σr (X⋆ ) − e♯ r r 3e♯ 10 µr µr T , 1+ ≤ ei U 2 ≤ n σr (X⋆ ) − e♯ 7 n p eTj E T + V⋆T ej 2 V⋆T E T op 3 µr/n e♯ T T 2 ej (I − V⋆ V⋆ )V 2 ≤ ≤ , σr (X⋆ ) − e♯ σr (X⋆ ) − e♯ r r 3e♯ 10 µr µr T , ej V 2 ≤ 1+ ≤ n σr (X⋆ ) − e♯ 7 n r r µr 10 µr <2 . max{kU k2,∞ , kV k2,∞ } ≤ 7 n n eTi (I − U⋆ U⋆T )U 2 ≤
eTi E 2 + U⋆T ei 2 U⋆T E op
For two orthogonal projectors of the same rank, the operator norm of their difference equals the largest sine of the principal angles. Hence (68a), (68b), and e♯ ≤ σr (X⋆ )/8 give e♯ , σr (X⋆ ) − e♯ e♯ V V T − V⋆ V⋆T op = (I − V⋆ V⋆T )V op ≤ , σr (X⋆ ) − e♯ n o 8 e ♯ . max U U T − U⋆ U⋆T op , V V T − V⋆ V⋆T op ≤ 7 σr (X⋆ ) U U T − U⋆ U⋆T
op
= (I − U⋆ U⋆T )U op ≤
For the rowwise projector bounds, the twop projector identities, e♯ ≤ σr (X⋆ )/8, and max{kU k2,∞ , kV k2,∞ } ≤ (10/7) µr/n give
Assumption 2.1,
U U T − U⋆ U⋆T = (I − U⋆ U⋆T )U U T − U⋆ U⋆T (I − U U T ), p r 4 µr/n e♯ 32 µr e♯ T T U U − U⋆ U⋆ 2,∞ ≤ ≤ , σr (X⋆ ) − e♯ 7 n σr (X⋆ )
V V T − V⋆ V⋆T = (I − V⋆ V⋆T )V V T − V⋆ V⋆T (I − V V T ), p r 4 µr/n e♯ 32 µr e♯ T T ≤ , V V − V⋆ V⋆ 2,∞ ≤ σr (X⋆ ) − e♯ 7 n σr (X⋆ ) r n o 32 µr e ♯ max U U T − U⋆ U⋆T 2,∞ , V V T − V⋆ V⋆T 2,∞ ≤ . 7 n σr (X⋆ ) 36
(69)
Finally, the tangent-projector identity and (69) give (PTX − PTX⋆ )Z = (U U T − U⋆ U⋆T )Z(I − V V T ) + (I − U⋆ U⋆T )Z(V V T − V⋆ V⋆T ), 16 e♯ PTX − PTX⋆ F→F ≤ . 7 σr (X⋆ ) This proves the lemma. We next control the normal component and the associated graph factors. C.2.2
Proof of Lemma 5.2
Proof. Let E := X⋆ − X. Weyl’s inequality and kGX − Skop ≤ e♯ give 3 σmin (GX ) ≥ σr (X⋆ ) − 2e♯ ≥ σr (X⋆ ), 4
G−1 X op ≤
4 . 3σr (X⋆ )
Relative to the orthogonal decompositions generated by U and V , the upper-left, lower-left, and upper-right blocks of X⋆ are GX , LX , and RTX , respectively. Since rank(X⋆ ) = r and GX is invertible, the Schur complement of GX vanishes. Therefore, the normal component is exactly T (70) NX = LX G−1 X RX . Moreover, we have max{kLX kop , kRX kop } ≤ e♯ ,
kNX kop ≤
2 4 e♯ . 3 σr (X⋆ )
For the rowwise graph factors, the definition of the sharp norm and the incoherence bound in Lemma 5.1 give, for every i, we have r µr T T T T e♯ , ei LX 2 ≤ ei E 2 + ei U 2 U E op ≤ 4 n r µr T T T T T T (71) e♯ , ei RX 2 ≤ ei E 2 + ei V 2 V E op ≤ 4 n r µr max{kLX k2,∞ , kRX k2,∞ } ≤ 4 e♯ . n Applying (71) to the left and right factors in (70) yields n o 16 r µr e2 ♯ T max kNX k2,∞ , NX 2,∞ ≤ . 3 n σr (X⋆ )
For every i, j, using the rowwise bounds on both graph factors gives T eTi NX ej ≤ eTi LX 2 G−1 X op ej RX 2 ≤
kNX k∞ ≤
2 64µr e♯ , 3n σr (X⋆ )
2 64µr e♯ . 3n σr (X⋆ )
T , and kNX k∞ gives Applying (2) to the bounds for kNX kop , kNX k2,∞ , NX 2,∞ 2 16 e♯ , kNX k♯ ≤ 3 σr (X⋆ )
which proves the sharp-norm bound in Lemma 5.2. We next use the graph factorization to verify exactness of the population tangent correction. 37
C.2.3
Proof of Lemma 5.3
Proof. Equation (4) gives the tangent coordinates M = GX − Σ,
B = LX ,
C = RX .
Moreover, Weyl’s inequality gives σmin (GX ) ≥ σr (X⋆ ) − 2 kX − X⋆ k♯ ≥
1 σr (X⋆ ) > 0, 2
so the updated core is invertible. Since rank(X⋆ ) = r, the Schur complement of GX in the block representation of X⋆ vanishes. Using tX = PTX (X⋆ − X) in (67) gives T RetrX (tX ) = X + tX + LX G−1 X RX , T PT ⊥ (X⋆ − X) = LX G−1 X RX , X
RetrX (tX ) = X + PTX (X⋆ − X) + PT ⊥ (X⋆ − X) = X⋆ . X
This proves the population identity. For the proofs of Lemmas 5.4 and 5.5, let X = U ΣV T ∈ Mr and η ∈ TX Mr , and define Mη := U T ηV ,
Bη := (I − U U T )ηV ,
Cη := (I − V V T )η T U ,
Sη := GX + Mη . (72) These coordinates describe the perturbation of the population tangent correction. C.2.4
Proof of Lemma 5.4
Proof. By the definition of the sharp norm and (72), we have n o max kMη kop , kBη kop , kCη kop ≤ kηk♯ , r n o µr kηk♯ . max kBη k2,∞ , kCη k2,∞ ≤ 4 n
By Lemma 5.2, we have G−1 X op ≤ 4/(3σr (X⋆ )). Since σmin (GX ) ≥ 3σr (X⋆ )/4 and kηk♯ ≤ σr (X⋆ )/16, we have 16 . (GX + Mη )−1 op ≤ 11σr (X⋆ )
Since e♯ ≤ σr (X⋆ )/8 < σr (X⋆ )/4, Lemma 5.3 applies. It follows from Lemma 5.3 and (67) that RetrX (tX + η) − X⋆ − η
−1 T = Bη Sη−1 RTX + LX Sη−1 CηT + Bη Sη−1 CηT − LX G−1 X Mη Sη RX .
For F , H ∈ Rn×r and K ∈ Rr×r , we have r n 1 n T F KH ♯ ≤ kKkop max kF kop kHkop , kF k2,∞ kHkop , 2 µr r o n 1 n kF kop kHk2,∞ , kF k2,∞ kHk2,∞ . 2 µr 4µr
(73)
(74)
Applying (74) to the four terms in (73), and using (26), e♯ ≤ σr (X⋆ )/8, and kηk♯ ≤ σr (X⋆ )/16, gives e♯ kηk♯ kηk2♯ kRetrX (tX + η) − X⋆ − ηk♯ ≤ 13 +6 . σr (X⋆ ) σr (X⋆ ) This proves the nonlinear remainder bound. 38
C.3
Proof of Lemma 5.5
We next prove the Frobenius estimates used in the local RGD argument and in the quadratic RGN continuation. Proof of Lemma 5.5. We first prove part (i). Let X = U ΣV T be a compact singular value decomposition, and use the graph factors in (23). Weyl’s inequality and kGX − Σkop ≤ kX − X⋆ kF give 1 σmin (GX ) ≥ σr (X⋆ ) − 2 kX − X⋆ kF ≥ σr (X⋆ ), 2
G−1 X op ≤
2 . σr (X⋆ )
Since rank(X⋆ ) = r and GX is invertible, the Schur complement of GX in the block representation of X⋆ vanishes. Moreover, kLX kF ≤ kX − X⋆ kF and kRX kop ≤ kX − X⋆ kF . Therefore, we have T PT ⊥ (X⋆ − X) = LX G−1 X RX , X
kX − X⋆ k2F PT ⊥ (X⋆ − X) ≤ 2 . X σr (X⋆ ) F
For the tangent-space motion, Weyl’s inequality gives σr (X) ≥ σr (X⋆ ) − kX − X⋆ kF . Hence the equal-rank projector identity gives U U T − U⋆ U⋆T V V T − V⋆ V⋆T
kX − X⋆ kF 4 kX − X⋆ kF ≤ , σr (X⋆ ) − kX − X⋆ kF 3σr (X⋆ ) 4 kX − X⋆ kF kX − X⋆ kF ≤ . ≤ op σr (X⋆ ) − kX − X⋆ kF 3σr (X⋆ )
op
≤
Using the tangent-projector identity, we consequently obtain PTX − PTX⋆ F→F ≤
8 kX − X⋆ kF kX − X⋆ kF ≤3 . 3σr (X⋆ ) σr (X⋆ )
This proves part (i). We next prove part (ii). Let X = U ΣV T be a compact singular value decomposition, and use the graph factors in (23). Weyl’s inequality and the tangent-coordinate bound kMη kop ≤ kηkF give G−1 X op ≤
1 , σr (X⋆ ) − 2 kX − X⋆ kF
Sη−1 op ≤
1 . σr (X⋆ ) − 2 kX − X⋆ kF − kηkF
In particular, the two cores are invertible under the hypotheses. Since rank(X⋆ ) = r, the Schur T complement identity gives PT ⊥ (X⋆ − X) = LX G−1 X RX . Applying (67) to tX + η and using X
−1 −1 Sη−1 − G−1 X = −GX Mη Sη , we have
RetrX (tX + η) − X⋆ − η
−1 Mη Sη−1 RTX . = Bη Sη−1 RTX + LX Sη−1 CηT + Bη Sη−1 CηT − LX GX
(75)
Moreover, max{kLX kF , kRX kF } ≤ kX − X⋆ kF and max{kMη kF , kBη kF , kCη kF } ≤ kηkF . Combining these bounds with (75) gives kRetrX (tX + η) − X⋆ − ηkF ≤
2 kX − X⋆ kF kηkF + kηk2F σr (X⋆ ) − 2 kX − X⋆ kF − kηkF +
≤
kX − X⋆ k2F kηkF σr (X⋆ ) − 2 kX − X⋆ kF σr (X⋆ ) − 2 kX − X⋆ kF − kηkF
11 kX − X⋆ kF kηkF 11 kηk2F + . 5 σr (X⋆ ) 10 σr (X⋆ )
This proves part (ii). 39
C.4
Comparison of two nearby matrices
We first compare aligned singular factors and graph cores at two nearby matrices. C.4.1
Variation of the aligned factors
The next lemma compares the singular factors and graph cores of two nearby matrices. We use the projector and Procrustes estimates in [31, Lemmas 4.1 and 4.5], aligning the left and right bases separately. Lemma C.1. Let X, Y ∈ Mr . For Z ∈ {X, Y }, represent Z = UZ SZ VZT TU T in independently aligned Procrustes gauges so that UX Y 0 and VX VY 0, and use the graph factors from (23) for each Z ∈ {X, Y }. Suppose that
n o σr (X⋆ ) . max kX − X⋆ k♯ , kY − X⋆ k♯ ≤ ρ ≤ 1000µr
Then T max kUX − UY kF , kVX − VY kF , UX UX − UY UYT
F
Moreover,
max {kLX − LY kF , kRX − RY kF } ≤
, VX VXT − VY VYT
ρ 1+3 σr (X⋆ )
and kGX − GY kF ≤ 5dX,Y ,
max
Z∈{X,Y }
G−1 Z op ≤
F
≤
3 dX,Y . 2 σr (X⋆ )
dX,Y ,
1 . σr (X⋆ ) − 2ρ
Finally, −1 G−1 X − GY F ≤
5dX,Y σr (X⋆ ) − 2ρ
2 .
T )X = 0, Y V S −1 = Proof. By Weyl’s inequality, min{σr (X), σr (Y )} > 0. Since (I − UX UX Y Y UY , and the transposed identities hold on the right, the Procrustes estimates in [31, Lemmas 4.1 and 4.5] give
3 dX,Y , 2 σr (X⋆ ) 3 dX,Y , ≤ F 2 σr (X⋆ )
T − UY UYT UX UX
3 dX,Y , 2 σr (X⋆ ) 3 dX,Y ≤ . F 2 σr (X⋆ )
kVX − VY kF ≤
kUX − UY kF ≤
VX VXT − VY VYT
This proves the factor and projector estimates. In the Procrustes gauges of the statement, set CV := VXT VY 0.
T UY 0, CU := UX
√ Since rank(X − Y ) ≤ 2r and max{kX − X⋆ k♯ , kY − X⋆ k♯ } ≤ ρ, we have dX,Y ≤ 2 2r ρ. Combining this with ρ ≤ σr (X⋆ )/(1000µr) and the projector estimates obtained from [31, Lemmas 4.1 and 4.5] shows that CU and CV are invertible. Their definitions give T I − CU2 = UYT (I − UX UX )UY ,
I − CV2 = VYT (I − VX VXT )VY .
40
Consequently, we have T (I − CU )SY = (I + CU )−1 UYT (I − UX UX )(Y − X)VY ,
SY (I − CV ) = UYT (Y − X)(I − VX VXT )VY (I + CV )−1 .
T XV T Since SX = UX X and SY = UY Y VY , substituting these identities gives T (X − Y )VX + (CU − I)SY CV + SY (CV − I), SX − SY = U X
kSX − SY kF ≤ 3dX,Y .
Using the definitions of the graph factors together with the factor and projector bounds, we have ρ dX,Y , σr (X⋆ ) ρ dX,Y . kRX − RY kF ≤ dX,Y + 3 σr (X⋆ ) kLX − LY kF ≤ dX,Y + 3
Moreover, we have T UX (X⋆ − X)VX − UYT (X⋆ − Y )VY
F
≤ dX,Y + 3
ρ dX,Y . σr (X⋆ )
Combining the preceding estimate with kSX − SY kF ≤ 3dX,Y yields kGX − GY kF ≤ 5dX,Y . For Z ∈ {X, Y }, the definition of GZ and Weyl’s inequality give G−1 Z op ≤
σmin (GZ ) ≥ σr (X⋆ ) − 2ρ,
1 . σr (X⋆ ) − 2ρ
The resolvent identity then gives −1 G−1 X − GY F ≤
This completes the proof.
5dX,Y σr (X⋆ ) − 2ρ
2 .
We now apply the factor estimates to the graph-retraction identity. The proof also verifies that both updated cores are invertible, so both retractions are well defined. C.4.2
Proof of Lemma 7.6
Proof. We align the singular factors as in Lemma C.1. For Z ∈ {X, Y }, write MZ := UZT ηZ VZ ,
BZ := (I − UZ UZT )ηZ VZ ,
and T CZ := (I − VZ VZT )ηZ UZ ,
Sη,Z := GZ + MZ .
sη := max{kηX kF , kηY kF },
∆ := dη + 3
Set
41
sη dX,Y . σr (X⋆ )
Using the aligned bases, we have MX − MY = (UX − UY )T ηX VX + UYT (ηX − ηY )VX + UYT ηY (VX − VY ), T BX − BY = (UY UYT − UX UX )ηX VX + (I − UY UYT )(ηX − ηY )VX
+ (I − UY UYT )ηY (VX − VY ),
T CX − CY = (VY VYT − VX VXT )ηX UX + (I − VY VYT )(ηX − ηY )T UX
+ (I − VY VYT )ηYT (UX − UY ).
Therefore, Lemma C.1 gives max {kMX − MY kF , kBX − BY kF , kCX − CY kF } ≤ ∆.
(76)
Moreover, the definitions of the tangent and graph coordinates imply o n max max {kMZ kF , kBZ kF , kCZ kF } ≤ sη , max max kLZ kop , kRZ kop ≤ ρ. Z∈{X,Y }
Z∈{X,Y }
The correction-size hypothesis gives sη ≤ 9µrρ2 /σr (X⋆ ).
(77) Furthermore, Weyl’s inequality
and (77) give
−1 Sη,Z
σmin (Sη,Z ) ≥ σr (X⋆ ) − 2ρ − sη > 0,
op
≤
1 , σr (X⋆ ) − 2ρ − sη
Z ∈ {X, Y }.
Using the resolvent identity, Lemma C.1, and (76), we also obtain −1 −1 − Sη,Y Sη,X
F
≤
5dX,Y + ∆ σr (X⋆ ) − 2ρ − sη
2 .
(78)
Since ρ < σr (X⋆ )/4, Lemma 5.3 applies at both base points. Hence (67) gives, for Z ∈ {X, Y }, we have RetrZ (tZ + ηZ ) − X⋆ − ηZ
−1 −1 −1 −1 T T T CZ − LZ G−1 CZ + BZ Sη,Z RTZ + LZ Sη,Z = BZ Sη,Z Z MZ Sη,Z RZ .
(79)
We compare the four terms in (79). In each product, we use the Frobenius norm for a correction or a difference and the operator norm for the remaining graph factors. For the first term, we have −1 −1 RTY RTX − BY Sη,Y BX Sη,X
−1 −1 −1 −1 (RX − RY )T . )RTX + BY Sη,Y − Sη,Y RTX + BY (Sη,X = (BX − BY )Sη,X −1 T . For CZ Interchanging the left and right factors gives the corresponding expansion for LZ Sη,Z the third term, we have T −1 −1 CYT CX − BY Sη,Y BX Sη,X
−1 −1 −1 −1 T T (CX − CY )T . )CX + BY Sη,Y − Sη,Y CX + BY (Sη,X = (BX − BY )Sη,X
Finally, the five-factor difference has the telescoping expansion −1 −1 −1 T T LX G−1 X MX Sη,X RX − LY GY MY Sη,Y RY
−1 −1 −1 −1 T T = (LX − LY )G−1 X MX Sη,X RX + LY (GX − GY )MX Sη,X RX
−1 −1 −1 −1 T T + LY G−1 Y (MX − MY )Sη,X RX + LY GY MY (Sη,X − Sη,Y )RX
T −1 + LY G−1 Y MY Sη,Y (RX − RY ) .
42
Applying Lemma C.1 and (76)–(78) termwise in these expansions, and then using ∆ = dη + 3sη dX,Y /σr (X⋆ ), sη ≤ 9µrρ2 /σr (X⋆ ), and ρ ≤ σr (X⋆ )/(1000µr), yields kRetrX (tX + ηX ) − RetrY (tY + ηY ) − (ηX − ηY )kF ≤
µrρ2 1 dη + 19 dX,Y . 4 σr (X⋆ )2
Adding kηX − ηY kF = dη proves (54).
D
Sharp-norm estimates for spectral initialization
We prove the deterministic reconstruction estimate and the finite approximation used in Section 6. The estimates concern the reconstructed matrix and do not require separation between adjacent singular values.
D.1
Proof of Lemma 6.2
Proof of Lemma 6.2. Suppose first that Z + 6= 0. With the dilation in (36), set " # b U b 1 U b −Σ), b Qb := √ , D := diag(Σ, P := Q⋆ Q⋆T . 2 Vb −Vb
The columns of Qb are orthonormal, and D −1 op ≤ 8/τ . The one-sided singular equation in (38a) and the residual bound in (38b) give r µr τ 0 X⋆ b b . (80) + W Q = QD + F , kF kop ≤ T 0 X⋆ 64 n The coefficient of the true singular spaces is 0 Σ⋆ b −1 . C := Q⋆T QD Σ⋆ 0 By (80), we have
b −1 + F D −1 , Q⋆ C = Qb − W QD
8 kCkop ≤ 1 + τ
τ τ + 512 64
r
µr n
≤
5 . 4
In particular, this estimate does not use the ratio σ1 (X⋆ )/τ . Since kW kop D −1 op ≤ 1/64, the same equation has the convergent expansion Qb =
For j ≥ 1, (37) implies
X j≥0
W j Q⋆ CD −j −
X
W j F D −j−1 .
(81)
j≥0
r r r µr τ j µr µr τ j j + . kW kop ≤ 3 (I − P )W Q⋆ 2,∞ ≤ 2 n 64 n n 64 j
Also, (I − P )W j F 2,∞ ≤ kW kjop kF kop . Applying these estimates to (81) gives r
r τ 15τ 3τ µr µr ≤ + , ≤ n 256(1 − 1/8) 64(1 − 1/64) 32 n 2,∞ r 9τ µr 3 µr b QbT (I − P ) b ≤ (I − P )QD , . ≤ (I − P )Q 4 n 128n ∞ 2,∞
b (I − P )QD
43
(82) (83)
b QbT is the symmetric dilation of Z + . By (38b), we have The matrix QD Z + − X⋆ op ≤
τ 113τ 7τ + = . 32 512 512
Decompose its error on the left by P and I − P . Since kQ⋆ k2,∞ ≤ +
+
T
max{ Z − X⋆ 2,∞ , (Z − X⋆ )
2,∞
}≤
113 3 + 512 32
p
(84)
µr/n, (82)–(84) yield
r r 161τ µr µr = . τ n 512 n
For the entrywise norm, decompose the error on both sides by the same projections. The truespace term is bounded by (µr/n) kZ + − X⋆ kop , the two mixed terms by 3τ µr/(32n) each, and the remaining term by (83). Therefore, 113 τ µr 6 9 245τ µr + Z − X⋆ ∞ ≤ + + = . 512 32 128 n 512n Combining these three bounds in (2) gives kZ + − X⋆ k♯ ≤ τ /4. If Z + = 0, then σ1 (X⋆ ) ≤ kY kop + kW kop ≤ 113τ /512. Incoherence gives kX⋆ k♯ = σ1 (X⋆ ), so the same conclusion follows.
D.2
Proof of Lemma 6.3
Proof of Lemma 6.3. We condition on the fixed matrix Y . Let d1 ≥ · · · ≥ dn > 0 be the eigenvalues of Y Y T + τ 2 I/4096, with corresponding orthonormal eigenvectors u1 , . . . , un . The hypothesis on σr+1 (Y ) gives dr+1 ≤
τ2 τ2 τ2 < . + 5122 4096 4000
(85)
Let U1 contain the eigenvectors whose eigenvalues are at least τ 2 /128, and let D1 be the diagonal matrix of those eigenvalues. By (85), their number is at most r. Let U2 and D2 contain the remaining eigenvectors and eigenvalues, so that YYT +
τ2 I = U1 D1 U1T + U2 D2 U2T , 4096
kD2 kop ≤
τ2 . 128
These decompositions are used only in the proof, not in the computation of T . ⊥ ], write the Step 1: the initial subspace. Put U[r] = [u1 , . . . , ur ]. In the eigenbasis [U[r] , U[r] Gaussian initial matrix as [GT1 , GT2 ]T , where G1 ∈ Rr×r . Orthogonal invariance preserves the standard Gaussian distribution. Since Pr{|g| ≤ t} ≤ t for g ∼ N (0, 1), −20 Pr min dist (G1 ):,j , span{(G1 ):,k : k 6= j} < n ≤ n−19 , 1≤j≤r
and therefore
≤ G−1 1 op Moreover,
E kG2 k2F ≤ n2 ,
√ 20 rn .
Pr{kG2 kF > n11 } ≤ n−20 .
Hence, with probability at least 1 − 2n−19 ≥ 1 − n−12 /48, ≤ n32 . G2 G−1 1 op
44
(86)
Step 2: weighted subspace estimates. Fix a realization satisfying (86). Thin QR does not change the column space, so (11) gives t ! 2 τ I G . range(Qt ) = range YYT + 4096 ⊥ , this space has graph matrix Relative to U[r] ⊕ U[r] −t −t diag(dtr+1 , . . . , dtn )G2 G−1 1 diag(d1 , . . . , dr ).
Since range(U1 ) ⊆ range(U[r] ), (85) and dj ≥ τ 2 /128 on range(U1 ) imply, for t ≥ 1 and s ∈ {0, 1/2, 1}, we have 2 s τ T s (I − Qt Qt )U1 D1 op ≤ n32 30−t . (87) 128 For t = ⌈12 log n⌉ and n ≥ 2, we have 32
−t
n 30
≤n
−7
1 ≤ 64
r
µr . n
Let Q denote the final factor, and put P = QQT and E = (I − P )U1 . Then r s 1 µr τ 2 s , s ∈ {0, 1/2, 1}. (88) kED1 kop ≤ 64 n 128 √ If U1 has no columns, then kY kop < τ /(8 2) < τ /8. No singular value is retained and the conclusion is immediate. We henceforth suppose that U1 is nonempty. Step 3: the retained singular factors. An orthonormal basis for range(P U1 ) is e 1 = (U1 − E)(I − E T E)−1/2 , Q
(I − E T E)−1/2
op
≤
4 . 3
e 2 . Then Q e T U1 = 0. Since U T E = E T E, direct substitution Complete it inside range(Q) by Q 2 1 gives e 1 = ED1 − ED1 E T E − (I − P )U2 D2 U T E (I − E T E)−1/2 , (I − P )Y Y T Q 2 eT Y Y T Q e 1 = −Q e T U2 D2 U T E(I − E T E)−1/2 . Q 2
2
2
Using (88) and µr/n ≤ 1, these identities imply r µr τ2 T e , ≤ (I − P )Y Y Q1 2048 n op r τ2 µr T T e e ≤ Q2 Y Y Q1 , 4096 n op
(89a) 2
eT Y Y T Q e 2 τ I. Q 2 128
(89b)
b = Vb Σ b and σmin (Σ) b ≥ τ /8. Its retained left factors The compressed SVD in (12) gives Y T U T b e e are Ritz vectors of Y Y in range(Q). Write U = Q1 Z1 + Q2 Z2 . The lower block of the Ritz equation is b 2 − (Q eT Y Y T Q e 2 )Z2 = Q eT Y Y T Q e 1 Z1 . Z2 Σ 2 2
The spectra on the two sides are separated by at least τ 2 /128. The integral solution of this Sylvester equation, together with (89b) and kZ1 kop ≤ 1, therefore gives r 1 µr 128 e T T e ≤ . kZ2 kop ≤ 2 Q2 Y Y Q1 Z1 τ 32 n op 45
e 2 is orthogonal to U1 , also Y Y T Q e2 Since Q
op
≤ τ 2 /128. Combining this with (89a) yields
τ2 (I − P )Y Y U ≤ 1024 op
b Σ. b Consequently, Moreover, P Y Vb = U bΣ b Y Vb − U
op
T b
T b b −1
= (I − P )Y Y U Σ
r
µr . n
τ ≤ 128 op
r
τ µr ≤ n 64
r
µr . n
If no singular value is retained, this residual estimate is unnecessary. Step 4: the unreconstructed part. The s = 1/2 estimate in (88) controls the part of Y in p range(U1 ), since σj (Y ) ≤ dj . On its orthogonal complement the singular values are smaller √ than τ /(8 2). Thus r τ 3τ τ µr + √ < . k(I − P )Y kop ≤ 512 n 32 8 2 The singular values discarded from P Y are smaller than τ /8. Therefore, Y − Z + op ≤ k(I − P )Y kop + P Y − QH≥τ /8 (QT Y ) op ≤
3τ τ 7τ + = . 32 8 32
This includes the case Z + = 0. Thus all approximation conclusions hold on (86). Since Y Y T + τ 2 I/4096 ≻ 0 and G has full column rank almost surely, every thin QR step is well defined. A fixed orthonormal completion defines the procedure on the remaining null event without changing its operation count.
E
Leave-one-out estimates for the initial RGN iterates
We prove the sampling and correction estimates used to control the RGN iterates for 0 ≤ k ≤ K̄. The concentration and geometric tools are given in Sections B and C. The passage to the quadratic Frobenius recurrence is proved in Subsection 7.4.
E.1
Sampling estimates for the completed operators
We use the completed sampling operators Rα and the index set A defined in Subsection 7.1. For α ∈ A, let Aα be the entrywise nonnegative square root of Rα , so (Aα )∗ Aα = Rα . We first establish the simultaneous sampling isometries for the actual and completed operators. E.1.1
Proof of Lemma 7.1
T Proof. Let δij := 1{(i,j)∈Ω} b and aij = PTX⋆ (ei ej ), where (a ⊗ a)(Z) := ha, ZiF a.
kaij k2F
By
Assumption 2.1, ≤ 2µr/n. Since the matrices ei eTj form an orthonormal basis and PTX⋆ is an orthogonal projector, we have X i,j
aij ⊗ aij = ITX⋆ ,
(aij ⊗ aij )2 = kaij k2F (aij ⊗ aij ), 2µr 1X kaij k2F (aij ⊗ aij ) IT . q nq X⋆ i,j
46
For the actual operator, PTX⋆ R0 PTX⋆ − PTX⋆ is the centered sum of (δij /q − 1)(aij ⊗ aij ) over all coordinates. For Rr,i , the centered sum is restricted to coordinates (a, b) with a 6= i; for Rc,j , it is restricted to coordinates with b 6= j. For the actual centered sum, each summand and the variance operator satisfy 2µr , nq X 1−q X 2µr kaij k2F (aij ⊗ aij ) IT . E (δij /q − 1)2 (aij ⊗ aij )2 = q nq X⋆ k(δij /q − 1)(aij ⊗ aij )kF→F ≤
i,j
i,j
Restricting the sums to a completed row or column preserves these two inequalities. Since q ≥ c1 µr log n/n and c1 = 215 , the Bernstein exponent at threshold 1/16 is at least 24 log n. Since dim(TX⋆ ) ≤ 2nr and |A| = 2n + 1, Lemma B.1 followed by a union bound gives 1 n−10 α Pr max PTX⋆ R PTX⋆ − PTX⋆ F→F > . ≤ 4nr(2n + 1)n−24 ≤ α∈A 16 24 We next record the weighted product estimate used to transfer these sampling bounds to nearby tangent spaces. E.1.2
A weighted product estimate
b q). Then, for Lemma E.1. Suppose that the conclusion of Lemma 7.2 holds with (Λ, pΛ ) = (Ω, every α ∈ A and all matrices F , H with n rows and the same number of columns, n o √ Aα (F H T ) F ≤ 3n min kF kF kHk2,∞ , kF k2,∞ kHkF .
Proof. For the actual operator, weighted row and column sums are the degrees divided by q, hence at most 2n. Completing one row or column adds at most one to every opposite weighted degree and replaces the completed weighted degree by n, so all weighted sums are at most 3n. α are the weights of Rα , then If fiT and hTj are the rows of F and H, and wij 2
Aα (F H T ) F = ≤
X
α wij (fiT hj )2
i,j
X
kfi k22
X
α wij khj k22
j 2 ≤ 3n kF kF kHk22,∞ . i
Interchanging rows and columns gives 2
Aα (F H T ) F ≤ 3n kF k22,∞ kHk2F . Taking square roots in the two preceding inequalities and minimizing the resulting upper bounds proves the lemma.
E.2
Proof of Lemma 7.3
We first prove the conditional row estimate; the corresponding column estimate follows by transposition in the subsequent leave-one-out argument. The Proof of Lemma 7.3 is after the following lemma.
47
Lemma E.2. Let δ1 , . . . , δn be independent Bernoulli(pΛ ) variables, and let z ∈ Rn and V ∈ Rn×r be deterministic. Set δj n − 1 zj . w = (wj )j=1 , wj = pΛ p If V T V = I and kV k2,∞ ≤ 2 µr/n, then, with probability at least 1 − (n + r + 2)n−24 , s p r r 64 µr/n log n µr 3 log n µr T 2 kwk2 + V w 2 ≤ 8 kzk2 + kzk∞ + kzk∞ . 4 n pΛ n pΛ
P T ≤ kzk2∞ /pΛ , Proof. For Zj := (p−1 j E(Zj Zj ) Λ δj − 1)zj ej , we have kZj k2 ≤ kzk∞ /pΛ , op P and j E kZj k22 ≤ kzk22 /pΛ . Thus the variance parameter in Lemma B.1 satisfies v ≤ kzk22 /pΛ , and applying the lemma with t = 24 log n gives s kzk22 kzk∞ 16 log n 3 log n , v≤ , kwk2 ≤ 4 kzk2 + kzk∞ . L≤ pΛ pΛ pΛ pΛ p ej := (p−1 δj − 1)zj V T ej , the incoherence bound kV k ≤ 2 For Z µr/n gives 2,∞ Λ p 2 µr/n kzk∞ /pΛ . Moreover, X j
ej Z ejT ) E(Z
≤ kzk2∞ /pΛ
and
X j
op
ej E Z
2 2
ej Z
2
≤
≤ 4µr kzk22 /(npΛ ).
Thus we may take v ≤ kzk2∞ /pΛ + 4µr kzk22 /(npΛ ) in Lemma B.1. Applying the lemma with t = 24 log n gives s p r 32 µr/n log n 3 log n µr T V w 2≤4 kzk2 + kzk∞ + kzk∞ . 2 pΛ n pΛ The Bernstein bounds for kwk2 and V T w 2 fail with probability at most (n + 1)n−24 and (r + 1)n−24 , respectively. Therefore, by a union bound, both estimates hold with probability at least 1 − (n + r + 2)n−24 . On the intersection of the two Bernstein events, adding the two estimates and collecting the kzk2 and kzk∞ terms gives s s p r µr/n log n 64 µr 3µr log n 3 log n kwk2 + V T w 2 ≤ 16 kzk2 + 4 kzk∞ + kzk∞ 2 n npΛ pΛ pΛ s p r 64 µr/n log n 3 log n µr ≤8 kzk2 + kzk∞ + kzk∞ , 4 pΛ n pΛ which is the claimed inequality. corr defined in We use the stopped auxiliary sequences and the quantities ρk , dloo k , and hk Subsection 7.1. With the fixed compact-SVD convention in Lemma 7.3, write
Xkα = Ukα Skα (Vkα )T , and set Nkα := NXkα as in (24). For each comparison, the left and right singular bases are aligned as in Section C.Now we give the proof of Lemma 7.3.
48
Proof of Lemma 7.3. Fix a row i and condition on X0 and the Bernoulli variables outside row i. The stopped row-i trajectory, its residual row, and its right singular factors are then fixed, whereas the row-i Bernoulli variables remain independent. Before stopping, the initial error bound and the accepted updates give Xkr,i − X⋆ ≤ ρk . ♯
Therefore, the proof of Lemma 5.1 gives, for every such k < K̄, σr (X⋆ ) . ≤ ρk ≤ ρ0 ≤ 8 ♯ r r 10 µr µr r,i ≤ <2 . Vk 7 n n 2,∞
Xkr,i − X⋆
p After stopping, the sequence is equal to X⋆ , and Assumption 2.1 gives kV⋆ k2,∞ ≤ µr/n < p 2 µr/n. Thus Lemma E.2 applies at every k < K̄. Applying Lemma E.2 to the transposed columnj process, conditional on the variables outside column j, gives the column inequalities in Lemma 7.3. Taking expectations over the conditioned variables and a union bound over at most 2nK̄ pairs gives failure probability at most n−10 /48 for K̄ ≤ n. The row and column bounds in Lemma 7.3 are obtained before alignment and are invariant under the subsequent left and right orthogonal Procrustes transformations. Consequently, conditional on every admissible realization of X0 , all row and column inequalities in the lemma hold simultaneously with failure probability at most n−10 /48.
E.3
Local solvability and correction size
We next use the simultaneous sampling events to prove local invertibility and bound the tangent correction. Proof of Lemma 7.4. Since µr ≥ 1 and e♯ ≤ σr (X⋆ )/(1000µr), Lemmas 5.1 and E.1 and (4) give r 64 µr e♯ T T T T U U − U⋆ U⋆ 2,∞ + V V − V⋆ V⋆ 2,∞ ≤ , 7 n σr (X⋆ ) √ e♯ 64 3 √ (90) , µr Aα (PTX − PTX⋆ ) F→F ≤ 7 σr (X⋆ ) 16 e♯ . PTX − PTX⋆ F→F ≤ 7 σr (X⋆ ) Fix ζ ∈ TX Mr . By Lemma 7.1, the triangle inequality, and (90), we have kAα ζkF ≥ Aα PTX⋆ ζ F − Aα (PTX − PTX⋆ )ζ F r √ e♯ 64 3 √ 15 ≥ PTX⋆ ζ F − kζkF µr 16 7 σr (X⋆ ) r 9 ≥ kζkF . 10 Using the upper inequality in Lemma 7.1, the triangle inequality, and (90), we also have kAα ζkF ≤ Aα PTX⋆ ζ F + Aα (PTX − PTX⋆ )ζ F r √ e♯ 17 64 3 √ µr PTX⋆ ζ F + kζkF ≤ 16 7 σr (X⋆ ) r 11 kζkF . ≤ 10 49
Here we used PTX⋆ ζ F ≥ kζkF − (PTX − PTX⋆ )ζ F in the lower bound and PTX⋆ ζ F ≤ kζkF in the upper bound. Squaring the inequalities p p kAα ζkF ≥ 9/10 kζkF and kAα ζkF ≤ 11/10 kζkF , and using (Aα )∗ Aα = Rα , proves (49). α , and η α give We next prove the correction estimate. The definitions of tX , ξX X α PTX Rα PTX ηX = PTX Rα NX .
Write X = U ΣV T as a compact singular value decomposition. By (49), Lemmas E.1 and 5.2, and Weyl’s inequality, we obtain 9 σmin (GX ) ≥ σr (X⋆ ) − 2e♯ ≥ σr (X⋆ ), 10 r 10 11 α kηX kF ≤ kAα NX kF 9 10 r 10 11 √ 3n LX G−1 ≤ X F kRX k2,∞ 9 10 e2♯ ≤ 9µr . σr (X⋆ ) Therefore, (50) follows.
E.4
Comparison of the tangent corrections
We first bound the variation in the sampled normal component between two nearby matrices. Lemma E.3. Let α ∈ A and X, Y ∈ Mr , and use dX,Y from Lemma C.1. Suppose that Assumption 2.1 holds and o n σr (X⋆ ) max kX − X⋆ k♯ , kY − X⋆ k♯ ≤ ρ ≤ . 1000µr p If dX,Y ≤ 2 µr/n ρ, then, on the joint events in Lemmas 7.1 and 7.2, it holds r ρ2 µr α kPTX R (NX − NY )kF ≤ 32 µr . n σr (X⋆ ) Proof. By the graph-factor representation in Lemma 5.2, we have −1 −1 T T NX − NY = (LX − LY )G−1 X RX + LY (GX − GY )RX T + LY G−1 Y (RX − RY ) .
(91)
Applying Lemmas C.1, E.1 and 5.2 termwise in (91), and using ρ ≤ σr (X⋆ )/(1000µr), gives √ ρ dX,Y kAα (NX − NY )kF ≤ 15 µr . σr (X⋆ ) p Since dX,Y ≤ 2 µr/n ρ, it follows that α
kA (NX − NY )kF ≤ 30µr
Finally, (49) and (Aα )∗ Aα = Rα give α
kPTX R (NX − NY )kF ≤
r
r
µr ρ2 . n σr (X⋆ )
r ρ2 11 µr α kA (NX − NY )kF ≤ 32 µr . 10 n σr (X⋆ )
This proves the claim. 50
We now estimate the tangent transport, which bound the error incurred by projecting a tangent correction onto a nearby tangent space. Lemma E.4. Let α ∈ A, h ≥ 0, X, Y ∈ Mrp , and Z ∈ TY Mr . Suppose that X and Y satisfy the hypotheses of Lemma E.3, that dX,Y ≤ 2 µr/n ρ, and that r n o ρ2 µr ρ2 T , max kZk2,∞ , Z 2,∞ ≤ 2h + 38µr . kZkF ≤ 9µr σr (X⋆ ) n σr (X⋆ ) Then, on the joint events in Lemmas 7.1 and 7.2, 1 kPTX R (PTY − PTX )ZkF ≤ µr 2 α
r
1 µr ρ2 + h. n σr (X⋆ ) 40
Proof. Write Y = UY SY VYT and use the tangent decomposition Z = UY M VYT + BVYT + UY C T ,
UYT B = 0,
VYT C = 0.
Since the three terms are mutually orthogonal, we have max {kM kF , kBkF , kCkF } ≤ kZkF . Moreover, B = ZVY − UY M and C = Z T UY − VY M T . Therefore, Lemma 5.1 and the assumed bounds on Z give r o n µr ρ2 . max kBk2,∞ , kCk2,∞ ≤ 2h + 56µr n σr (X⋆ ) Define T ∆U := (I − UX UX )UY ,
∆V := (I − VX VXT )VY .
Since we have ∆V = (I − VX VXT )(Y − X)T UY SY−T ,
T ∆U = (I − UX UX )(Y − X)VY SY−1 ,
Weyl’s inequality gives max {k∆U kF , k∆V kF } ≤
dX,Y . σr (X⋆ ) − ρ
(92)
T )(Y − X)V S −1 gives For every i, the identity ∆U = (I − UX UX Y Y
eTi ∆U 2 ≤
eTi (Y − X) 2 + eTi UX 2 kY − Xkop σr (Y )
r ρ µr ≤9 , n σr (X⋆ )
where the last inequality follows from the definition of the sharp norm, Lemma 5.1, and σr (Y ) ≥ σr (X⋆ ) − ρ. Applying the same calculation to ∆V = (I − VX VXT )(Y − X)T UY SY−T gives r n o ρ µr . (93) max k∆U k2,∞ , k∆V k2,∞ ≤ 9 n σr (X⋆ ) Furthermore, since UYT B = 0 and VYT C = 0, we have T )B 2,∞ ≤ kBk2,∞ + kUX k2,∞ (UX − UY )T B F , (I − UX UX
(I − VX VXT )C 2,∞ ≤ kCk2,∞ + kVX k2,∞ (VX − VY )T C F .
51
p Combining these inequalities with Lemma C.1, dX,Y ≤ 2 µr/n ρ, and the Frobenius bound on Z, we obtain r n o µr ρ2 T T . (94) max (I − UX UX )B 2,∞ , (I − VX VX )C 2,∞ ≤ 2h + 57µr n σr (X⋆ ) Since Z ∈ TY Mr , substituting its tangent decomposition into (4) gives T (PTY − PTX )Z =∆U M ∆TV + (I − UX UX )B∆TV + ∆U C T (I − VX VXT ).
(95)
Applying Lemma E.1 to the first term and using (92)–(93), we have p Aα (∆U M ∆TV ) F ≤ 9 3µr
ρ k∆U kF kM kF . σr (X⋆ )
Applying Lemma E.1 to the second and third terms in (95), and then using (92) and (94), we obtain T Aα ((I − UX UX )B∆TV ) F + Aα (∆U C T (I − VX VXT )) F r √ dX,Y µr ρ2 ≤ 2 3n 2h + 57µr . σr (X⋆ ) − ρ n σr (X⋆ ) p Substituting dX,Y ≤ 2 µr/n ρ, ρ ≤ σr (X⋆ )/(1000µr), kM kF ≤ kZkF , and µr ≤ n into the bounds for the three terms in (95) gives r 2 1 µr ρ2 α kA (PTY − PTX )ZkF ≤ µr + h. (96) 5 n σr (X⋆ ) 45
Finally, (49), (Aα )∗ Aα = Rα , and (96) give r r 1 1 11 µr ρ2 α α kA (PTY − PTX )ZkF ≤ µr + h. kPTX R (PTY − PTX )ZkF ≤ 10 2 n σr (X⋆ ) 40 This proves the result.
E.5
Proof of Lemma 7.5
We now combine the normal-motion and tangent-transport estimates to compare the actual and leave-one-out tangent corrections. Proof of Lemma 7.5. First estimate: coordinatewise bounds. Fix α and a row i. By the definicorr and by Lemma C.1 in the aligned gauges, we have tions of dloo k and hk ηkα − ηkr,i
F
≤ 2hcorr k ,
Xkα − Xkr,i
F
≤ 2dloo k ,
Vkα − Vkr,i
op
≤3
dloo k . σr (X⋆ )
Set X = Xkr,i = U SV T , N = Nkr,i , and η = ηkr,i . Since Rr,i is the identity on row i and N V = 0, the correction equation and, for every a ∈ Rr , its test against ei aT V T give PTX Rr,i (N − η) = 0, ei aT V T = U (U T ei )aT V T + (I − U U T )ei aT V T ∈ TX Mr , eTi (N − η)V = 0,
eTi ηkr,i Vkr,i = 0.
52
Since eTi ηkr,i Vkr,i = 0, the ith row of ηkr,i contains only its U C T component. Hence Lemma 5.1, (50), µr/n ≤ 1 from Assumption 2.1, and ρk /σr (X⋆ ) ≤ 1/(1000µr) give r r ρ2k µr r,i µr T r,i µr , ≤ 18 ηk ≤2 ei ηk n n σr (X⋆ ) F 2 r ρ2k dloo µr k eTi ηkα Vkα 2 ≤ 2hcorr + 54 µr k n σr (X⋆ ) σr (X⋆ ) r ρ2k 1 µr corr µr . (97) ≤ 2hk + 8 n σr (X⋆ ) Replacing (r, i, U , V ) by (c, j, V , U ) and transposing gives the column form of (97). Moreover, comparing each row with its row-completed correction, and each column with its columncompleted correction, gives r o n ρ2k µr α T corr α µr , α ∈ A. (98) max kηk k2,∞ , (ηk ) 2,∞ ≤ 2hk + 18 n σr (X⋆ ) For a tangent matrix η at X = U ΣV T , we have η = (ηV )V T + U (η T U )T − U (U T ηV )V T . Equation (97) and its transpose, together with the tangent decomposition, Lemma 5.1, and kηkF ≤ 9µrρ2k /σr (X⋆ ), give ρ2k , kηkop ≤ 9µr σr (X⋆ )
r ρ2k n corr kηk♯ ≤ 2 hk + 19µr . µr σr (X⋆ )
This proves (51). Second estimate: comparison of the corrections. Fix a row i, and set X = Xk0 ,
Y = Xkr,i ,
η = ηk0 ,
η = ηkr,i ,
Z = NY − η.
The two correction equations are PTX R0 PTX η = PTX R0 NX ,
(99a)
PTY Rr,i PTY η = PTY Rr,i NY .
(99b)
Since η ∈ TX and η ∈ TY , (99b) gives PTY Rr,i Z = 0. Subtracting (99b) from (99a) and adding and subtracting NY and Rr,i give (PTX R0 PTX )(η − PTX η) = PTX R0 (NX − NY ) + PTX (R0 − Rr,i )Z {z } | {z } | T1
r,i
T2 0
+ (PTX − PTY )R Z + PTX R (PTY − PTX )η . | {z } | {z } T3
(100)
T4
Replacing (r, i) by (c, j) and transposing gives the column form of (100). By Lemma 7.4, the inverse of the left-hand operator has norm at most 10/9. We bound T1 , . . . , T4 separately, as in [21, Lemma 11]. r ρ2k µr For T1 , Lemma E.3 gives kT1 kF ≤ 32 µr . n σr (X⋆ ) For T2 , write the completed row of Z as ziT , set δij := 1{(i,j)∈Ω} b , and write wj = (δij /q −
1)(zi )j . Conditional on X0 and the Bernoulli variables outside row i, zi and Vkr,i are fixed, 53
whereas {δij }nj=1 remain independent Bernoulli(q) variables. Hence the conclusion of Lemma 7.3 gives (47) with z = zi and V = Vkr,i . Write Xkr,i = U ΣV T , N = Nkr,i , and η = ηkr,i . Testing the correction equation against T ei a V T ∈ TX r,i Mr and using N V = 0 give k
eTi (N − η)V = 0,
eTi ηV = 0.
Set Coff := (I − V V T )η T U . Using Lemma 5.2 for Nkr,i , the identity eTi ηV = 0 and Lemma 7.4 for eTi η, and the transposed form of (97) for Coff , the tangent decomposition of η gives r 16 µr ρ2k r,i T , ≤ ei Nk 3 n σr (X⋆ ) 2 16 4µr ρ2k ≤ Nkr,i , ∞ 3 n σr (X⋆ ) r r ρ2k µr µr T kηkF ≤ 18 µr , ei η 2 ≤ 2 n n σr (X⋆ ) eTj Coff 2 ≤ eTj η T U 2 + eTj V 2 U T ηV op r r ρ2k ρ2k 1 µr µr corr ≤ 2hk + µr + 18 µr , 8 n σr (X⋆ ) n σr (X⋆ ) r µr ρ2k , kzi k2 ≤ 24µr n σr (X⋆ ) r 60(µr)2 ρ2k µr corr kzi k∞ ≤ +4 h . n σr (X⋆ ) n k
(101a) (101b)
Substituting (101a) and (101b) into (47), and using q ≥ c4 µr log n/n with c4 = 220 , gives r r ρ2k 17 µr µr r,i T kwk2 + (Vk ) w ≤ 13 µr + hcorr . 2 2 n n σr (X⋆ ) 50 k Equation (4), Lemma C.1, and ρk ≤ σr (X⋆ )/(1000µr) give 2
2
2
2
PTY (ei w T ) F = UYT ei 2 kwk22 + (I − UY UYT )ei 2 VYT w 2 , r dloo 3ρk µr T k kwk2 ≤ kwk2 + VY w 2 , 2 3 σr (X⋆ ) σr (X⋆ ) n r 7 µr ρ2k + hcorr . kT2 kF ≤ 14µr n σr (X⋆ ) 20 k For T3 , we use the row-coordinate estimate (97) and its transpose. Together with the tangent decomposition and the Frobenius bound on η, they give r (µr)2 ρ2k µr corr hk + 37 . kηk∞ ≤ 8 n n σr (X⋆ ) p Consequently, (99b), Lemmas C.1, 5.2 and 7.2, dloo ≤ 2 µr/n ρk , and ρk ≤ σr (X⋆ )/(1000µr) k
54
give r µr corr (µr)2 ρ2k hk + 64 , kZk∞ ≤ 8 n n σr (X⋆ )
UYT Rr,i Z = 0,
Rr,i ZVY = 0,
Rr,i Z op ≤ 3n kZk∞ ,
dloo k Rr,i Z op , σr (X⋆ ) r ρ2k 3 6 µr . µr + hcorr kT3 kF ≤ 5 n σr (X⋆ ) 20 k
(PTX − PTY )Rr,i Z F ≤ 3
For T4 , (50) and (98), followed by Lemma E.4, give kηkF ≤ 9µr max{kηk2,∞ , η
T
ρ2k , σr (X⋆ )
r
µr ρ2k , n σr (X⋆ ) r ρ2k 1 1 µr µr + hcorr . kT4 kF ≤ 2 n σr (X⋆ ) 40 k + 38µr } ≤ 2hcorr k 2,∞
To compare the two corrections themselves, we also account for the component of η normal to TX . By Lemmas C.1 and 7.4 and (100), we have 4
kη − ηkF ≤
3dloo 10 X k kTj kF + kηkF . 9 σr (X⋆ ) j=1
Substituting the four estimates for kT1 kF , . . . , kT4 kF derived in this proof, using p loo dk ≤ 2 µr/n ρk , and applying their transposed counterparts to each completed column yield r ρ2k 7 ρk µr corr . µr + hcorr hk ≤ 53 + 54 σr (X⋆ ) n σr (X⋆ ) 12 k Since ρk /σr (X⋆ ) ≤ 1/(1000µr) and µr ≥ 1, rearranging gives r ρ2k µr corr hk ≤ 128 µr . n σr (X⋆ ) This proves (52). Together with (51), this completes the proof.
F
Matrix-free implementation and computational complexity
We describe the matrix-free realization of the RGD and RGN steps and record the resulting complexity. Let X = U SV T ∈ Mr be a compact singular value decomposition. Every ξ ∈ TX Mr can be written as ξ = U M V T + BV T + U C T ,
U T B = 0,
V T C = 0;
see [28, Sec. 2.1]. Define ΦX (M , B, C) := U M V T + BV T + U C T . Then ΦX is an isometry onto TX Mr with Φ∗X (Z) = U T ZV , (I − U U T )ZV , (I − V V T )Z T U . 55
For the RGD direction, orthogonality gives kξk2F = kM k2F + kBk2F + kCk2F . Thus the numerator in (6) costs O(nr) operations. Its denominator is obtained by first forming U M , then evaluating ΦX (M , B, C) on Ω and summing the squared entries, which costs O(|Ω|r + nr 2 ) operations. The graph retraction costs O(nr 2 + r 3 ) operations. Since r ≤ n, the exact line search and the graph retraction therefore preserve the O(|Ω|r + nr 2 ) cost per RGD iteration. For RGN, writing z = (M , B, C), the tangent normal equation (8) is equivalent to HX z = −gX ,
HX := Φ∗X q −1 PΩb ΦX ,
gX := Φ∗X q −1 PΩb (X − X⋆ ).
(102)
Thus no r(2n − r) × r(2n − r) normal matrix is formed. Both the formation of gX and one b + nr 2 ) operations, and the graph retraction costs O(nr 2 + r 3 ) application of HX cost O(|Ω|r operations; see [1]. This gives the standard orders for fixed-rank Riemannian matrix completion; see, for example, [31]. By Lemma 7.4, HX is positive definite along the RGN iterates covered by Theorem 4.2. Therefore, exact CG applied to (102) is well defined and terminates in at most dim(TX Mr ) = r(2n − r) iterations [11]. In particular, 1 ≤ Jk ≤ r(2n − r) whenever an RGN correction is computed, and the kth RGN iteration costs b + nr 2 ) . O Jk (|Ω|r
F.1
Proof of Proposition 2.2
Proof. Set Y = Z + q −1 PΛ (X⋆ − Z) = U ΣV T + S,
where rank(Z) ≤ r and S has at most m nonzero entries. For Q ∈ Rn×r , Y T Q = V Σ(U T Q) + S T Q,
Y Q = U Σ(V T Q) + SQ.
Thus the two matrix products in each step of (11), together with the thin QR factorization, require O(mr + nr 2 ) operations. Since there are ⌈12 log n⌉ such steps, the subspace iteration 2 costs O (mr + nr ) log n . For the final spectral reconstruction, compute e Y T Q = QR,
Then and hence
RT = UR ΣR VRT .
e R )T , H≥τ /8 (QT Y ) = UR H≥τ /8 (ΣR )(QV e R )T . Tτ (Y ) = QUR H≥τ /8 (ΣR )(QV
The thin QR factorization and the r × r singular value decomposition require O(mr + nr 2 ) additional operations. Evaluating the entries of Z on Λ and summing the squared residuals costs O(mr) operations. This proves the result.
F.2
Computational complexity of RGD and RGN
Theorem F.1. Let 0 < ε ≤ 1/4. Under the conditions of Theorem 3.1, RGD attains relative Frobenius accuracy ε with complexity 1 2 O µnr log n log(nκ) log n + log . ε 56
Under the conditions of Theorem 3.2, RGN attains relative Frobenius accuracy ε with arithmetic complexity 1 , O µnr 2 log n log n log(2µrκ) + J log log ε where Jk ≤ J ≤ r(2n − r) for all RGN iterations used to attain the prescribed accuracy. Proof. Initialization. Let mℓ = |Ω(ℓ) | and let B = K + 1 for RGD and B = K + 2 for RGN. By Proposition 2.2, at most K reconstruction steps require # ! " K X 2 O r mℓ + Knr log n ℓ=1
P operations. The residual tests cost at most O(r K and have lower order. Since ℓ=1 mℓ ) operations PB−1 the observation components are independent Bernoulli(q) sets, ℓ=0 mℓ = O(Bqn2 ) with high probability. Moreover, Bq = O(p) and qn & µr log n, so Knr 2 = O(Bqn2 r). Therefore the initialization cost is O(pn2 r log n). (103) RGD. By the matrix-free implementation above, one RGD iteration costs O(|Ω|r + nr 2 ). At the sampling rate of Theorem 3.1, |Ω| = O(µnr log n log(nκ)) with high probability. Combining the resulting per-iteration cost, the linear convergence in Theorem 3.1, and (103) gives 1 O µnr 2 log n log(nκ) log n + log . ε RGN. By the matrix-free implementation above, the kth RGN iteration costs 2 b b = O(µnr log n) with high O Jk (|Ω|r + nr ) . At the sampling rate of Theorem 3.2, |Ω| probability. Combining this estimate, the convergence rate in Theorem 3.2, and (103) gives 1 2 . O µnr log n log n log(2µrκ) + J log log ε
References [1] P.-A. Absil and Ivan V. Oseledets. Low-rank retractions: A survey and new results. Computational Optimization and Applications, 62(1):5–29, 2015. [2] Jeffrey D. Blanchard, Jared Tanner, and Ke Wei. CGIHT: Conjugate gradient iterative hard thresholding for compressed sensing and matrix completion. Information and Inference: A Journal of the IMA, 4(4):289–327, 2015. [3] Hanqin Cai, Jian-Feng Cai, and Juntao You. Structured gradient descent for fast robust low-rank Hankel matrix completion. SIAM Journal on Scientific Computing, 45(3):A1172– A1198, 2023. [4] HanQin Cai, Longxiu Huang, Xiliang Lu, and Juntao You. Accelerating ill-conditioned Hankel matrix recovery via structured Newton-like descent. Inverse Problems, 41(7):075015, 2025. [5] Jian-Feng Cai, Tong Wu, and Ruizhe Xia. Fast non-convex matrix sensing with optimal sample complexity. In Proceedings of the Forty-first Conference on Uncertainty in Artificial Intelligence, pages 497–520, 2025. 57
[6] Emmanuel J. Candès and Benjamin Recht. Exact matrix completion via convex optimization. Foundations of Computational Mathematics, 9(6):717–772, 2009. [7] Emmanuel J. Candès and Terence Tao. The power of convex relaxation: Near-optimal matrix completion. IEEE Transactions on Information Theory, 56(5):2053–2080, 2010. [8] Yudong Chen. Incoherence-optimal matrix completion. IEEE Transactions on Information Theory, 61(5):2909–2923, 2015. [9] Yudong Chen, Srinadh Bhojanapalli, Sujay Sanghavi, and Rachel Ward. Completing any low-rank matrix, provably. Journal of Machine Learning Research, 16(94):2999–3034, 2015. [10] Lijun Ding and Yudong Chen. Leave-one-out approach for matrix completion: Primal and dual analysis. IEEE Transactions on Information Theory, 66(11):7274–7301, 2020. [11] Anne Greenbaum. Iterative Methods for Solving Linear Systems, volume 17 of Frontiers in Applied Mathematics. Society for Industrial and Applied Mathematics, Philadelphia, PA, 1997. [12] Nathan Halko, Per-Gunnar Martinsson, and Joel A. Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review, 53(2):217–288, 2011. [13] Moritz Hardt. Understanding alternating minimization for matrix completion. In 2014 IEEE 55th Annual Symposium on Foundations of Computer Science, pages 651–660, 2014. [14] Moritz Hardt and Mary Wootters. Fast matrix completion without the condition number. In Proceedings of the 27th Conference on Learning Theory, pages 638–678, 2014. [15] Prateek Jain, Raghu Meka, and Inderjit S. Dhillon. Guaranteed rank minimization via singular value projection. In Advances in Neural Information Processing Systems, volume 23, pages 937–945, 2010. [16] Prateek Jain and Praneeth Netrapalli. Fast exact matrix completion with finite samples. In Proceedings of the 28th Conference on Learning Theory, pages 1007–1034, 2015. [17] Prateek Jain, Praneeth Netrapalli, and Sujay Sanghavi. Low-rank matrix completion using alternating minimization. In Proceedings of the Forty-Fifth Annual ACM Symposium on Theory of Computing, pages 665–674, 2013. [18] Raghunandan H. Keshavan, Andrea Montanari, and Sewoong Oh. Matrix completion from a few entries. IEEE Transactions on Information Theory, 56(6):2980–2998, 2010. [19] Christian Kümmerle and Claudio M. Verdun. A scalable second order method for illconditioned matrix completion from few samples. In Proceedings of the 38th International Conference on Machine Learning, pages 5872–5883, 2021. [20] Zhang Liu and Lieven Vandenberghe. Interior-point method for nuclear norm approximation with application to system identification. SIAM Journal on Matrix Analysis and Applications, 31(3):1235–1256, 2009. [21] Cong Ma, Kaizheng Wang, Yuejie Chi, and Yuxin Chen. Implicit regularization in nonconvex statistical estimation: Gradient descent converges linearly for phase retrieval, matrix completion, and blind deconvolution. Foundations of Computational Mathematics, 20(3):451–632, 2020.
58
[22] Thanh T. Ngo and Yousef Saad. Scaled gradients on Grassmann manifolds for matrix completion. In Advances in Neural Information Processing Systems, volume 25, pages 1412–1420, 2012. [23] Ruoyu Sun and Zhi-Quan Luo. Guaranteed matrix completion via non-convex factorization. IEEE Transactions on Information Theory, 62(11):6535–6579, 2016. [24] Jared Tanner and Ke Wei. Normalized iterative hard thresholding for matrix completion. SIAM Journal on Scientific Computing, 35(5):S104–S125, 2013. [25] Tian Tong, Cong Ma, and Yuejie Chi. Accelerating ill-conditioned low-rank matrix estimation via scaled gradient descent. Journal of Machine Learning Research, 22(150):1–63, 2021. [26] Joel A. Tropp. User-friendly tail bounds for sums of random matrices. Foundations of Computational Mathematics, 12(4):389–434, 2012. [27] Eilon Vaknin Laufer and Boaz Nadler. RGNMR: A Gauss–Newton method for robust matrix completion with theoretical guarantees. In Advances in Neural Information Processing Systems, volume 38, pages 66931–66973, 2025. [28] Bart Vandereycken. Low-rank matrix completion by Riemannian optimization. SIAM Journal on Optimization, 23(2):1214–1236, 2013. [29] Lingxiao Wang, Xiao Zhang, and Quanquan Gu. A unified computational and statistical framework for nonconvex low-rank matrix estimation. In Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, pages 981–990, 2017. [30] Tianming Wang and Ke Wei. Leave-one-out analysis for nonconvex robust matrix completion with general thresholding functions. Numerical Algorithms, 2026. [31] Ke Wei, Jian-Feng Cai, Tony F. Chan, and Shingyu Leung. Guarantees of Riemannian optimization for low rank matrix completion. Inverse Problems and Imaging, 14(2):233–265, 2020. [32] Xingyu Xu, Yandi Shen, Yuejie Chi, and Cong Ma. The power of preconditioning in overparameterized low-rank matrix sensing. In Proceedings of the 40th International Conference on Machine Learning, pages 38611–38654, 2023. [33] Xiaojing Zhu and Fengyi Yuan. A Riemannian regularized Gauss–Newton method for low-rank matrix completion. AIMS Mathematics, 10(12):28556–28582, 2025. [34] Pini Zilber and Boaz Nadler. GNMR: A provable one-line algorithm for low rank matrix recovery. SIAM Journal on Mathematics of Data Science, 4(2):909–934, 2022.
59