arXiv CSopen access

Improved Guarantees for Heterogeneous Treatment-Effect Estimation via Matrix Completion

Improved Guarantees for Heterogeneous Treatment-Effect Estimation via Matrix Completion

arXiv:2605.30319v1 [stat.ML] 28 May 2026

Anay Mehrotra Stanford University

Phuc Tran Vin University

Van H. Vu The University of Hong Kong

Manolis Zampetakis Yale University

Abstract A central goal of modern causal inference is estimating heterogeneous treatment effects to answer questions like “how does an intervention affect each unit,” rather than only on average. We study this problem with panel-data where we observe n units across m times under unknown, non-uniform treatment assignments. The data in this setting is naturally represented as a matrix of all unit–time treatment effects. Estimating heterogeneous treatment effects can then be expressed as obtaining a good estimation of each row’s average in this matrix. This allows us to formulate the problem as matrix completion, which can be solved under natural low-rankness assumptions. However, existing matrix-completion guarantees are not powerful enough to get meaningful bounds for the per-row guarantee required for estimating the heterogeneous treatment effect; roughly speaking, they are only useful for estimating average treatment effect bounds, as also illustrated in a recent line of work. We give a simple, computationally efficient estimator that, without knowledge of the propensities and under standard low-rankness and √ e ( 1/n + n/(m2 )). Technically, our regularity assumptions, achieves a row-wise ℓ2 error of O analysis establishes the first sharp row-wise ℓ2 -perturbation bound for low-rank approximation, complementing existing spectral-, Frobenius-, and entrywise perturbation theory.

1

Introduction

A central question in modern causal inference is to understand for whom and when a treatment works, not just whether it works on average [AI16]. As Kravitz, Duan, and Braslow [KDB04] put it, “the benefit or harm of [average] treatments in clinical trials can be misleading and fail to reveal the potentially complex mixture of substantial benefits for some, little benefit for many, and harm for a few.” Estimating such heterogeneous effects requires data that observes the same units (representing, e.g., users, patients, customers, or regions) across many randomized occasions, and a range of modern experimental designs produces data of precisely this form, including mobile-health trials, sequential A/B tests, contextual bandit experiments, and switchback experiments [LCLS10; JKPW22; BSZ23]. In all of these settings, we observe a panel of n units across m times, where m is typically larger than n (i.e., m ≥ n). Each unit–time pair (i, j) has two potential outcomes: Yij (1) which arises under treatment and Yij (0) which arises in the absence of treatment, i.e., in control. We observe only one outcome, the one corresponding to the assignment that was actually

made [IR15]. Stacking across units and times yields two outcome matrices Y (1), Y (0) ∈ Rn×m in
which exactly one entry from each pair Yij (1), Yij (0) is observed. The other entry is missing and is the counterfactual we would like to estimate. We focus on the randomized setting, in which the treatment indicator Dij ∈ {0, 1} is an independent Bernoulli draw with probability pij ∈ (0, 1). These probabilities can vary across units and times, reflecting different protocols, contexts, or eligibility rules. Further, they are, in general, unknown to the analyst: this is because, in many logged or retrospective experiments the assignment matrix is recorded but the platform’s specific traffic splits, randomization rules, and ramp-up schedules are not [SJ15; SSSC+16]. A canonical example is the family of mobile-health trials that randomize each participant many times over the course of an experiment, known as micro-randomized trials [KHSB+15; LKTM16; BAWM18; QWCK+22]. As a concrete example, consider the HeartSteps study [KSSL+19]. Here, researchers followed n = 44 participants for six weeks, randomizing each at up to five points per day among three options (no notification, a walking suggestion, or a stand-up suggestion) and recorded the participant’s step count over the following thirty minutes as the outcome; this yielded m = 210 times per participant. The average treatment effect, aggregated across all participants and notification types, was not statistically significant. However, a finer analysis revealed striking heterogeneity: for instance, the lift from walking suggestions, which initially more than doubled step count, decayed substantially over the six weeks, while stand-up suggestions had no detectable effect throughout. Such findings translate directly into actionable interventions, for example, reducing the frequency of walking suggestions as a trial progresses, and they are visible only through estimates of heterogeneous treatment effects. Signal–Noise Model. A substantial line of work in causal inference is devoted to estimating heterogeneous treatment effects [AI16; CCDD+18; WA18; KSBY19]. To formalize the treatment effect, we adopt the standard signal-plus-noise model Y ( a) = A( a) + E( a) for a ∈ {0, 1}, where E[ E( a)] = 0, so that A( a) := E[Y ( a)] is the expected potential outcome under action a for each unit at each time. The expectation is taken over noise from the environment and other unobserved factors. The natural object for heterogeneous-effect estimation is then the unit–time treatment-effect matrix M := A(1) − A(0) = E[Y (1)] − E[Y (0)] , whose row Mi = ( Mi1 , Mi2 , . . . , Mim ) records how unit i’s response to the treatment varies over time. For each unit i, one can define several heterogeneous treatment effects by computing averages Avgi (S) := ∑ j∈S Mij / |S| over subsets S ⊆ [m]. For instance, if S = [m], then Avgi (S) is the treatment effect specifically for unit i averaged over all times. One can also select S to be other meaningful subsets such as the first week of a trial, weekends, or high-engagement days. Matrix completion for causal inference. Estimating the matrix M naturally decomposes into two estimation problems, one for E[Y (1)] and one for E[Y (0)]. Fix any a ∈ {0, 1}: estimating E[Y ( a)] is a matrix-completion problem, where the goal is to recover the missing entries from those observed. Without further structure on the signal matrices A(0) and A(1) this task is impossible; however, real-world matrices have natural structural constraints that make the task tractable. In particular, Athey, Bayati, Doudchenko, Imbens, et al. [ABDI+21] proposed utilizing (approximate) low-rankness of the potential output matrices. Intuitively, low-rankness amounts to assuming that unit responses are governed by a small number of latent factors (such as baseline activity level,

receptiveness to prompts, time-of-day patterns, habituation, or seasonality). This assumption is well grounded in the causal panel data literature, where a long line of work builds on precisely this latent-factor structure [AG03; ADH10; Xu17; AAHI+21; BN21]. Casting the problem as matrix completion has the further benefit of bringing to bear a broad and well-developed algorithmic toolkit, e.g., [CR09; CT10; MHT10]. Limitations of prior guarantees. Following Athey, Bayati, Doudchenko, Imbens, et al. [ABDI+21], a growing line of work has applied matrix completion to treatment effect estimation under lowrankness and other standard regularity requirements [Xu17; AAHI+21; BN21; ADSS23]. The guarantees produced by these works, however, are not strong enough for our purpose: √ b − M∥F / nm, as in Athey, Bayati, • They either control the error in Frobenius norm, ∥ M Doudchenko, Imbens, et al. [ABDI+21]. This is insufficient to control the above heterogeneous effects Avgi (S) because, e.g., a single unit i with constant row-normalized error √ b i − Mi ∥ / m = Ω(1) contributes only O(1/√n) to the normalized Frobenius error, and so ∥M 2 leaves Avgi (S) uncontrolled at that unit. b ij − Mij |, as in Agarwal, Dahleh, Shah, • Or they control the worst-case entrywise error maxi,j | M and Shen [ADSS23]. While this is strong enough to bound Avgi (S), the best known bounds on √ √ this entrywise error decay only as 1/ log m + 1/ log n [ADSS23], which is prohibitively slow.1 Indeed, achieving an error of, e.g., 0.1 requires an astronomically large n, m ≳ e100 . b be an Thus, in our setting, none of the existing guarantees yields a meaningful bound. Let M estimator of M. A natural way to estimate Avgi (S) for each unit i and subsets S of interest is to satisfy: √ b ∥2,∞ := (1/√m) · max1≤i≤n ∥ Mi − M b i ∥2 ≤ poly(1/n, 1/m) . (1/ m) · ∥ M − M (1.1) d i (S) := ∑ j∈S M b ij / |S|. This immediately gives bounds on the estimation error of Avgi (S). Let Avg Now Cauchy–Schwarz inequality shows that, for every unit i and subset S ⊆ [m], p p d i (S) − Avg (S)| ≤ b ∥2,∞ ≤ 1/|S| · ∥ M − M m/|S| · poly(1/n, 1/m) . |Avg i d i (S) − Avg (S)| ≤ poly(1/n, 1/m). Hence, for any unit i and every subset S with |S| = Ω(m), |Avg i This leads us to the central question in the paper: Can such a row-wise guarantee be achieved from the observed entries, under existing structural assumptions? A bit more formally: b that, under existing Question. Given the observed outcomes, can one efficiently compute a matrix M √ b ∥2,∞ ≤ poly(1/n, 1/m)? low-rankness and regularity conditions, satisfies (1/ m)∥ M − M

1.1

Our contributions √

We answer the above question affirmatively, providing an estimator that achieves (1/ m)∥ M − √ b ∥2,∞ ≲ 1/n + 1/m with high probability under standard low-rankness and regularity conditions. M Concretely, we work under the following regularity conditions on the signal A and noise E. Informal Assumption 1.1 (Signal and noise regularity). For each action a ∈ {0, 1}, let σ1 ( a) ≥ σ2 ( a) ≥ . . . denote the singular values of A( a), and let U ( a), V ( a) collect the left and right singular vectors corresponding to its leading r singular components. The following hold for each a ∈ {0, 1}: 1 They incur this slow rate because their rate depends on the size of the largest all-1s sub-matrix of D which, in our

setting, has sides O(

p

log n) and O(

p

log m) with high probability.

√ • (Approximate low rankness) A( a) has at most r “large” singular values: σr+1 ( a) ≲ K m + n. • (Bounded signal and noise) The signal and noise are bounded: ∥ A( a)∥∞ , ∥ E( a)∥∞ ≤ K. • (Independent and mean-zero noise) The entries of E( a) are independent and mean-zero. • (Row and column incoherence) The leading left and right singular spaces are “incoherent:”

√

n max1≤i≤n ∥ei⊤ U ( a)∥2 ≤ µ

√

and

m max1≤ j≤m ∥e⊤ j V ( a )∥2 ≤ µ .

√ • (Signal-to-noise ratio) The first singular value of A is “large:” σ1 ( a) ≳ rK n + m. As mentioned, approximate low-rankness is standard and well motivated [ABDI+21]. Further, boundedness and mean-zero noise are mild regularity requirements satisfied in most practical settings. The incoherence and signal-to-noise assumptions are regularity assumptions that can also be shown to be necessary. Row incoherence rules out the degenerate case in which the signal is concentrated on a few units; and, due to this, the information from these units’ rows is not reliably available from the rest of the rows. Column incoherence is the analogous requirement across time. Finally, the lower bound on the signal-to-noise ratio is mild and needed for the algorithm to reliably estimate the rank r from the observed data; it is widely used in matrix perturbation work [SS90]. We need some additional notation to state our result. For each action a ∈ {0, 1}, define pij ( a) := Pr( Dij = a), so that pij (1) = pij and pij (0) = 1 − pij . Let pi ( a) := (1/m) ∑m j=1 pij ( a ) be the average probability of observing Y ( a) for unit i. We define q := mina,i pi ( a),

r p := maxa,i,j

pij ( a) , pi ( a )

P( a)ij :=

and

pij ( a) −1. pi ( a )

(1.2)

Here, q is the smallest average observation rate for any unit i and action a. The parameter r p measures how uneven the observation probabilities can be within a row after normalizing by the row average. Finally, P( a) measures the remaining within-row nonuniformity of the design. To gain some intuition, consider the special case where all observation probabilities are constant, pij = c. In this case, q = c, r p = 1, and P( a) = 0 for each a. Theorem 1.2 (Main guarantee, informal). Under Informal Assumption 1.1, there is a polynomial-time algorithm which, given only the observed outcomes Yobs and the assignment matrix D, and without any e satisfying, with high probability, knowledge of the propensities pij , outputs an estimate M 1 e Kr3/2 µ b ∥2,∞ ≤ O √ · ∥M − M m

"s

rp q



1 n + 2 n m



maxa ∥ P( a)∥op +p m · min{m, n}

#! .

In the simplest case where pij = c (for each i, j) and K, r, µ = O(1), the main theorem yields ! r 1 1 n e b ∥2,∞ ≤ O √ · ∥M − M + , n m2 m which achieves the desired poly(1/n, 1/m) rate in the standard regime where m ≳ n. In contrast, in the same setting the prior works either do not provide any non-trivial bound [ABDI+21] or only achieve a poly(1/log n, 1/log m) rate [ADSS23].

Another useful special case is the row-homogeneous design, where pij = pi for all j. This models experiments in which different units may be randomized at different rates, but each unit’s randomization rate is stable over time. Here, we obtain the following corollary Corollary 1.3 (Row-homogeneous propensities, informal). Suppose Informal Assumption 1.1 holds, √ pij = pi with pi ∈ [Ω(1), 1 − Ω(1)], and K, r, µ = O(1). Then with high probability, (1/ m) · ∥ M − √ e ( 1/n + n/m2 ). b ∥2,∞ ≤ O M Technical Novelty (also see Section 4). One important advantage of the estimator that we propose is that it is very simple and easy to implement: it is a simple row-scaled truncated-SVD algorithm (Algorithm 1). The analysis of this simple algorithm in our setting though, requires several new tools in the spirit of modern results from matrix perturbation theory. Existing work on perturbation bounds of low-rank approximations has primarily focused either on bounds for different norms e.g., [EY36; DTTZ14; Cha15; BV24; MV25; TV25b; TVV25; TV25a], which are not applicable in our setting, or bounds that do not gracefully improve for low-rank matrices [CTP19; AFWZ20], and hence, are not useful for analyzing our algorithm. Our main technical contribution is a perturbation theory for truncated SVD in the ∥·∥2,∞ norm, which directly controls the reconstructed error of our algorithm. To obtain the required sharp bound, we need to carefully adapt the contour expansion method, introduced in [TV25b; TVV25; TV26a] and applied in the more challenging norm ∥·∥2,∞ . In particular, the contour bootstraping argument used in [TV25b; TVV25; TV26a] does not apply in our case and we need a new idea to control the error that arises from the contour expansion method. We refer the reader to Section 4 for more details. A second difficulty comes from the fact that the probabilities pij are unknown and nonuniform. If the entrywise probabilities pij ( a) were known, one could form an unbiased inverse-propensityweighted matrix by scaling each observed entry by 1/pij ( a). Our algorithm does not have this information and instead scales row i by its empirical observation frequency. The population analogue of this row scaling is unbiased when the design is row-homogeneous (i.e., pij = pik for each j, k), but under nonuniform propensities its expectation satisfies pij ( a) Aij ( a)/pi ( a) = Aij ( a) + Pij ( a) Aij ( a). Thus the perturbation is not purely mean-zero noise: it also contains a deterministic bias induced by nonuniformity within one row. Our analysis incorporates this bias directly into the perturbation argument, yielding a bound whose additional design-dependent p term is governed by ∥ P( a)∥op / m · min(m, n), which decays as m, n → ∞ and vanishes when P( a) = 0.

1.2

Additional Related Works

Our work is broadly related to two lines of work: (1) work utilizing low-rank approximation results in causal inference and (2) work on matrix perturbation bounds for low-rank approximations. We have already mentioned key works from both lines in the introduction. Here we give a more detailed comparison with the first line. A more technical comparison with the second line appears in Technical Overview (Section 4). A growing literature uses low-rank structure to estimate missing counterfactual outcomes in panel data. As we have already mentioned, closest to our work are Athey, Bayati, Doudchenko, Imbens, et al. [ABDI+21] and Agarwal, Dahleh, Shah, and Shen [ADSS23], and their guarantees are insufficient for the heterogeneous effects we study. Several other works also exploit low-rankness in panel outcomes but their goal is to estimate (different

types of) average treatment effects [ASS18; AMSS19; AAHI+21; BN21; FFW21]. In contrast, we estimate heterogeneous treatment effects. A different line of work imposes low-rankness on the covariate matrix (rather than on the potential-outcome matrices) and use it to impute or de-noise covariates [ASSS21; ASS25a], while they also use low-rankness, their estimation targets and setting is quite different from our work. Finally, Agarwal, Shah, and Shen [ASS25b] extend the low-rank framework to settings with more than two treatment conditions and give row-wise guarantees for estimating Avgi ([m]). Our row-wise control of M is stronger: it yields bounds on Avgi (S) for every sufficiently large S ⊆ [m], not just S = [m].

2

Preliminaries

In this section, we introduce basic notation and preliminaries. Notation. For constants N and 1 ≤ k ≤ N, let e N,k be the kth standard basis vector in R N ; when p the dimension is clear, we write ek . For a vector v and p ≥ 1, define ∥v∥ p as (∑ j v j )1/p and ∥v∥∞ as max j v j . For a matrix B ∈ Rn×m , we write Bi for its ith row, B· j for its jth column, ∥ B∥∞ for maxi,j Bij , and ∥ B∥op and ∥ B∥F for its operator and Frobenius norms respectively. The central object for our row-wise guarantees is the row-wise ℓ2 -norm ∥ B∥2,∞ := max1≤i≤n ∥ Bi ∥2 . Potential Outcomes and the Heterogeneous Treatment Effect. We observe n units (e.g., users or patients) across m times (e.g., decision points or days). For each unit–time pair (i, j) and action a ∈ {0, 1}, Yij ( a) denotes the potential outcome of unit i at time j under action a. Following the standard potential-outcomes framework [IR15; HR23], we adopt the signal-plus-noise model Y ( a) = A( a) + E( a) with E[ E( a)] = 0, so that A( a) := E[Y ( a)] is the mean potential-outcome matrix. Our goal is to estimate the heterogeneous treatment-effect matrix M := A(1) − A(0), whose entry Mij is the mean effect of treating unit i at time j. Unless otherwise stated, expectations are taken over the outcome noise E, not the treatment assignments. Assignment and Observed Data. We make standard assumptions on the assignment mechanism [IR15; HR23]. For each (i, j), the treatment indicator Dij ∈ {0, 1} is an independent Bernoulli draw, Dij ∼ Ber( pij ), where the probabilities pij may vary across units and times and are not assumed to be known. The assignments D = ( Dij ) are independent of the potential outcomes {Y (0), Y (1)}, and the observed outcome at (i, j) is Yijobs := Dij Yij (1) + (1 − Dij )Yij (0). To treat the two actions symmetrically, we introduce, Dij (1) := Dij ,

Dij (0) := 1 − Dij ,

pij (1) := pij ,

pij (0) := 1 − pij ,

so that Dij ( a) = 1 when Yij ( a) is observed and E[ Dij ( a)] = pij ( a). The corresponding partially e( a) is defined entrywise by observed matrix Y eij ( a) := Yij ( a) Y

if

Dij ( a) = 1

and

eij ( a) := ⋆ Y

otherwise ,

with ⋆ denoting a missing entry. The analyst equivalently sees the pair ( D, Y obs ), or the pair of e (0), Y e(1)) together with their observation masks. partially observed matrices (Y Row Propensities and Non-uniformity. The rates we obtain depend on the parameters q, r p , and P( a) defined in Equation (1.2). The parameter q controls how often the rarest action is

observed at the rarest unit. Requiring q ≥ Ω(1) is substantially milder than the standard overlap condition [IR15; HR23], which requires every pij to be bounded away from 0 and 1; our condition instead constrains only the row averages, leaving individual pij free to be arbitrarily close to 0 or 1 at many time points within a unit. The parameters r p and P( a) quantify within-row nonuniformity. By construction, each row of P( a) averages to zero, and P(0) = P(1) = 0 exactly when the design is row-homogeneous, i.e., pij = pi for all i, j. SVD, Low-Rankness, and Noise. For a rank-r matrix A ∈ Rn×m , write A = UΣV ⊤ = n×r and V = [ v , . . . , v ] ∈ Rm×r have orthonormal ∑rℓ=1 σℓ uℓ v⊤ r 1 ℓ , where U = [ u1 , . . . , ur ] ∈ R columns and the singular values σ1 ≥ σ2 ≥ · · · ≥ σr are arranged in non-increasing order. We write As for the best rank-s approximation of A, and δs ( A) := σs − σs+1 for its sth singular value gap; when A is clear from context we abbreviate δs := δs ( A). The rowand column-incoherence parameters of A are µ R ( A) :=

√

⊤ U ∥2 n max1≤i≤n ∥en,i

and

µC ( A) :=

√

⊤ V ∥2 , m max1≤ j≤m ∥em,j

(2.1)

and we set µ( A) := max{µ R ( A), µC ( A)}. These quantities are small when the singular spaces spread across units and time points, and large when either side concentrates on few coordinates. When both potential-outcome signals are under consideration we write µ := maxa∈{0,1} µ( A( a)). Definition 2.1 ((K, σ)-Bounded Random Matrix). A random matrix ER ∈ Rn×m is called (K, σ)bounded if its entries satisfy E ER,ij = 0, E[| ER,ij |2 ] ≤ σ2 , and E[| ER,ij |ℓ ] ≤ K ℓ−2 σ2 for every ℓ ≥ 2 and all i ∈ [n], j ∈ [m]. This condition fixes the variance scale and supplies the higher-moment control used in Bernsteintype matrix concentration [Tro12; Tro15].

3

Our Main Results

In this section, we present our main results. We begin with the formal version of Informal Assumption 1.1, followed by our estimator, and then its guarantees. Recall the parameters q, r p , and P( a) from Equation (1.2). For each a ∈ {0, 1}, write s m + rpn T ( a) := log(m + n) + ∥ P( a)∥op . q We impose the following standard regularity conditions Assumption 3.1 (Regularity conditions). Fix constants K A , KE , µ > 0 and r ∈ N. Set K := K A + KE . For each a ∈ {0, 1}, the following conditions hold: √ 1. (Low-Rankness) ∥ A( a) − Ar ( a)∥op ≲ K m + n, where Ar ( a) is A( a)’s best rank-r approximation 2. (Bounded Signal and Noise) ∥ A( a)∥∞ ≤ K A and E is (KE , KE )-bounded (Definition 2.1). 3. (Independent and mean-zero) E( a)’s entries are independent of each other and D, and mean-zero

4. (Row and column incoherence) µC ( A( a)), µ R ( A( a)) ≤ µ, where µ(·) is defined in (2.1). 5. (Signal-to-noise ratio) The leading singular value satisfies σ1 ( a) ≳ KrT ( a). We refer the reader to Section 1.1 and Section A.1 for discussion of this assumption. b(0) and A b(1) separately and Next, we present our estimator (see Algorithm 1). It computes A b b b := A(1) − A(0). To estimate A( a), we fill unobserved entries with zeros, rescale each returns M row by its empirical observation frequency, and keep the largest singular block separated by a spectral gap. Algorithm 1 Row-scaled spectral estimator 1: Input: Observed outcomes Y obs , matrix D, r from Assumption 3.1, and thresholds τ0 , τ1 . 2: for each a ∈ {0, 1} do 3: 4: 5: 6: 7: 8:

−1 Initialize the estimates pbi ( a) ← max{m−1 ∑m j=1 Dij ( a ), m } for all i ∈ [ n ]. eub ( a) ← Zij ( a)/ pbi ( a). For each (i, j) ∈ [n] × [m], set Zij ( a) ← Dij ( a)Y obs and Y

ij ij ub ub ub e Let e σ1 ( a) ≥ e σ2 ( a) ≥ . . . be the singular values of Y ( a). Let S a ← {1 ≤ s ≤ r : e σsub ( a) − e σsub +1 ( a ) ≥ τa }.

Set b s( a) ← max S a if S a ̸= ∅, and b s( a) ← 0 otherwise. ub b e Set A( a) ← [Y ( a)]bs(a) , where [ B]s denotes the best rank-s approximation of B

b(1) − A b(0). b ←A 9: Return M Running time of Algorithm 1. The dominant cost is computing the top r + 1 singular values eub ( a) for each a ∈ {0, 1}, which suffice to determine b eub ( a)]bs(a) . This and vectors of Y s( a) and [Y e (nmr ) time via Lanczos’s method or randomized SVD [GV13]. All remaining steps run in takes O e (nmr ). O(nm + nr ) time, giving a total running time of O Main Result. Next, we state our main result, which bounds the error in estimation of A(0) and A(1); combining the two bounds immediately implies a bound on the estimation error for M. Theorem 3.2 (Upper bound on Error). Suppose Assumption 3.1 holds. Run Algorithm 1 with thresholds b(0) τa = 96KT ( a) for a ∈ {0, 1}. Then, with probability at least 1 − O(1/(m+n)), the output matrices A b(1) satisfy, for each a ∈ {0, 1}, and A "r # P ( a )∥ ∥ √ r r op p p 3/2 b( a) − A( a)∥ + + √ . µ ( m + n) log4 (m + n) ∥A 2,∞ ≲ Kr mq nq mn Now, applying Theorem 3.2 to both a = 0, 1 and using the triangle inequality yields the following b row-wise guarantee for M. Corollary 3.3 (Row-wise recovery of the treatment-effect matrix). Under the assumptions of Theorem 3.2, with probability at least 1 − O(1/(m+n)), "s  #  b i − Mi ∥ maxa∈{0,1} ∥ P( a)∥op rp 1 ∥M n 4 3/2 2 √ ≲ Kr µ log (m + n) max + + p . q n m2 1≤ i ≤ n m m min{m, n}

Thus, when m ≥ n, K, r, µ, r p and q−1 are constants, and the design is row-homogeneous (so e (n−1/2 ). Next, to build some intuition, we highlight the P(0) = P(1) = 0), the error scales as O bounds we obtain under two regimes that are common in panel-data applications. Special Case I (Row-homogeneous design): Here, pij may differ across units i but, for each i, is invariant in j: pij = pi for every i ∈ [n] and j ∈ [m]. (3.1) This setting captures stratified and covariate-adaptive randomization protocols, and also includes the Bernoulli design pij = c as a sub-case. It is useful because under any row-homogeneous design, P(0) = P(1) = 0 and r p = 1, leading to the following bound. Corollary 3.4 (Row-homogeneous design). Consider the row-homogeneous design above (Equation (3.1)). Suppose the assumptions of Theorem 3.2 hold, K, r, µ = O(1), q = Ω(1), and m ≥ n. Then, with √
e 1 √ b − M∥ probability 1 − O(1/(m+n)) it holds that 1/ m · ∥ M 2,∞ ≤ O ( / n). In the balanced regime m ≍ n, this matches the best known guarantees for matrix completion under the Bernoulli design, despite the fact that our algorithm does not have access to the propensities pi . Special Case II (Spectrally small within-row non-uniformity): Our next regime is more general: it allows pij to vary in time within rows but requires the variation to be spectrally small. Concretely, √ √ e ( m + n) . for some ν ≥ 0 , maxa∈{0,1} ∥ P( a)∥op ≤ ν O (3.2) For example, this holds whenever P( a)’s entries are independent and sub-Gaussian, as in contextual-bandit-style adaptive designs in which propensities depend sufficiently mildly on observed covariates; standard random-matrix bounds then give the spectral condition with high probability. Corollary 3.5 (Spectrally small within-row nonuniformity). Consider the setting above (Equation (3.2)). Suppose the assumptions of Theorem 3.2 hold, K, r, µ = O(1), q = Ω(1), and m ≥ n. Then, with √ e −1/2 ). b − M∥ probability at least 1 − O(1/(m+n)), it holds that (1/ m) · ∥ M 2,∞ ≤ ν O ( n e (n−1/2 ) row-normalized rate of the row-homogeneous case, Thus, Corollary 3.5 preserves the O with ν inflating only the constant. Crucially, the algorithm itself does not change between the two regimes; the analyst does not need to verify the spectral bound or estimate ν.

4

Technical Overview

In this section, we sketch the key ideas behind the proof of Theorem 3.2. In this section, we focus on a = 1 and, hence, omit it from the notation, writing, e.g., A and P for A(1) and P(1) respectively. Now, our goal is to prove that with probability at least 1 − O(1/(m+n)), r  √ √ ∥ P∥op rp rp 4 3/2 b − A∥2,∞ ≲ Kr µ log (m + n)( m + n) ∥A + + √ . (4.1) mq nq mn Where the constants K, r, and µ are from Assumption 3.1 and the parameters (rq , q, P) are as defined in Equation (1.2). Before presenting the proof, we need to set up some notation.

eub is from Algorithm 1. A eub , where Y Notation and basic observations. We use Aub to denote E Y
direct computation shows Aub := pij aij /pi ij . For any X ∈ Rn×m , we denote the corresponding

e = xeij and the corresponding “scaled” matrix by X e ub := xeij /pi . With “observed” matrix by X ij

ij

eub = A eub + E eub . Finally, we define this notation, since Y = A + E, we can write Y    pij eub − Aub ) + E eub . E0 := Aub − A = and ER := (Y − 1 aij pi ij Next, we observe that Assumption 3.1(5) and the choice that e σsub − e σsub +1 ≳ KT ( a ) are equivalent ub ub σs − e σs+1 ≥ 8(∥ ER ∥op + ∥ E0 ∥op ) respectively. (Where we to σ1 > 10r (∥ ER ∥op + ∥ E0 ∥op ), and e selected the constants to simplify exposition.) b − A∥2,∞ into two parts: Our approach. At a high-level, to prove (4.1), we split ∥ A 1. Term 1 (Tail-error of low-rank approximation):

∥ As − A∥2,∞ ; and

2. Term 2 (Perturbation of low-rank approximations):

b − As ∥2,∞ ∥A

One can bound the first term, ∥ As − A∥2,∞ , by combining the singular decomposition of A and the definition of the threshold position s. The key difficulty is bounding the second term, which measures the perturbation of low-rank approximations in ∥ · ∥2,∞ . Concretely, the triangle inequality implies b − A∥2,∞ = ∥( A b − As ) + ( As − A)∥2,∞ ≤ ∥ A b − As ∥2,∞ + ∥ As − A∥2,∞ . ∥A Step 1 (Bounding Term 1). Since A = ∑ri=1 σi ui vi⊤ , A − As = ∑ri=s+1 σi ui vi⊤ . Therefore,

∥ As − A∥2,∞ =

⊤ max ∥en,k ( As − A)∥2

1≤ k ≤ n

max ∥ ∑r≥i>s σi uik vi⊤ ∥2 q r = max ∑i=s+1 σi2 u2ik 1≤ k ≤ n q r ≤ max ∥ui ∥∞ ∑i=s+1 σi2 s <i ≤r √ ≤ max ∥ui ∥∞ rσs+1 .

=

1≤ k ≤ n

s <i ≤r

We claim that σs+1 ≤ 10r (∥ ER ∥op + ∥ E0 ∥op ). Suppose, for contradiction, that σs+1 > 10r (∥ ER ∥op + ∥ E0 ∥op ). Since rank( A) ≤ r (i.e., σr+1 = 0), it follows that s + 1 ≤ r. Arguing as in Remark A.2 (with σs+1 in place of σ1 ), there exists some j with s + 1 ≤ j ≤ r such that e σjub − e σjub +1 > 8(∥ ER ∥op + ∥ E0 ∥op ). This contradicts the definition of s as the largest index satisfying e σsub − e σsub +1 > 8(∥ ER ∥op + ∥ E0 ∥op ). Hence, σs+1 ≤ 10r (∥ ER ∥op + ∥ E0 ∥op ). Therefore,

∥ As − A∥2,∞ ≤ 10r3/2 (∥ ER ∥op + ∥ E0 ∥op ) · maxs<i≤r ∥ui ∥∞ . To bound ∥ ER ∥, we use the following lemma. Its proof will be presented later in Section C.3. p p Lemma 4.1. With high probability, ∥ ER ∥op ≤ 12K log(m + n) · r p (m+n)/q.

√

√

We also have ∥ E0 ∥op ≤ K A ∥ P∥op . By the definition of µ, maxs<i≤r ∥ui ∥∞ ≤ µ(1/ m + 1/ n). Combining all estimates, we obtain  r q √ √
∥ P∥op rp rp 3/2 ∥ As − A∥2,∞ ≤ 12 log(m + n)Kr µ m+ n · + + √ . (4.2) mq nq mn Step 2 (Bounding Term 2). Our estimate is based on the following theorem, which is a key part of our technical contribution. We compare its proof techniques with prior work in Section 4.1. The discussion of the perturbation of low-rank approximations in ∥ · ∥2,∞ and the detailed proof of Theorem 4.2 appear in Section B. e = A + ER + E0 , where ER is (K, σ)-bounded. There is a universal constant C > 0 Theorem 4.2. Let A satisfying: If δs ≥ 6(∥ ER ∥op + ∥ E0 ∥op ), then with probability 1 − O(1/(m+n)), !
 √ es − As ∥2,∞ σs m + n σ + ∥ E0 ∥op µ ∥A K log4 (m + n) µ 2 √ √ √ +√ ≤ log (m + n) + . δs C r m n m+n eub = A + ER + E0 , we have the following observations. Since Y • Since e σsub − e σsub +1 ≥ 8(∥ ER ∥op + ∥ E0 ∥op ), by Weyl’s inequality (Theorem A.3), we have σs − σs+1 ≥ e σsub − e σsub +1 − 2(∥ ER ∥op + ∥ E0 ∥op ) ≥ 6(∥ ER ∥op + ∥ E0 ∥op ). • By the definition of ER = ( ERij )ij , we have (

(1− pij ) aij/pi + ε ij/pi with probability pij ,

ERij =

− pij aij/pi with probability 1 − pij .

For each ℓ ≥ 2, a direct moment calculation gives h

E | ERij |

ℓ

i

≤

pij (K A + KE )ℓ piℓ

=

pij K ℓ piℓ

.

eub , A) with the noises ER , E0 . These observations allow us to apply Theorem 4.2 on the pair (Y p Indeed, the parameters K and σ in Theorem 4.2 are replaced respectively by K/q and r p K/q in our setting. Thus, with probability at least 1 − O((m + n)−1 ), there is a constant C > 0 such that b − As ∥2,∞ esub − As ∥2,∞ √ √ σs ∥A ∥Y √ √ = ≲ log4 (m + n)Kµ( m + n) · δs C r C r

r

∥ P∥op rp rp + + √ mq nq mn



r (∥ E ∥ +∥ E ∥ )

Since σs /δs = 1 + (σs+1 /δs ) ≤ 1 + 8(∥E R∥opop+∥E 0∥opop∥) ≤ 2r, we further have R

0

 r  √ √ ∥ P∥op rp rp 4 3/2 b √ ∥ A − As ∥2,∞ ≤ O log (m + n)Kr µ( m + n) · + + . mq nq mn b − As ∥2,∞ , ∥ As − A∥2,∞ , we finally obtain Combining all estimates (4.2), (4.3) on ∥ A   r √ √ ∥ P∥op rp rp 4 3/2 b − A∥2,∞ = O log (m + n)Kr µ( m + n) · ∥A + + √ . mq nq mn

(4.3)

4.1

High-level Proof Sketch of Theorem 4.2

To prove Theorem 4.2, we utilize a carefully designed adaptation of the contour integral approach that has been explored heavily in many recent works; e.g., [KX16; OVW23; TV25b; TV25c; TV26a]. es − As directly, we work with their symmetrized versions: Instead of analyzing A " # " # " # e 0 A 0 E 0 A A := , E := ⊤ , Ae := e⊤ . A⊤ 0 E 0 A 0 es − As by The symmetrized versions are related to A e es − As ∥2,∞ := max1≤k≤n ∥e⊤ ( A es − As )∥2 = max1≤k≤n ∥e⊤ ∥A n,k m+n,k ( A2s − A2s ) 2 . Thus, to prove Theorem 4.2 it suffices to bound ∥e1⊤ (Ae2s − A2s )∥2 . Next, using the Cauchy integral theorem (Theorem A.4), we obtain 1 e1⊤ (Ae2s − A2s ) = 2πi

R

Γ

e −1 )−1 − (zI − A)−1 ] dz, z · e1⊤ [(zI − A)

where Γ is a contour in C that encloses ±σ1 , ±σ2 , . . . , ±σs and excludes ±σs+1 , ±σs+2 , . . . , ±σr . Challenge (Prior techniques are insufficient to bound ∥·∥2,∞ norm). If instead of the ∥·∥2,∞ , we wanted to bound the ∥·∥2 or ∥·∥op norms, then one could straightforwardly bound ∥e1⊤ (Ae2s − A2s )∥2 by ∥Ae2s − A2s ∥op , and then apply the existing spectral-norm bounds on perturbations of low-rank approximations (e.g., [EY36; TVV25; TV25a]). However, to obtain meaningful bounds on the Heterogenous treatment effect, we need to focus on the ∥·∥2,∞ norm and, here, using the aforementioned bounds yields a suboptimal bound which can off by a large factor, of up to √ n + m. Thus, obtaining a sharp bound with respect to the row-wise ℓ2 -norm, which is crucial to prove Theorem 4.2, remains a formidable analytical challenge that requires some new ideas. Ideas. To obtain the sharp bound, we carefully adapt the contour expansion method, introduced in [TV25b; TVV25; TV26a]. In particular, we repeatedly apply the Sherman–Morrison–Woodbury e = A + E, to get formula M−1 − ( M + N )−1 = ( M + N )−1 N M−1 [HJ12] and A ⊤ e1⊤ (Ae2s − A2s ) = ∑∞ k =1 e1 Hk ,

where

1 Hk := 2πi

R

Γ

z · (zI − A)−1 [E (zI − A)−1 ]k dz.

Intuitively, [TVV25] used a contour bootstrapping argument to show that ∥Ae2s − A2s ∥op is of the same order as ∥ H1 ∥op . For entrywise control, [TV25b] bounded ∥ Hk ∥∞ for all k ≥ 1. In our row-wise ℓ2 setting, we instead bound ∥e1⊤ Hk ∥2 for 1 ≤ k ≤ ℓ log n, with a suitable constant ℓ, and show that the remaining tail is negligible. Further challenges. Bounding ∥e1⊤ Hk ∥2 presents several further challenges. Sharp estimates ⊤ ℓ ℓ require controlling the interactions e1⊤ E ℓ (U V ) , ∀1 ≤ ℓ ≤ k. Naive bounds based on ∥ e1 E ∥2 , ∥E ∥op , ℓ or E ℓ (U V ) 2 are suboptimal. Since E is random, one expects E to spread mass across the entries ⊤

⊤ ℓ U of (U V ) , making the incoherence parameter µ essential. In [TV25b], the authors handled e1 E ( V ) when E is random and mean-zero. In our setting, E = E R + E0 contains both random and deterministic components, so separating their contributions is highly nontrivial; see Section C.

5

Conclusion

In this work, we study heterogeneous treatment-effect estimation in panel experiments where each unit is randomized many times, the propensities are unknown and may vary across units and time, and the potential-outcome matrices are approximately low-rank. We propose a simple row-scaled truncated-SVD estimator that uses only the observed outcomes and the assignment matrix, and we show that it recovers each unit’s treatment-effect trajectory in a row-wise ℓ2 sense. This is enough to estimate treatment-effect averages for any individual unit over reasonably large subsets of times, which is the kind of guarantee one needs for the heterogeneous-effect questions that motivated the problem. At the heart of the analysis is a new perturbation bound for truncated SVD in the ∥·∥2,∞ norm, which also makes precise how nonuniformity in the design affects estimation. Many natural questions remain. It would be interesting to relax the low-rank, incoherence, and signal-to-noise conditions, and to handle adaptive or dependent assignment mechanisms that arise in sequential experiments. Finally, sharpening the dependence on the design-nonuniformity term, especially in highly heterogeneous designs, is an interesting direction as well.

References [AAHI+21]

[ABDI+21]

[ADH10]

[ADSS23]

[AFWZ20]

[AG03] [AI16] [AMSS19]

[ASS18]

Dmitry Arkhangelsky, Susan Athey, David A. Hirshberg, Guido W. Imbens, and Stefan Wager. “Synthetic Difference-in-Differences”. In: American Economic Review 111.12 (2021), pp. 4088–4118 (cit. on pp. iii, vi). Susan Athey, Mohsen Bayati, Nikolay Doudchenko, Guido Imbens, and Khashayar Khosravi. “Matrix Completion Methods for Causal Panel Data Models”. In: Journal of the American Statistical Association 116.536 (2021), pp. 1716–1730 (cit. on pp. ii–v). Alberto Abadie, Alexis Diamond, and Jens Hainmueller. “Synthetic Control Methods for Comparative Case Studies: Estimating the Effect of California’s Tobacco Control Program”. In: Journal of the American Statistical Association 105.490 (2010), pp. 493–505 (cit. on p. iii). Anish Agarwal, Munther Dahleh, Devavrat Shah, and Dennis Shen. “Causal Matrix Completion”. In: Proceedings of the Thirty Sixth Conference on Learning Theory. Vol. 195. Proceedings of Machine Learning Research. PMLR, 2023, pp. 3821–3826. url: https://proceedings. mlr.press/v195/agarwal23c.html (cit. on pp. iii–v). Emmanuel Abbe, Jianqing Fan, Kaizheng Wang, and Yiqiao Zhong. “Entrywise eigenvector analysis of random matrices with low expected rank”. In: Annals of statistics 48.3 (2020), p. 1452 (cit. on p. v). Alberto Abadie and Javier Gardeazabal. “The Economic Costs of Conflict: A Case Study of the Basque Country”. In: American Economic Review 93.1 (2003), pp. 113–132 (cit. on p. iii). Susan Athey and Guido Imbens. “Recursive Partitioning for Heterogeneous Causal Effects”. In: Proceedings of the National Academy of Sciences 113.27 (2016), pp. 7353–7360 (cit. on pp. i, ii). Muhammad Amjad, Vishal Misra, Devavrat Shah, and Dennis Shen. “mRSC: Multi-dimensional Robust Synthetic Control”. In: Proceedings of the ACM on Measurement and Analysis of Computing Systems 3.2 (2019), pp. 1–27 (cit. on p. vi). Muhammad Amjad, Devavrat Shah, and Dennis Shen. “Robust Synthetic Control”. In: Journal of Machine Learning Research 19.22 (2018), pp. 1–51. url: https://jmlr.org/papers/ v19/17-777.html (cit. on p. vi).

[ASS25a]

[ASS25b]

[ASSS21]

[BAWM18]

[BN21]

[BSZ23] [BV16]

[BV24]

[CCDD+18]

[Cha15] [CR09] [CT10]

[CTP19]

[DTTZ14]

[EY36] [FFW21]

[GV13] [Hig08]

Anish Agarwal, Devavrat Shah, and Dennis Shen. “On Model Identification and Out-ofSample Prediction of PCR with Applications to Synthetic Controls”. In: Journal of Machine Learning Research 26.117 (2025), pp. 1–58. url: https://jmlr.org/papers/v26/230102.html (cit. on p. vi). Anish Agarwal, Devavrat Shah, and Dennis Shen. “Synthetic Interventions: Extending Synthetic Controls to Multiple Treatments”. In: Operations Research 74.2 (2025), pp. 840–859 (cit. on p. vi). Anish Agarwal, Devavrat Shah, Dennis Shen, and Dogyoon Song. “On Robustness of Principal Component Regression”. In: Journal of the American Statistical Association 116.536 (2021), pp. 1731–1745 (cit. on p. vi). Audrey Boruvka, Daniel Almirall, Katie Witkiewitz, and Susan A. Murphy. “Assessing TimeVarying Causal Effect Moderation in Mobile Health”. In: Journal of the American Statistical Association 113.523 (2018), pp. 1112–1121 (cit. on p. ii). Jushan Bai and Serena Ng. “Matrix Completion, Counterfactuals, and Factor Analysis of Missing Data”. In: Journal of the American Statistical Association 116.536 (2021), pp. 1746–1763 (cit. on pp. iii, vi). Iavor Bojinov, David Simchi-Levi, and Jinglong Zhao. “Design and Analysis of Switchback Experiments”. In: Management Science 69.7 (2023), pp. 3759–3777 (cit. on p. i). Afonso S Bandeira and Ramon Van Handel. “Sharp nonasymptotic bounds on the norm of random matrices with independent entries”. In: Ann. Probab. 44.4 (2016), pp. 2479–2506 (cit. on pp. xxvi, xxviii). Abhinav Bhardwaj and Van Vu. “Matrix perturbation: Davis-Kahan in the infinity norm”. In: Proceedings of the 2024 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA) (2024), pp. 880–934 (cit. on p. v). Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins. “Double/Debiased Machine Learning for Treatment and Structural Parameters”. In: The Econometrics Journal 21.1 (2018), pp. C1–C68 (cit. on p. ii). Sourav Chatterjee. “Matrix Estimation by Universal Singular Value Thresholding”. In: The Annals of Statistics 43.1 (2015), pp. 177–214 (cit. on p. v). Emmanuel J. Candès and Benjamin Recht. “Exact Matrix Completion via Convex Optimization”. In: Foundations of Computational Mathematics 9.6 (2009), pp. 717–772 (cit. on p. iii). Emmanuel J. Candès and Terence Tao. “The Power of Convex Relaxation: Near-Optimal Matrix Completion”. In: IEEE Transactions on Information Theory 56.5 (2010), pp. 2053–2080 (cit. on p. iii). Joshua Cape, Minh Tang, and Carey E Priebe. “The two-to-infinity norm and singular subspace geometry with applications to high-dimensional statistics”. In: Annals of Statistics 47.5 (2019), pp. 2405–2439 (cit. on p. v). C. Dwork, K. Talwar, A. Thakurta, and L. Zhang. “Analyze Gauss: Optimal bounds for privacy-preserving principal component analysis”. In: Proceedings of the Forty-Sixth Annual ACM Symposium on Theory of Computing. 2014, pp. 11–20 (cit. on p. v). G. Eckart and G. Young. “The approximation of one matrix by another of lower rank”. In: Psychometrika 1 (1936), pp. 211–218 (cit. on pp. v, xii, xvii, xix, xxi). Iván Fernández-Val, Hugo Freeman, and Martin Weidner. “Low-Rank Approximations of Nonseparable Panel Models”. In: The Econometrics Journal 24.2 (2021), pp. C40–C77 (cit. on p. vi). Gene H. Golub and Charles F. Van Loan. Matrix Computations. 4th ed. Johns Hopkins University Press, 2013 (cit. on p. viii). Nicholas J. Higham. Functions of Matrices: Theory and Computation. See §1.3 and §7.6 for the Cauchy–Dunford integral. SIAM, 2008 (cit. on p. xx).

[HJ12] [HR23]

[IR15] [JKPW22]

[JSS22] [Kat80] [KDB04]

[KHSB+15]

[KSBY19]

[KSSL+19]

[KX16] [LCLS10]

[LKTM16]

[MHT10]

[MV25]

[OVW23]

[QWCK+22]

R. A. Horn and C. R. Johnson. Matrix Analysis. Cambridge University Press, 2012 (cit. on pp. xii, xxi). Miguel A. Hernan and James M. Robins. Causal Inference: What If. Chapman & Hall/CRC monographs on statistics & applied probability. Taylor & Francis, 2023. url: https:// books.google.co.in/books?id=FPkN0AEACAAJ (cit. on pp. vi, vii). Guido W. Imbens and Donald B. Rubin. Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction. Cambridge University Press, 2015 (cit. on pp. ii, vi, vii). Ramesh Johari, Pete Koomen, Leonid Pekelis, and David Walsh. “Always Valid Inference: Continuous Monitoring of A/B Tests”. In: Operations Research 70.3 (2022), pp. 1806–1821 (cit. on p. i). V. Jain, A. Sah, and M. Sawhney. “On the smoothed analysis of the smallest singular value with discrete noise”. In: Bull. Lond. Math. Soc. 54.2 (2022), pp. 369–388 (cit. on p. xxiv). Tosio Kato. Perturbation Theory for Linear Operators. Classics in Mathematics. Springer, 1980 (cit. on p. xx). Richard L. Kravitz, Naihua Duan, and Joel Braslow. “Evidence-Based Medicine, Heterogeneity of Treatment Effects, and the Trouble with Averages”. In: The Milbank Quarterly 82.4 (2004), pp. 661–687 (cit. on p. i). Predrag Klasnja, Eric B. Hekler, Saul Shiffman, Audrey Boruvka, Daniel Almirall, Ambuj Tewari, and Susan A. Murphy. “Microrandomized Trials: An Experimental Design for Developing Just-in-Time Adaptive Interventions”. In: Health Psychology 34.S (2015), pp. 1220– 1228 (cit. on p. ii). Sören R. Künzel, Jasjeet S. Sekhon, Peter J. Bickel, and Bin Yu. “Metalearners for Estimating Heterogeneous Treatment Effects Using Machine Learning”. In: Proceedings of the National Academy of Sciences 116.10 (2019), pp. 4156–4165 (cit. on p. ii). Predrag Klasnja, Shawna Smith, Nicholas J. Seewald, Andy Lee, Kelly Hall, Brook Luers, Eric B. Hekler, and Susan A. Murphy. “Efficacy of Contextually Tailored Suggestions for Physical Activity: A Micro-randomized Optimization Trial of HeartSteps”. In: Annals of Behavioral Medicine 53.6 (2019), pp. 573–582 (cit. on p. ii). V. Koltchinskii and D. Xia. “Perturbation of linear forms of singular vectors under Gaussian noise”. In: 71 (2016), pp. 397–423 (cit. on pp. xii, xx). Lihong Li, Wei Chu, John Langford, and Robert E. Schapire. “A Contextual-Bandit Approach to Personalized News Article Recommendation”. In: Proceedings of the 19th International Conference on World Wide Web. WWW ’10. ACM, 2010, pp. 661–670 (cit. on p. i). Peng Liao, Predrag Klasnja, Ambuj Tewari, and Susan A. Murphy. “Sample Size Calculations for Micro-randomized Trials in mHealth”. In: Statistics in Medicine 35.12 (2016), pp. 1944–1971 (cit. on p. ii). Rahul Mazumder, Trevor Hastie, and Robert Tibshirani. “Spectral Regularization Algorithms for Learning Large Incomplete Matrices”. In: Journal of Machine Learning Research 11.80 (2010), pp. 2287–2322 (cit. on p. iii). Oren Mangoubi and Nisheeth K. Vishnoi. “Private Low-Rank Approximation for Covariance Matrices, Dyson Brownian Motion, and Eigenvalue-Gap Bounds for Gaussian Perturbations”. In: J. ACM 72.2 (Mar. 2025). url: https://doi.org/10.1145/3716496 (cit. on pp. v, xix). Sean O’Rourke, Van Vu, and Ke Wang. “Matrices with Gaussian noise: Optimal estimates for singular subspace perturbation”. In: IEEE Transactions on Information Theory (2023) (cit. on pp. xii, xx). Tianchen Qian, Ashley E. Walton, Linda M. Collins, Predrag Klasnja, Stephanie T. Lanza, Inbal Nahum-Shani, Mashfiqui Rabbi, Michael A. Russell, Maureen A. Walton, Hyesun Yoo, and Susan A. Murphy. “The Micro-Randomized Trial for Developing Digital Interventions:

[SJ15]

[SS90] [SSSC+16]

[Tro12] [Tro15] [TV08] [TV25a] [TV25b]

[TV25c]

[TV25d]

[TV26a]

[TV26b] [TVV25]

[Vu07] [WA18]

[Wey12] [Xu17]

Experimental Design and Data Analysis Considerations”. In: Psychological Methods 27.5 (2022), pp. 874–894 (cit. on p. ii). Adith Swaminathan and Thorsten Joachims. “Batch Learning from Logged Bandit Feedback through Counterfactual Risk Minimization”. In: Journal of Machine Learning Research 16.52 (2015), pp. 1731–1755 (cit. on p. ii). G.W. Stewart and J. Sun. Matrix Perturbation Theory. Academic Press, 1990 (cit. on pp. iv, xx). Tobias Schnabel, Adith Swaminathan, Ashudeep Singh, Navin Chandak, and Thorsten Joachims. “Recommendations as Treatments: Debiasing Learning and Evaluation”. In: Proceedings of the 33rd International Conference on Machine Learning. Vol. 48. Proceedings of Machine Learning Research. PMLR, 2016, pp. 1670–1679 (cit. on p. ii). Joel A. Tropp. “User-Friendly Tail Bounds for Sums of Random Matrices”. In: Foundations of Computational Mathematics 12.4 (2012), pp. 389–434 (cit. on p. vii). Joel A. Tropp. “An Introduction to Matrix Concentration Inequalities”. In: Foundations and Trends in Machine Learning 8.1–2 (2015), pp. 1–230 (cit. on p. vii). T. Tao and V. Vu. “Random matrices: the circular law”. In: Commun. Contemp. Math. 10.02 (2008), pp. 261–307 (cit. on p. xxiv). Phuc Tran and Van Vu. “New perturbation bounds for low rank approximation of matrices: Beyond Eckart-Young-Mirsky”. In: (2025). arXiv:2511.08875 (cit. on pp. v, xii, xix, xxi). Linh Tran and Van Vu. “Fast exact recovery of noisy matrix from few entries: the infinity norm approach”. In: Proceedings of the 39th Conference on Neural Information Processing Systems (NeurIPS 2025). 2025 (cit. on pp. v, xii, xxi, xxvi). Phuc Tran and Van Vu. “Davis–Kahan Theorem under a moderate gap condition”. In: Communications in Contemporary Mathematics (2025). World Scientific, doi: 10.1142/S021919972550035X (cit. on pp. xii, xx, xxi). Phuc Tran and Nisheeth K. Vishnoi. “Perturbation Bounds for Low-Rank Inverse Approximations under Noise”. In: Proceedings of the 39th Conference on Neural Information Processing Systems (NeurIPS 2025). 2025 (cit. on p. xx). Phuc Tran and Van Vu. “New matrix perturbation bounds with relative norm: Perturbation of eigenspaces”. In: ArXiv preprint: 2409.20207 (2026). eprint: 2409.20207 (cit. on pp. v, xii, xx, xxi). Phuc Tran and Van Vu. “Eigenvalue Stability and New Perturbation Bounds for the extremal eigenvalues of a matrix”. In: (2026). arXiv:2603.19758 (cit. on pp. xxi, xxvi, xxvii). Phuc Tran, Nisheeth K. Vishnoi, and Van H. Vu. “Spectral Perturbation Bounds for Low-Rank Approximation with Applications to Privacy”. In: Proceedings of the 39th Conference on Neural Information Processing Systems (NeurIPS 2025). 2025 (cit. on pp. v, xii, xix–xxi). Van Vu. “Spectral norm of random matrices”. In: Combinatorica 27.6 (2007), pp. 721–736 (cit. on p. xxvi). Stefan Wager and Susan Athey. “Estimation and Inference of Heterogeneous Treatment Effects Using Random Forests”. In: Journal of the American Statistical Association 113.523 (2018), pp. 1228–1242 (cit. on p. ii). Hermann Weyl. “Das asymptotische Verteilungsgesetz der Eigenwerte linearer partieller Differentialgleichungen”. In: Mathematische Annalen 71.4 (1912), pp. 441–479 (cit. on p. xvii). Yiqing Xu. “Generalized Synthetic Control Method: Causal Inference with Interactive Fixed Effects Models”. In: Political Analysis 25.1 (2017), pp. 57–76 (cit. on p. iii).

A

Additional Discussion and Preliminaries

In this section, we collect remarks on regularity conditions and some additional preliminaries.

A.1

Additional Remarks on Informal Assumption 1.1 and Assumption 3.1

Remark A.1 (Thresholds). In the theorem below we take τa = 96KT ( a) for a ∈ {0, 1}. Any threshold of the same order gives the same bound after changing constants. The truncation by m−1 in the definition of pbi ( a) is only to avoid division by zero; under Assumption 3.1, it is inactive with high probability. Remark A.2 (Why the empirical rank is well-defined). The signal-to-noise condition ensures that the gap-selection step in Algorithm 1 is nonempty with high probability. Since rank( A( a)) ≤ r, we have σ1 ( a) = ∑rℓ=1 (σℓ ( a) − σℓ+1 ( a)), with σr+1 ( a) = 0. Hence σ1 ( a) > 120KrT ( a) implies that eub ( a) − A( a)∥op ≤ σs ( a) − σs+1 ( a) > 120KT ( a) for some s ≤ r. On the high-probability event ∥Y 12KT ( a), Weyl’s inequality gives e σsub ( a) − e σsub +1 ( a ) ≥ 96KT ( a ), so the algorithm selects at least s( a) is therefore an empirical effective rank, not one admissible truncation level. The selected b necessarily the algebraic rank of A( a).

A.2

Some Classical Results

In this section, we recall standard results used in Section 3, Section 4, and Section B.1. e := A + E. Then, Theorem A.3 (Weyl’s inequality [Wey12]). Let A, E be n × m matrices, and define A for any 1 ≤ i ≤ min{m, n},

|e λi − λi | ≤ ∥ E∥op

and

|e σi − σi | ≤ ∥ E∥op ,

e and σi , e where λi , e λi are the ith eigenvalues of A and A, σi are the corresponding singular values. Theorem A.4 (Cauchy’s integral theorem). Let Γ be a simple closed contour, and let f be an analytic function in the whole simply connected domain S containing Γ. Then  Z  f ( a), a inside Γ, 1 f (z) dz = . 0, 2πi Γ z − a a outside Γ In particular, for f (z) = z, one has 1 2πi

Z

  a,

z dz = 0, Γ z−a

a inside Γ, a outside Γ

.

e ∈ Rn×m , and let As , A es denote their Theorem A.5 (Eckart–Young–Mirsky bound [EY36]). Let A, A e − A. Then, respective best rank-s approximations. Set E := A
es − As ∥op ≤ 2 σs+1 + ∥ E∥op , ∥A where σs+1 is the (s + 1)th singular value of A.

B

Perturbation of Low-Rank Approximations in Row-Wise ℓ2 Norm

Let us recall the formal definition of row-wise ℓ2 norm. For each given pair of natural numbers k ≤ N, denote the standard basis of R N by {e N,k for 1 ≤ k ≤ N }. Given an n × m matrix M with n rows r1 , r2 , · · · , rn . Define ⊤ ∥ M∥2,∞ := max ∥ri ∥2 = max ∥en,i M ∥2 . 1≤ i ≤ n

1≤ i ≤ n

Back to our setting of perturbation of low-rank approximations. Let A be n × m matrix, admitting the singular decomposition: r

A = ∑ σi ui vi⊤ , i =1

where σ1 ≥ σ2 ≥ · · · ≥ σr are singular values with the corresponding pair of left/right singular vectors (ui , vi )1≤i≤r . For each 1 ≤ i ≤ r − 1, the ith singular gap is δi := σi − σi+1 . Let ER be an n × m random matrix of mean zero and E0 be an n × m deterministic matrix. We denote e = A + E. E = ER + E0 and A For a given rank-parameter s, our goal is to bound es − As ∥2,∞ . ∥A Before stating our main bound, we recall the following notions and definitions. Incoherence number. Define µ R :=

√

⊤ U ∥2 n max ∥en,i 1≤ i ≤ n

µC :=

and

√

m · max ∥Vem,j ∥2 . 1≤ j ≤ m

The the incoherence number of A is µ = max{µ R , µC }. Definition B.1. A random matrix ER is called (K, σ)-bounded if ER has entries satisfying E ER,ij = 0, E[| ER,ij |2 ] ≤ σ2 , E[| ER,ij |l ] ≤ K l −2 σ2 for all l ≥ 2 ∈ N and i ∈ [n], j ∈ [m]. We restate our key result - Theorem 4.2: Theorem 4.2. Let ER be a (K, σ)-bounded random matrix. There is a universal constant C > 0 such that the following holds. If δs ≥ 6(∥ ER ∥op + ∥ E0 ∥op ), then with probability at least 1 − O(1/(m+n)), es − As ∥2,∞ ∥A √ ≤ C r



K log4 (m+n) log (m + n) + √m+n 2



µ µ √ + √ m n

   √ σs · m + n δs σ + ∥ E0 ∥op .

(B.1)

To ease presentation, let us now focus on the case that n = Θ(m), in which the RHS simplifies to 

K log4 m log m + √ m 2



σs σ + ∥ E0 ∥op µ . δs

We obtain the following corollaries in this setting. The first one is when E0 = 0 (purely random noise), and the second one is when ∥ E0 ∥ is large.

Corollary B.2 (Random perturbation). Under the above setting, let ER be a (K, σ)-bounded random matrix. There is a universal constant C > 0 such that the following holds. If δs ≥ 4∥ ER ∥op , then with probability at least 1 − O(m−1 ),

√  K log4 m  σs ∥( A + ER )s − As ∥2,∞ ≤ C r log2 m + √ (B.2) µσ . δs m
e (1), then RHS becomes O e σs . Moreover, in many widely studied structured If r, K, σ, µ are all O δs models (e.g., spiked covariance, stochastic block, and graph Laplacian models), one typically e (1). In comparison to the existing bounds on has σs = O(δs ), yielding the clean bound O √ ∥( A + ER )s − As ∥ e.g., [EY36; MV25; TVV25; TV25a], which are at best O(∥ ER ∥) = O( m), our √ ∥ · ∥2,∞ -bound is smaller by a factor m. Since the optimal bound for ∥( A + ER )s − As ∥op is O(∥ ER ∥op ), and √ √ ∥( A + ER )s − As ∥op ≤ n∥( A + ER )s − As ∥2,∞ = Θ( m∥( A + ER )s − As ∥2,∞ ) , our bound is sharp up to a logarithmic factor. Corollary B.3. Under the above setting, let ER be a (K, σ)-bounded random matrix and ∥ E0 ∥ ≫ σ. There is a universal constant C > 0 such that the following holds. If δs ≥ 4(∥ ER ∥op + ∥ E0 ∥op ), then with probability at least 1 − O(m−1 ),

√  K log4 m  σs ∥( A + ER + E0 )s − As ∥2,∞ ≤ C r log2 m + √ µ ∥ E0 ∥op . δs m

(B.3)

Similar to the discussion after Corollary B.2, in many typical settings the right-hand side simplifies e (∥ E0 ∥op ), which is sharp up to logarithmic factors. For example, let m = n, ER = 0, A be to O diagonal, and E0 = cIm . Then

∥( A + ER + E0 )s − As ∥2,∞ = σ1 + c − σ1 = c = ∥ E0 ∥op . Unlike random noise ER , whose effect spreads across all directions, the deterministic component E0 can have a localized impact on the perturbation.

B.1

Proof of Theorem 4.2

In this subsection, we present the full detailed proof of Theorem 4.2. Recall the definition that for a given vector u ∈ R N , ∥u∥∞ := max1≤i≤ N |ui |. ! 0 A e as follows. Set A := Step 1 - Symmetrization. We symmetrize A and A , E := A⊤ 0 ! ! e 0 E 0 A , and then Ae = e⊤ . Given the singular decomposition of A = UΣV ⊤ , it is easy E⊤ 0 A 0 to see that A admits the following spectral decomposition

A=

U √ 2 V √ 2

U √ 2 − √V2

!

Σ 0 0 −Σ

!

⊤ U √ 2 ⊤ U √ 2

⊤ V √ 2 ⊤ − V√2

! ,

that is, A has the eigenvalues ±σi with the corresponding eigenvector u±i := √1 2

ui ± vi

! . Similarly,

e ±i . Ae has the eigenvalues ±e σi with the corresponding eigenvector u The best rank-2s approximation of A is

A2s :=

U √ 2 V √ 2

U √ 2 − √V2

!

Σs 0 0 −Σs

⊤ U √ 2 ⊤ U √ 2

!

⊤ V √ 2 ⊤ − V√2

! ,

which, in fact, equals !

A2s =

0 UΣs V ⊤ VΣs U ⊤ 0

!

Ae2s =

eΣ e⊤ e sV 0 U eΣ e⊤ e sU V 0

.

Similarly, ,

and then,

Ae2s − A2s =

eΣ e ⊤ − UΣs V ⊤ e sV 0 U ⊤ ⊤ eΣ e − VΣs U e sU V 0

! .

It yields es − As ) = e1⊤ (U eΣ e ⊤ − UΣs V ⊤ ) = e1⊤ (Ae2s − A2s ) , e sV e1⊤ ( A 2 2 2 and then e es − As ∥2,∞ := max ∥e⊤ ( A es − As )∥2 = max ∥e⊤ ∥A n,k m+n,k ( A2s − A2s ) 2 . 1≤ k ≤ n

1≤ k ≤ n

Step 2 - Contour representation of perturbation. Let Γ be a contour in C that encloses ±σ1 , ±σ2 , . . . , ±σs and excludes ±σs+1 , ±σs+2 , . . . , ±σr . The well-known contour integral representation [Kat80; SS90; Hig08] gives us: 1 2πi

Z Γ

s

s

i =1

i =1

⊤ z(zI − A)−1 dz = ∑ σi u+i u⊤ +i + ∑ (− σi )u−i u−i = A2s .

e with the corresponding eigenvectors ue1 , ue2 , . . . , uen . Let e λ1 ≥ · · · ≥ e λn denote the eigenvalue of A The construction of Γ (presented later) and the gap assumption 4∥ E∥ < δs ensure that the eigenvalues e λi lie inside Γ if and only if |e λi | ∈ {±e σ1 , ±e σ2 , . . . , ±e σs }. Then, similarly, we have R s s 1 ⊤ ⊤ − 1 e e e +i u e +i + ∑i=1 (−e e −i u e −i = A2s . dz = ∑i=1 e σi u σi )u 2πi Γ z ( zI − A) Thus, we obtain the following contour identity for the perturbation: 1 Ae2s − A2s = 2πi Therefore, e1⊤



Z Γ

e −1 )−1 − (zI − A)−1 ] dz. z[(zI − A)

(B.4)

Z  1 e e −1 )−1 − (zI − A)−1 ] dz. A2s − A2s = z · e1⊤ [(zI − A) 2πi Γ

This contour method is widely used in functional perturbation analysis, e.g., [Kat80; Hig08;KX16;  ⊤ e OVW23; TV25d; TVV25; TV25c; TV26a]. Note that one can trivially bound ∥e A2s − A2s ∥2 by 1

∥Ae2s − A2s ∥op , and then apply the existing bounds (e.g., [EY36; TVV25; TV25a]) in perturbation theory. However, this approach only gives a suboptimal bound that is off by a factor of up √ to n + m. Thus, obtaining a sharp bound with respect to the ℓ1 -norm remains a formidable analytical challenge. Step 3 - Contour expansion. We adapt the contour expansion method, introduced in [TV25b; TVV25; TV26a]. Using the Sherman-Morrison-Woodbury formula M−1 − ( M + N )−1 = ( M + e = A + E, we obtain N )−1 N M−1 [HJ12] and the fact that A e −1 = (zI − A)−1 E (zI − A) e −1 , (zI − A)−1 − (zI − A) and hence e −1 = (zI − A)−1 E (zI − A)−1 + (zI − A)−1 E [(zI − A) e −1 − (zI − A)−1 ]. (zI − A)−1 − (zI − A) e)−1 − (zI − By repeatedly applying the resolvent formula L times, one can expand ze1⊤ [(zI − A A)−1 ] into L

e −1 − (zI − A)−1 ]. ∑ z · e1⊤ (zI − A)−1 [E (zI − A)−1 ]k + e1⊤ (zI − A)−1 [E (zI − A)−1 ] L−1 · [(zI − A)

k =1

Therefore, 



e1⊤ Ae2s − A2s =

L

∑ Fk

!

+ Ftail , where

k =1

1 z · e1⊤ (zI − A)−1 [E (zI − A)−1 ]k dz and 2πi Γ Z 1 e −1 − (zI − A)−1 ] dz. Ftail := e⊤ (zI − A)−1 [E (zI − A)−1 ] L−1 · [(zI − A) 2πi Γ 1 Fk :=

Z

(B.5)

Thus, by the triangle inequality, we have 

e1⊤ Ae2s − A2s

L

 2

≤

∑ ∥ Fk ∥2

!

+ ∥ Ftail ∥2 .

k =1

We set L = 10ℓ log(m + n) (the constant ℓ will be chosen later). The remaining duty is to bound ∥ Fk ∥2 for each 1 ≤ k ≤ L and ∥ Ftail ∥2 . Indeed, we are going to show that the main part is ∑kL=1 ∥ Fk ∥2 and ∥ Ftail ∥2 is negligible. Inspired by the construction of Γ in [TVV25; TV25c; TV26b], we set Γ as as the union of two disjoint rectangles Γ+ ∪ Γ− , where: • Γ+ isolates {σ1 , σ2 , . . . , σs }, with (1) its left vertical edge intersecting the real axis at a0 := σs − δs /2; (2) its right vertical edge intersecting the real axis at a1 := 2σ1 ; (3) its height (from the real line) is T = 2σ1 . • Γ− isolates {−σ1 , −σ2 , . . . , −σs }, with (1) its left vertical edge intersecting the real axis at b1 := −2σ1 ; (2) its right vertical edge intersecting the real axis at b0 := −σs + δs /2; (3) its height (from the real line) is T = 2σ1 .

See the figure below for an illustration. Γ−

Γ+

−σs + δ2s −σ1

−σs

σs − δ2s

−σs+1

0

σs+1

σs

σ1

Step 4 - Bounding ∥ Fk ∥ for a natural number k ≤ L. We split the spectral decomposition of u u⊤

u u⊤

I

(zI − A)−1 = ∑ri=1 z+−i σ+i i + ∑ri=1 z−+i σ−i i + m+n

⊤ −(∑ri=1 u+i u⊤ + i +u− i u− i ) into P + Q, where z

:=

⊤ r u u⊤ u+ i u⊤ ∑ri=1 u+i u⊤ −i −i + i + u− i u− i +i + − ∑ z − σi ∑ z + σi z i =1 i =1

=

∑ (z − σi )z u+i u⊤+i + ∑ (z + σi )z u−i u⊤−i ,

r

P

r

r

σi

i =1

−σi

(B.6)

i =1

Im+n . and Q = z We can rewrite 2πiFk as what follows Z Γ

z · e1⊤ P[E (zI − A)−1 ]k dz +

Z Γ

z · e1⊤ Q[E (zI − A)−1 ]k dz

1 z · e1⊤ E (zI − A)−1 · [E (zI − A)−1 ]k−1 dz z Γ Γ Z Z Z 1 1 = z · e1⊤ P[E (zI − A)−1 ]k dz + z · e1⊤ E P · [E (zI − A)−1 ]k−1 dz + z · 2 e1⊤ E · [E (zI − A)−1 ]k−1 dz z z Γ Γ Γ = ··· " Z # Z k 1 ⊤ l 1 = ∑ z · l e1 E P[E (zI − A)−1 ]k−l dz + z · k+1 e1⊤ Ek dz. z z Γ l =0 Γ

=

Z

z · e1⊤ P[E (zI − A)−1 ]k dz +

Z

(B.7) Note that z = 0 is outside of the contour Γ. Thus, the last term k

1 2πi l =0

Fk := ∑

Z Γ

z·

R

Γ

z · zk1+1 e1⊤ Ek dz is zero, and hence

1 ⊤ l e E P[E (zI − A)−1 ]k−l dz. zl 1

By the triangle inequality, we have k

∥ Fk ∥2 ≤ ∑ Gk,l

where

l =0

1 Gk,l := · 2π

Z Γ

1 ⊤ l e E P[E (zI − A)−1 ]k−l 2 | dz|. l z −1 1

If l ≥ 1, splitting out the factor [E (zI − A)−1 ]k−l , we have Z k−l Z maxz∈Γ ∥E (zI − A)−1 ∥op 1 ⊤ l 1 1 ⊤ l · e E P | dz | ≤ · e E P 2 | dz|. 1 2 l − 1 k − l + 1 l 2π 2 π Γ z −1 1 Γ z The last inequality is true by our gap assumption that 4∥E ∥op < δs and the construction of Γ, which implies minz∈Γ,i∈[n] |z − λi | ≥ δs /2 ≥ 2∥ E∥, and hence

Gk,l ≤

max ∥E (zI − A)−1 ∥op ≤ max ∥(zI − A)−1 ∥op · ∥E ∥op = z∈Γ

z∈Γ

∥E ∥op ∥E ∥op 1 ≤ = . minz∈Γ,i∈[n] |z ± σi | 2∥E ∥op 2

For the second factor, we use the following lemma. Its proof is delayed to the next section.

Lemma B.4. Given a natural number l ≥ 1. Under the assumption of Theorem 4.2, with probability at least 1 − O((m + n)−2 ), there is a universal constant C > 0 such that 1

√

σs · (σ m + n ) l − 1 δs 3



R

Γ

1 e⊤ E l P z l −1 1 √

2

| dz|

2πC r

is at most

!  4 ∥E ∥ log ( m + n ) K log ( m + n ) op 0 √ √ log2 (m + n) + max ∥u±i ∥∞ + . 1≤ i ≤r m+n σ m+n

Thus, by Lemma (B.4), we further obtain that Gk,l , l ≥ 1 is at most

√ √ C r 2
l −1 σs · (σ m + n ) k 3 δs 2

!  4 ∥E ∥ log ( m + n ) K log ( m + n ) op 0 √ √ , max ∥u±i ∥∞ + log2 (m + n) + 1≤ i ≤r m+n σ m+n (B.8) with probability at least 1 − O((m + n)−2 ) for some universal constant C. 

Next, for l = 0, splitting out the factor [E (zI − A)−1 ]k−1 , we have k −1 Z maxz∈Γ ∥E (zI − A)−1 ∥op ze1⊤ PE (zI − A)−1 2 | dz| · 2π Γ Z 1 ze1⊤ PE (zI − A)−1 2 | dz|. ≤ k −1 · 2 π Γ

Gk,0 ≤

Similarly, for the second factor, we use the following technical lemma, whose proof will be presented in the next section. Lemma B.5. Under the assumption of Theorem 4.2, with probability at least 1 − O((m + n)−2 ), there is a universal constant C > 0 such that Z Γ

√ √ σs ze1⊤ PE (zI − A)−1 2 | dz| ≤ C r log2 (m + n) · (∥E0 ∥op + σ)( m + n ) · max ∥u±i ∥∞ . δs 1≤i≤r

By Lemma B.5, with probability at least 1 − O((m + n)−2 ), Z Γ

√ √ σs ze1⊤ PE (zI − A)−1 2 | dz| ≤ C r log2 (m + n) · (∥E0 ∥op + σ)( m + n ) · max ∥u±i ∥∞ . δs 1≤i≤r

Thus,

√ √ C r σs Gk,0 ≤ k log2 (m + n) · (∥E0 ∥op + σ) · ( m + n ) · max ∥u±i ∥∞ . δs 1≤i≤r 2

(B.9)

The estimates on Gk,l and the fact that max1≤i≤r ∥u±i ∥∞ ≥ √m1+n , imply that with probability at least 1 − O(k (m + n)−2 ), there is a universal constant C > 0, such that k

∥ Fk ∥ ≤ ∑ Gk,l l =0

√   √ C r σs K log4 (m + n) 2 √ ≤ k · (∥E0 ∥op + σ)( m + n ) log (m + n) + max ∥u±i ∥∞ . δs 1≤ i ≤r 2 m+n

(B.10)

Step 5 - Bounding ∥ Ftail ∥. We have 1 e −1 − (zI − A)−1 ] dz e1⊤ (zI − A)−1 [E (zI − A)−1 ] L−1 · [(zI − A) · 2π Γ 2   10σ1 −1 −1 L −1 − 1 e ∥op + ∥(zI − A)−1 ∥op ≤ · max ∥(zI − A) [E (zI − A) ] ∥op · ∥(zI − A) 2π z∈Γ   10σ1 L −1 L e −1 ∥op + ∥(zI − A)−1 ∥op . · ∥E ∥op · max ∥(zI − A)−1 ∥op · ∥(zI − A) ≤ 2π z∈Γ (B.11)

∥ Ftail ∥2 =

Z

The first inequality is obtained by the fact that the length of the contour Γ is at most 10σ1 . Moreover, maxz∈Γ ∥(zI − A)−1 ∥op = min |1z−λ | = δs1/2 = δ2s , the RHS is at most j

j

  2 10σ1 L −1 − 1 e ∥op + · (2∥E ∥op /δs ) · ∥(zI − A) . πδs δs

(B.12)

e −1 ∥op . For z = a + ib, we have Next, we upper bound ∥(zI − A) e = σmin (zI − A)

min

1≤ j ≤ m + n

q

b2 + ( a − e λ j )2 ,

here e λ j are eigenvalues of Ae. Let ℓ be the natural number such that σ1 ≤ nℓ . We have the following cases. e ≥ (m + n)−ℓ , and hence ∥(zI − A) e −1 ∥op ≤ (m + n)ℓ . • Case 1: b ≥ (m + n)−ℓ . Thus, σmin (zI − A) • Case 2: b < nℓ . So z must be on the vertical sides of Γ, i.e. Rez is either ±2σ1 or ±(σs − δs /2) If Rez = 2σ1 , then e ≥ σ1 − ∥E ∥op ≥ σ1 /2, and hence ∥(zI − A) e −1 ∥op ≤ 2 . σmin (zI − A) σ1 If Rez = σs − δs /2, then

∥zI − A − E0 ∥op ≤ n−ℓ + δs /2 + ∥E0 ∥op ≤ δs ≤ σ1 . e with respect to the random By applying [TV08; JSS22] on the pair of (zI − A − E0 , zI − A) noise E R , with probability at least 1 − (m + n)−C2 1

e ≥ (m + n)−2(C2 +2)ℓ+ 2 +o(1) . σmin (zI − A) It means, with probability at least 1 − (m + n)−2 , 1

e −1 ∥op ≤ (m + n)8ℓ+ 2 +o(1) . ∥(zI − A) All cases imply that with probability at least 1 − (m + n)−2 , e −1 ∥op ≤ (m + n)8ℓ+ 12 +o(1) . ∥(zI − A)

e −1 ∥, we obtain that with probabilCombining (B.11) with (B.12) and the upper bound of ∥(zI − A) ity at least 1 − (m + n)−2 ,     ∥E ∥op · σs ∥E ∥op · σs 1 1 (m + n)9ℓ+1+o(1) (m + n)9ℓ+1+o(1) √ ∥ Ftail ∥2 ≤ ≤ . ·√ · · · L − 1 δs δs (δs /2∥E ∥op ) 2 L −1 m+n m+n Thus, by setting L = 10ℓ · log(m + n), we have   ∥E ∥op · σs 1 ·√ , ∥ Ftail ∥2 ≤ o δs m+n which is negligible in comparison to ∑kL=1 ∥ Fk ∥2 . We finally obtain that e1⊤



Ae2s − A2s



  K log4 (m + n) σs 2 √ ≤ C r · (∥E0 ∥op + σ)( m + n ) log (m + n) + max ∥u±i ∥∞ . 2 δs 1≤ i ≤r m+n

√

√

This completes our proof. Remark B.6. We can slightly improve the upper bound on ∥ Fk ∥2 by combining exactly the upper √ bounds of Gk,0 , Gk,l without replacing √m1+n by max1≤i≤r ∥u±i ∥∞ . Indeed, ∥ Fk ∥2 /(C r ) is at most log2 (m + n) +

K log4 (m+n) √ m+n

2k

√

  ∥E0 ∥op σs , · ( m + n ) (X0 + σ) · max ∥u±i ∥∞ + √ δs 1≤ i ≤r m+n

where X0 := max1≤i≤r,1≤ j≤m+n |u⊤ ±i E0 v j | with v j , 1 ≤ j ≤ m + n eigenvectors of A. Define X0 := ⊤ e − A∥2,∞ is at most, max1≤i,j≤r |ui E0 v j |. As a result, we obtain that ∥ A     √ ∥ E0 ∥op K log4 (m + n) σs 1
1 2 √ log (m + n) + · ( m + n ) µ √ + √ ( X0 + σ ) + √ , δs m n m+n m+n which improves the RHS of Theorem 4.2.

C

Proofs of the Technical Lemmas

C.1

Bounding Gk,l , l ≥ 1 (Proof of Lemma B.4)

  −σi r ⊤ , which can be rearWe expand zl1−1 e1⊤ E l P into zl1−1 e1⊤ E l ∑ri=1 (z−σσi )z u+i u⊤ + u u ∑ − i i = 1 +i −i (z+σi )z i ranged as h i h i σ (e⊤ E l u )

r ∑ri=1 i zl (1z−σi+) i · u⊤ + i + ∑ i =1

−σi (e1⊤ E l u−i ) · u⊤ −i zl (z+σi )

.

Since {u±i }1≤i≤r are orthonormal system, ∥ zl1−1 e1⊤ E l P∥2 equals s

r 1 · |z|l

∑

σi2 ·(e1⊤ E l u+i )2 r i =1 |z−σi |2

+∑

σi2 ·(e1⊤ E l u−i )2 r i =1 |z+σi |2

σ2

σ2

∑ri=1 |z−σi |2 +∑ri=1 |z+σi |2

≤ max1≤i≤r |e1⊤ E l u±i | ·

i

i

|z|l

(C.1)

And hence, s 1 2π ·

R

Γ

1 ⊤ l e E P 2 | dz| ≤ z l −1 1

max1≤i≤r |e1⊤ E l u±i | 2π

·

R

σ2

σ2

∑ri=1 |z−σi |2 +∑ri=1 |z+σi |2 i

Γ

i

|z|l

| dz|.

(C.2)

The second factor can be bounded as follows. r σ2

σ2

∑ri=1 |z−iσi |2 + ∑ri=1 |z+iσi |2

√

1 σi · max | dz| l i ∈[r ] | z ± σi | Γ Γ |z| √ Z √ Z σi σi 1 1 · max · max = 2r | dz | + 2r | dz|. l l + − | z ± σ | | z ± σi | | z | | z | i ∈[ r ] i ∈[ r ] Γ Γ i

Z

|z|l

| dz| ≤

2r

Z

(C.3)

Moreover,

√

√ 1 σs σi 2r | dz| ≤ 2r · · max l δs /2 i ∈[r ] | z ± σi | Γ+ | z | Z

Z

√ 2σs 1 1 · | dz| ≤ 2r · 3 log l δs (σs /2)l −1 Γ+ | z |



10σ1 σs

 . (C.4)

i s Here, the first follows the fact that |z−σiσ | ≤ |Rezσ− , while the last inequality is obtained ≤ δsσ/2 σi | i √ R by replacing δD in by [TV26b]Section 8.2 by σs /2. A similar upper bound for 2r Γ− |z1|l · maxi∈[r] |z±σiσ | | dz|. i

Next, we handle the first factor as follows. We split E = E R + E0 , where E R is (K, σ)-bounded random matrix. By the triangle inequality, we have

|e1⊤ E l ui | = |e1⊤ (E R + E0 )l ui | l −1   l ⊤ l t l −t ≤ | e 1 E R ui | + ∑ ∥E R ∥op ∥E0 ∥op . t t =0

(C.5)

For the first term, by [TV25b][Lemma 4.2], we have !  3/2 3 log ( m + n ) K log ( m + n ) √ √ max ∥u±i ∥∞ + max |e1⊤ E Rl u±i | ≤ (2σ m + n)l · log(m + n) + , 1≤ i ≤r 1≤ i ≤r m+n m+n (C.6) with probability at least 1 − O((m + n)−2 ). For other terms, using [Vu07; BV16], with probability at least 1 − O((m + n)−2 ), we have   √ l l l −t t l −t . (C.7) ∥E R ∥op ∥E0 ∥op ≤ · (2σ m + n)t ∥E0 ∥op t t

√



R

Combing (C.2), (C.6), (C.7), (C.4), and the line below it, we obtain that

Γ

1 e⊤ E l P z l −1 1 √

2πC r

2

| dz|

is at most

√   √ 4σ m + n
l −1 σs K log4 (m + n) 2 √ max ∥u±i ∥∞ · (σ m + n ) log (m + n) + σs δs 1≤ i ≤r m+n √   l −1   2∥E0 ∥op l −t−1 √ l 4σ m + n
t σs ∥E0 ∥ log(m + n) √ +∑ · ·( m+n )· σs σs δs m+n t =0 t !   √ ∥E0 ∥op log(m + n) 1 σs K log4 (m + n) 2 √ √ ≤ l −1 · ( σ m + n ) log (m + n) + max ∥u±i ∥∞ + , δs 1≤ i ≤r 3 m+n σ m+n with probability at least 1 − O((m + n)−2 ) for some universal constant C > 0. The last inequality  √  is obtained by the fact that σs ≥ δs ≥ 6(∥E R ∥op + ∥E0 ∥op ) ≥ 3 4σ m + n + 2∥E0 ∥op . This proves Lemma B.4.

Bounding Gk,0 (Proof of Lemma B.5)

C.2

Consider the spectral decomposition of A := ∑in=1 λ j v j v⊤ j , where λ1 ≥ λ2 ≥ · · · ≥ λm+n are its eigenvalues with the corresponding eigenvector v j , 1 ≤ j ≤ m + n. It yields another way of v j v⊤

j +n −1 as follows. ⊤ presenting (zI − A)−1 = ∑m j=1 z−˘i . Thus, we rewrite ze1 P E ( zI − A) ! ! r r m + n v j v⊤ σ u − σ u j i +i1 i −i1 ⊤ ze1⊤ PE (zI − A)−1 = ∑ u⊤ u E ∑ +i + ∑ ( z − σ ) ( z + σi ) −i z − ˘i i i =1 i =1 j =1 ! m+n r r v⊤ σi u+i1 −σi u−i1 ⊤ j ⊤ = ∑ ∑ (u+ i E v j ) + ∑ (u− i E v j ) · . ( z − σ ) ( z + σ ) z − ˘j i i j =1 i =1 i =1

(C.8)

Therefore, ∥ze1⊤ PE (zI − A)−1 ∥2 equals v u 2 r um+n r σi u+i1 −σi u−i1 ⊤ 1 t ⊤ ∑ ∑ (z − σi ) (u+i E vj ) + ∑ (z + σi ) (u−i E vj ) · |z − ˘j |2 , j =1 i =1 i =1 which, by Cauchy-Schwartz inequality, is at most v ! um+n 2 2 r r u σ u − σ u 1 i + i1 i − i1 ⊤ ⊤ t ∑ 2r ∑ (u+ i E v j ) + ∑ (u− i E v j ) · (z − σi ) (z + σi ) | z − ˘ j |2 j =1 i =1 i =1

≤

max

1≤i ≤r,1≤ j≤m+n

v um+n u  r ⊤ |u±i E v j | · max ∥u±i ∥∞ t ∑ 2r ∑ 1≤ i ≤r

j =1

r σi2 σi2  1 + . ∑ 2 2 |z − σi | |z + σi | |z − λ j |2 i =1 i =1

s i ≤ δsσ/2 , ∥ze1⊤ PE (zI − A)−1 ∥2 is at most Since |z−σiσ | ≤ |Rezσ− σ| i

i

√

4rσs m+n . max | u⊤ ±i E v j | · max ∥u±i ∥∞ · δs min j |z − λ j | 1≤ i ≤r 1≤i ≤r,1≤ j≤m+n Thus, Z Γ

√

ze1⊤ PE (zI − A)−1 2 | dz| ≤

max

1≤ i ≤r 1≤ j ≤ m + n

4rσs m + n | u⊤ ±i E v j | · max ∥u±i ∥∞ · δs 1≤ i ≤r

Z

1 | dz|. Γ min j | z − λ j |

By [TV26b] Section 8.2, the RHS is at most

√

4rσs m + n max | u⊤ · 4 log ±i E v j | · max ∥u±i ∥∞ · δs 1≤i ≤r,1≤ j≤m+n 1≤ i ≤r



10σ1 δs

 .

On the other hand, by the Bernstein inequality, with probability at least 1 − O((m + n)−2 ), max

1≤i ≤r,1≤ j≤m+n

| u⊤ ±i E v j | ≤

max

1≤i ≤r,1≤ j≤m+n

⊤ | u⊤ ±i E R v j | + |u±i E0 v j | = O ( σ log n + ∥E0 ∥op ).

Combining these estimates above, with probability 1 − O((m + n)−2 ), there is a universal constant C, such that Z √ √ σs ze1⊤ PE (zI − A)−1 2 | dz| ≤ C r log2 (m + n)(∥E0 ∥op + σ) · ( m + n ) · max ∥u±i ∥∞ . δs 1≤i≤r Γ This proves our lemma.

C.3

Bounding ∥ ER ∥op (Proof of Lemma 4.1)

eub − E A eub + E eub . Our goal is to bound ∥ ER ∥op . By the triangle Recall the definition that ER := A inequality, eub − E A eub ∥op + ∥ E eub − E A eub ∥op + ∥ E eub ∥op = ∥ A eub ∥op . ∥ ER ∥op ≤ ∥ A Moreover, using [BV16], we obtain r   p p ( 1 − p ) ij ij ub ub e −EA e ∥op ≤ 6 · max{ MR , MC } + maxi,j aij · log(m + n) , where ∥A p2i r r pij (1− pij ) p (1− p ) m 2 and MC := max1≤ j≤m ∑in=1 a2ij ij p2 ij . MR := max1≤i≤n ∑ j=1 aij p2 i

Similarly, we also have r max 1≤i≤n

1≤ j ≤ m

(C.9)

(C.10)

i

∥ Eeub ∥op is at most 6

2 pik (1− pik ) , ∑m k =1 E ε ik p2

r

i

plj (1− plj ) ∑nl=1 E ε2lj p2

r



+ maxi,j KE

l

pij (1− pij ) p log(m + n). p2i

By the definition of K A that ∥ A∥∞ = maxi,j | aij | ≤ K A , we have r MR ≤ K A · max1≤i≤n

∑m j =1

pij (1− pij ) ≤ K A · max1≤i≤n p2i

q

pij

∑m j =1 p2 = K A ·

q

i

m min1≤i≤n pi .

Similarly, r MC ≤ K A · max1≤ j≤m

pij (1− pij ) ≤ KA · ∑in=1 p2i

r

pij

max1≤ j≤m ∑in=1 p min1≤i≤n pi

i

.

Therefore, max{ MR , MC } ≤ K A · Since

r maxi,j aij

pij (1− pij ) p2i

≤ K A · maxi,j

q

q

m +r p n min pi .

pij /pi ≤ KA · pi

we obtain eub − E A eub ∥op ≤ 12K A ∥A

p

log(m + n) ·

q

q

rp min pi ,

m +r p n min pi .

(C.11)

By a similar argument, we also obtain

∥ Eeub ∥op ≤ 12KE

p

log(m + n) ·

q

m +r p n min pi .

Combining the inequalities (C.11), (C.12) with (C.9), we finally obtain   q p m +r p n ∥ ER ∥op ≤ 12(K A + KE ) log(m + n) · min pi . We complete the proof of Lemma 4.1.

(C.12)

(C.13)