John Ellipsoids via Lazy Updates

Authors: david p woodruff
taisuke  yasuda
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 20:08:34

Abstract: We give a faster algorithm for computing an approximate John ellipsoid around $n$ points in $d$ dimensions. The best known prior algorithms are based on repeatedly computing the leverage scores of the points and reweighting them by these scores (Cohen et al., 2019). We show that this algorithm can be substantially sped up by delaying the computation of high accuracy leverage scores by using sampling, and then later computing multiple batches of high accuracy leverage scores via fast rectangular matrix multiplication. We also give low-space streaming algorithms for John ellipsoids using similar ideas.

Main paper: 
 introduction: 
 The John ellipsoid problem is a classic algorithmic problem, which takes as input a set of 𝑛 points {a 1 , a 2 , . . . , a 𝑛 } in 𝑑 dimensions, and asks for the minimum volume enclosing ellipsoid (MVEE) of these 𝑛 points. If 𝑃 is the convex hull of these 𝑛 points, then a famous result of Fritz John [Joh48] states that such an ellipsoid satisfies 1 𝑑 𝑄 ⊆ 𝑃 ⊆ 𝑄, and if 𝑃 is furthermore symmetric, then 1 √ 𝑑 𝑄 ⊆ 𝑃 ⊆ 𝑄. Equivalently, we may consider the 𝑛 input points to be constraints of a polytope {x : ⟨a 𝑖 , x⟩ ≤ 1, 𝑖 ∈ [𝑛]}, in which case the problem is to compute a maximal volume inscribed ellipsoid (MVIE). These two problems are related by taking polars, which corresponds to inverting the quadratic form defining the ellipsoid. In this work, we focus on the symmetric case, so that the polytope 𝑃 may be written as 𝑃 = {x : ‖Ax‖ ∞ ≤ 1}, where A denotes the 𝑛 × 𝑑 matrix with the 𝑛 input points a 𝑖 in the rows, and our goal is to output an ellipsoid 𝑄 which approximately satisfies 𝑄 ⊆ 𝑃 ⊆ √ 𝑑 • 𝑄.The John ellipsoid problem has far-reaching applications in numerous fields of computer science.In statistics, the John ellipsoid problem is equivalent to the dual of the D-optimal experiment design problem [Atw69, Sil13] , in which one seeks weights for selecting a subset of a dataset to observe in an experiment. Other statistical applications of John ellipsoids include outlier detection [ST80] and pattern recognition [Ros65, Gli98] . In the optimization literature, computing John ellipsoids is a fundamental ingredient for the ellipsoid method for linear programming [Kha79] , cutting plane methods [TKE88] , and convex programming [Sho77, Vai96] . Other applications in theoretical computer science include sampling [Vem05, CDWY18, GN23] , bandits [BCK12, HK16] , differential privacy [NTZ13] , coresets [TMF20, TWZ + 22], and randomized numerical linear algebra [Cla05, DDH + 09, WY23, BMV23] .We refer to a survey of Todd [Tod16] for an account on algorithms and applications of John ellipsoids.Following a long line of work on algorithms for computing John ellipsoids [Wol70, Atw73, KT93, NN94, Kha96, STT78, KY05, SF04, TY07, AST08, Yu11, HFR20] based on convex optimization techniques, the work of [CCLY19] developed a new approach towards computing approximate John ellipsoids via a simple fixed point iteration approach. The approach of [CCLY19] starts with an observation on the optimality condition for the dual problem, given by minimize𝑛 ∑︁ 𝑖=1 w 𝑖 -log det(A ⊤ WA) -𝑑 subject to w 𝑖 ≥ 0, 𝑖 ∈ [𝑛]where W = diag(w). 1 The optimality condition requires that the dual weights w satisfyEQUATIONfor each 𝑖 ∈ [𝑛], where 𝜏 𝑖 (A) for a matrix A denotes the 𝑖-th leverage score of A.Definition 1.1 (Leverage score). Let A ∈ R 𝑛×𝑑 . Then, for each 𝑖 ∈ [𝑛], we define the 𝑖-th leverage score as𝜏 𝑖 (A) := a ⊤ 𝑖 (A ⊤ A) -a 𝑖 = sup Ax̸ =0 [Ax](𝑖) 2 ‖Ax‖ 2 2 .The optimality condition of (1) can be viewed as a fixed point condition, which suggests the following iterative algorithm for computing approximate John ellipsoidsEQUATIONwhere W (𝑡) := diag(w (𝑡) ). After repeating this update for 𝑇 = 𝑂(𝜀 -1 log(𝑛/𝑑)) iterations starting with w (0) = 𝑑/𝑛 • 1 𝑛 , it can be shown that the ellipsoid 𝑄 defined by the quadratic form Q = A ⊤ W (𝑇 ) A, i.e. 𝑄 = {x ∈ R 𝑑 : x ⊤ Qx ≤ 1}, satisfiesEQUATIONNote that the computation of leverage scores can be done in 𝑂(𝑛𝑑 𝜔-1 ) time, where 𝜔 ≤ 2.371552 is the current exponent of fast matrix multiplication [DWZ23, WXXZ24] . Thus, this gives an algorithm running in time 𝑂(𝜀 -1 𝑛𝑑 𝜔-1 log(𝑛/𝑑)) for outputting an ellipsoid with the guarantee of (3).It is known that input matrices with additional structure admit even faster implementation of this algorithm. For instance, [CCLY19, SYYZ22] give faster algorithms for sparse matrices A for which the number of nonzero entries nnz(A) is much less than 𝑛𝑑. The work of [SYYZ22] shows that this approach can also be sped up for matrices A with small treewidth.Our first main result is a substantially faster algorithm for computing John ellipsoids. In the typical regime where 𝑛 ≫ 𝑑, our algorithm runs in 𝑂(𝜀 -1 𝑛𝑑) log(𝑛/𝑑) time to output an ellipsoid 𝑄 satisfying (3), and 𝑂(𝜀 -1 𝑛𝑑 2 ) log(𝑛/𝑑) time to output an ellipsoid 𝑄 which approximates the maximal volume up to a (1 + 𝜀) factor. We will discuss our techniques for this result in Sections 1.2.1 and 1.2.2. EQUATIONWe start by showing how to approximate leverage scores up to (1 + 𝜀) factors in Õ(𝑛𝑑) time, which had not been observed before to the best of our knowledge. Note that if we compute exact leverage scores using fast matrix multiplication, then this takes time 𝑂(𝑛𝑑 𝜔-1 ). Alternatively, sketching-based algorithms for approximate leverage scores are known, which gives the following running time for sparse matrices A with nnz(A) nonzero entries. Theorem 1.2 ([SS11, DMMW12, CW13]). There is an algorithm which, with probability at least 1 -𝛿, outputs 𝜏 ′ 𝑖 for 𝑖 ∈ [𝑛] such that 𝜏 ′ 𝑖 = (1 ± 𝜀)𝜏 𝑖 (A) and runs in time 𝑂(𝜀 -2 nnz(A) log(𝑛/𝛿)) + poly(𝑑𝜀 -1 log(𝑛/𝛿)).If the goal is to compute (1 + 𝜀)-approximate leverage scores for a dense 𝑛 × 𝑑 matrix, then we are not aware of a previous result which does this in a nearly linear Õ(𝑛𝑑) time, which we now show: Theorem 1.3. There is an algorithm which, with probability at least 1 -𝛿, outputs𝜏 ′ 𝑖 for 𝑖 ∈ [𝑛] such that 𝜏 ′ 𝑖 = (1 ± 𝜀)𝜏 𝑖 (A) in time 𝑂(𝑛𝑑) poly log(𝜀 -1 log(𝑛/𝛿)) + 𝑂(𝑛) poly(𝜀 -1 log(𝑛/𝛿)).Our improvement comes from improving the running time analysis of a sketching-based algorithm of [DMMW12] by using fast rectangular matrix multiplication. We will need the following result on fast matrix multiplication: Theorem 1.4 ([Cop82, Wil11, Wil24]). There is a constant 𝛼 ≥ 0.1 and an algorithm for multiplying a 𝑚 × 𝑚 and a 𝑚 × 𝑚 𝛼 matrix in 𝑂(𝑚 2 poly log 𝑚) time, under the assumption that field operations can be carried out in 𝑂(1) time.By applying the above result in blocks, we get the following version of this result for rectangular matrix multiplication. Corollary 1.5. There is a constant 𝛼 ≥ 0.1 and an algorithm for multiplying a 𝑛 × 𝑑 for 𝑛 ≥ 𝑑 and a 𝑑 × 𝑡 matrix in 𝑂(𝑛𝑑 + 𝑛𝑡 1/𝛼+1 ) poly log 𝑡 time, under the assumption that field operations can be carried out in 𝑂(1) time.Proof. Let 𝑚 = 𝑡 1/𝛼 . If 𝑑 ≤ 𝑚, then matrix multiplication takes only 𝑂(𝑛𝑑𝑡) = 𝑂(𝑛𝑡 1/𝛼+1 ) time, so assume that 𝑑 ≥ 𝑚. We partition the first matrix into an 𝑂(𝑛/𝑚) × 𝑂(𝑑/𝑚) block matrix with blocks of size 𝑚 × 𝑚 and the second into an 𝑂(𝑑/𝑚) × 1 block matrix with blocks of size 𝑚 × 𝑚 𝛼 . By Theorem 1.4, each block matrix multiplication requires 𝑂(𝑚 2 poly log 𝑚) time, and we have 𝑂(𝑛𝑑/𝑚 2 ) of these to do, which gives the desired running time.That is, the above result shows that when multiplying an 𝑛 × 𝑑 matrix A with a 𝑑 × 𝑡 matrix for a much smaller 𝑡, then this multiplication can be done in roughly 𝑂(𝑛𝑑) poly log 𝑡 time. The work of [DMMW12] shows that the leverage scores of A can be written as the row norms of AR for a 𝑑 × 𝑡 matrix with 𝑡 = 𝑂(𝜀 -2 log(𝑛/𝛿)), and thus this gives us the result of Theorem 1.3.By using Theorem 1.3, we already obtain a John ellipsoid algorithm which runs in time 𝑂(𝜀 -1 𝑛𝑑) log(𝑛/𝑑) poly log(𝜀 -1 log 𝑛) time, which substantially improves upon prior algorithms for dense input matrices A. We now show how to obtain further improvements by using the idea of lazy updates. At the heart of our idea is to only compute the quadratic forms for the John ellipsoids for most iterations, and defer the computation of the weights until we have computed roughly 𝑂(log 𝑛) iterations. At the end of this group of iterations, we can then compute the John ellipsoid weights via fast matrix multiplication as used in Theorem 1.3, which allows us to remove the suboptimal poly log log 𝑛 terms in the dominating term of the running time. Theorem 1.6. Given A ∈ R 𝑛×𝑑 , let 𝑃 be the polytope defined by 𝑃 = {x ∈ R 𝑑 : ‖Ax‖ ∞ ≤ 1}. For 𝜀 ∈ (0, 1), there is an algorithm, Algorithm 3, that runs in time 𝑂(𝜀 -1 𝑛𝑑)(log(𝑛/𝑑) + poly log(𝜀 -1 log 𝑛)) + 𝑂(𝑛) poly(𝜀 -1 log 𝑛) + 𝑂(𝑛 0.1 )𝑑 𝜔+1 𝜀 -3 (log 𝑛) 2 and returns an ellipsoid 𝑄 such that 1√ 1+𝜀 • 𝑄 ⊆ 𝑃 ⊆ √ 𝑑 • 𝑄 with probability at least 1 -1/ poly(𝑛).The full proof of this result is given in Section 2.Let Q (𝑡) = A ⊤ W (𝑡) A, where the weights w (𝑡) are defined as (2). Note that with this notation, the update rule for the iterative algorithm of [CCLY19] can be written asEQUATIONThus, given high-accuracy spectral estimates to the quadratics Q (𝑡) , we can recover the weights w (𝑡) 𝑖 to high accuracy in 𝑂(𝑑 𝜔 ) time per iteration by evaluating the quadratic forms a ⊤ 𝑖 (Q (𝑡) ) -a 𝑖 and then multiplying them together. This approach is useful for fast algorithms if we only need to to do this for a small number of indices 𝑖 ∈ [𝑛]. This is indeed the case if we only need these weights for a row sample of √ W (𝑡) A, which is sufficient for computing a spectral approximation to the next quadratic form Q (𝑡) . Furthermore, we only need low-accuracy leverage scores (up to a factor of, say, 𝑛 0.1 ) to obtain a good row sample, which can be done quickly for all 𝑛 rows [LMP13, CLM + 15]. Thus, by repeatedly sampling rows of √ W (𝑡) A, computing high-accuracy weights on the sampled rows, and then building an approximate quadratic, we can iteratively compute high-accuracy approximations Q(𝑡) to the quadratics Q (𝑡) . More formally, our algorithm takes the following steps:• We first compute low-accuracy leverage scores of √ W (𝑡-1) A, which can be done in 𝑂(𝑛𝑑) time. This gives us the weights w (𝑡) to low accuracy, say u (𝑡) , for all 𝑛 rows.• We use the low-accuracy weights u (𝑡) to obtain a weighted subset of rows of √ U (𝑡) A which spectrally approximates Q (𝑡) . Note, however, that we do not yet have the sampled rows of √ W (𝑡) A to high accuracy, since we do not know the weights w (𝑡) to high accuracy.• If we only need the weights w (𝑡) for a small number of sampled rows 𝑆 ⊆ [𝑛], then we can explicitly compute these using (4), since we inductively have access to high-accuracy quadratics Q(𝑡 ′ ) for 𝑡 ′ = 0, 1, 2, . . . , 𝑡 -1. These can then be used to build Q(𝑡) .While this algorithm allows us to quickly compute high-accuracy approximate quadratics Q(𝑡) , this algorithm cannot be run for too many iterations, as the error in the low-accuracy leverage scores u (𝑡) grows to poly(𝑛) factors in 𝑂(log 𝑛) rounds. This is a problem, as this error factor directly influences the number of leverage score samples needed to approximate Q(𝑡) . We will now use the fast matrix multiplication trick from the previous Section 1.2.1 to fix this problem. Indeed, after 𝑂(log 𝑛) iterations, we will now have approximate quadratics Q(1) , Q(2) , . . . , Q(𝑡) for 𝑡 = 𝑂(log 𝑛). Now we just need to compute the 𝑛 John ellipsoid weights which are given byv (𝑡) 𝑖 = 𝑡 ∏︁ 𝑡 ′ =1 ‖e ⊤ 𝑖 A( Q(𝑡 ′ -1) ) -1/2 ‖ 2 2 .To approximate this quickly, we can approximate each term ‖e⊤ 𝑖 A( Q(𝑡 ′ -1) ) -1/2 ‖ 2 2 by ‖e ⊤ 𝑖 A( Q(𝑡 ′ -1) ) -1/2 G (𝑡 ′ ) ‖ 2 2 for a random Gaussian matrix G (𝑡 ′ ), by the Johnson-Lindenstrauss lemma [JL84] . Here, the number of columns of the Gaussian matrix can be taken to be poly(𝜀 -1 log 𝑛), so now all we need to compute is the matrix productA • [( Q(0) ) -1/2 G (0) , ( Q(1) ) -1/2 G (1) , . . . , ( Q(𝑡) ) -1/2 G (𝑡)] which is the product of a 𝑛 × 𝑑 matrix and a 𝑑 × 𝑚 matrix for 𝑚 = poly(𝜀 -1 log 𝑛). By Theorem 1.4, this can be computed in 𝑂(𝑛𝑑 poly log 𝑚) time. However, this resetting procedure is only run 𝑂(𝜀 -1 ) times across the 𝑇 = 𝑂(𝜀 -1 log 𝑛) iterations, so the running time contribution from the resetting is just 𝑂(𝜀 -1 𝑛𝑑) poly log(𝜀 -1 log 𝑛) + 𝑂(𝑛) poly(𝜀 -1 log 𝑛). Overall, the total running time of our algorithm is 𝑂(𝜀 -1 𝑛𝑑)(log(𝑛/𝑑) + poly log(𝜀 -1 log 𝑛)) + 𝑂(𝑛) poly(𝜀 -1 log 𝑛). Remark 1.7. In general, our techniques can be be viewed as a way of exploiting the increased efficiency of matrix multiplication when performed on a larger instance by delaying large matrix multiplication operations, so that the running time is 𝑂(𝜀 -1 𝑛𝑑 log(𝑛/𝑑)) + 𝑂(𝜀 -1 )𝑇 𝑟 where 𝑇 𝑟 is the time that it takes to multiply A by a 𝑑 × 𝑟 matrix for 𝑟 = poly(𝜀 -1 log 𝑛). While we have instantiated this general theme by obtaining faster running times via fast matrix multiplication, one can expect similar improvements by other ways of exploiting the economies of scale of matrix multiplication, such as parallelization. For instance, we recover the same running time if we can multiply 𝑟 vectors in parallel so that 𝑇 𝑟 = 𝑂(𝑛𝑑).The problem of computing John ellipsoids is also well-studied in the streaming model, where the input points a 𝑖 arrive one at a time in a stream [MSS10, AS15, WY22, MMO22, MMO23]. The streaming model is often considered when the number of points 𝑛 is so large that we cannot fit all of the points in memory at once, and the focus is primarily on designing algorithms with low space complexity. Our second result of this work is that approaches similar to the one we take to prove Theorem 1.6 in fact also give a low-space implementation of the iterative John ellipsoid algorithm of [CCLY19] . Theorem 1.8 (Streaming algorithms). Given A ∈ R 𝑛×𝑑 , let 𝑃 be the polytope defined by 𝑃 = {x ∈ R 𝑑 : ‖Ax‖ ∞ ≤ 1}. Furthermore, suppose that A is presented in a stream where the rows a 𝑖 ∈ R 𝑑 arrive one by one in a stream. For 𝜀 ∈ (0, 1), there is an algorithm, Algorithm 1, that makes 𝑇 = 𝑂(𝜀 -1 log(𝑛/𝑑)) passes over the stream, takes 𝑂(𝑑2 𝑇 ) time to update after each new row, and returns an ellipsoid 𝑄 such that 1√ 1+𝜀 • 𝑄 ⊆ 𝑃 ⊆ √ 𝑑 • 𝑄.Furthermore, the algorithm uses at most 𝑂(𝑑 2 𝑇 ) words of space.In Section 1.2.2, we showed that by storing only the quadratics Q(𝑡) and only computing the weights∏︀ 𝑡 𝑡 ′ =1 a 𝑖 ( Q(𝑡 ′ -1)) -a 𝑖 as needed, we could design fast algorithms for John ellipsoids. In fact, this idea is also useful in the streaming setting, since storing all of the weights w (𝑡) 𝑖 requires 𝑂(𝑛) space per iteration, whereas storing the quadratics Q(𝑡) requires only 𝑂(𝑑 2 ) space per iteration. Furthermore, in the streaming setting, we may optimize the update running time by instead storing the pseudo-inverse of the quadratics ( Q(𝑡) ) -and then updating them by using the Sherman-Morrison low rank update formula. 2 This gives the result of Theorem 1.8. for 𝑡 = 0 to 𝑇 do3: Let Q (𝑡) = 0 4: for 𝑖 = 1 to 𝑛 do 5: if 𝑡 = 0 then 6: Let Q (𝑡) ← Q (𝑡) + a 𝑖 a ⊤ 𝑖 7: else 8: Let w (𝑡) 𝑖 = ∏︀ 𝑡 𝑡 ′ =1 a ⊤ 𝑖 (Q (𝑡 ′ -1) ) -a 𝑖 9: Let Q (𝑡) ← Q (𝑡) + w (𝑡)𝑖 a 𝑖 a ⊤ 𝑖 10:return 1 𝑇 +1 ∑︀ 𝑇 𝑡=0 Q (𝑡)Throughout this paper, A will denote an 𝑛×𝑑 matrix whose 𝑛 rows are given by vectors a 𝑖 ∈ R 𝑑 . For positive numbers 𝑎, 𝑏 > 0, we write 𝑎 = (1±𝜀)𝑏 to mean that(1-𝜀)𝑏 ≤ 𝑎 ≤ (1+𝜀)𝑏. For symmetric positive semidefinite matrices Q, R, we write Q = (1±𝜀)R to mean that (1-𝜀)R ⪯ Q ⪯ (1+𝜀)R,where ⪯ denotes the Löwner order on PSD matrices.2 Fast algorithms 
approximating the quadratics: 
 We will analyze the following algorithm, Algorithm 2, for approximating the quadratics Q (𝑡) of the iterative John ellipsoid algorithm.Fix a row 𝑖. Note then that for each iteration 𝑡, ‖e ⊤ 𝑖 A( Q(𝑡-1) ) -1/2 G‖ 2 2 is distributed as an independent 𝜒 2 variable with 𝑘 degrees of freedom, scaled by ‖e ⊤ 𝑖 A( Q(𝑡-1) ) Let Q(0) = A ⊤ W(0) A for W(0) = diag( w(0) ).for 𝑡 = 1 to 𝑇 do 4:Compute A( Q(𝑡-1) ) -1/2 G for a 𝑑 × 𝑘 Gaussian matrix G.Let u(𝑡) 𝑖 = u (𝑡-1) 𝑖 • ‖e ⊤ 𝑖 A( Q(𝑡-1) ) -1/2 G‖ 2 2 for each 𝑖 ∈ [𝑛]. 6:Let S (𝑡) be a (1 + 𝜀)-approximate leverage score sample of √ U (𝑡) A (Thm. 2.3).For rows 𝑖 sampled by S (𝑡) , set v(𝑡) 𝑖 = ∏︀ 𝑡 𝑡 ′ =1 a ⊤ 𝑖 ( Q(𝑡 ′ -1) ) -a 𝑖 . 8: Let Q(𝑡) = (S (𝑡) √ V (𝑡) A) ⊤ S (𝑡) √ V (𝑡) A. 9: return { Q(𝑡) } 𝑇 𝑡=0 after 𝑇 iterations, u (𝑇 ) 𝑖 = 𝑇 ∏︁ 𝑡=1 ‖e ⊤ 𝑖 A( Q(𝑡-1) ) -1/2 G‖ 2 2 ,which is distributed asEQUATIONfor i.i.d. 𝜒 2 variables 𝑋 𝑡 with 𝑘 degrees of freedom. We will now bound each of the terms in this product.We need the following bound on products of 𝜒 2 variables, which generalizes [SSK17 E 𝑒 𝑠 log 𝑋𝑖 = E 𝑋 𝑠 𝑖 = 1 2 𝑘/2 Γ(𝑘/2) ∫︁ ∞ 0 𝑥 𝑠 𝑥 𝑘/2-1 𝑒 -𝑥/2 𝑑𝑥 = 2 𝑠 Γ(𝑠 + 𝑘/2) Γ(𝑘/2) = 𝐶 𝑠,𝑘 .Then by Chernoff bounds,Pr {︃ 𝑡 ∏︁ 𝑖=1 𝑋 𝑖 ≤ 1 𝑅 }︃ = Pr {︃ exp (︃ 𝑡 ∑︁ 𝑖=1 -𝑠 log 𝑋 𝑖 )︃ ≥ 𝑅 𝑠 }︃ ≤ E[𝑒 -𝑠 log 𝑋𝑖 ] 𝑡 𝑅 -𝑠 = 𝐶 𝑡 -𝑠,𝑘 𝑅 -𝑠 and Pr {︃ 𝑡 ∏︁ 𝑖=1 𝑋 𝑖 ≥ 𝑅 }︃ = Pr {︃ exp (︃ 𝑡 ∑︁ 𝑖=1 𝑠 log 𝑋 𝑖 )︃ ≥ 𝑅 𝑠 }︃ ≤ E[𝑒 𝑠 log 𝑋𝑖 ] 𝑡 𝑅 -𝑠 = 𝐶 𝑡 𝑠,𝑘 𝑅 -𝑠Using the above result, we can show that if the 𝜒 2 variables have 𝑘 = 𝑂(1/𝜃) degrees of freedom and the number of rounds 𝑇 is 𝑐 log 𝑛 for a small enough constant 𝑐, then the product of the 𝜒 2 variables ∏︀ 𝑇 𝑡=1 𝑋 𝑡 will be within a 𝑛 𝜃 factor for some small constant 𝜃 > 0. Lemma 2.2. Fix a small constant 𝜃 > 0 and let 𝑘 = 𝑂(1/𝜃). Let 𝑇 = 𝑐 log 𝑛 for a sufficiently small constant 𝑐 > 0. Let 𝑋 1 , 𝑋 2 , . . . , 𝑋 𝑇 be 𝑡 i.i.d. 𝜒 2 variables with 𝑘 degrees of freedom. Then, It follows that with probability at least 1 -1/ poly(𝑛), the products of 𝜒 2 variables appearing in (5) are bounded by 𝑛 𝜃 for every row 𝑖 ∈ [𝑛] and for every iteration 𝑡 ∈ [𝑇 ]. We will condition on this event in the following discussion.Pr {︃ 1 𝑛 𝜃 ≤ 𝑡 ∏︁ 𝑖=1 𝑋 𝑖 ≤ 𝑛 𝜃 }︃ ≥ 1 - 1 poly(𝑛) . Proof. Set 𝑠 = 𝑘/2.In the previous section, we have established that the products of 𝜒 2 variables in (5) are bounded by 𝑛 𝜃 . We will now use this fact to show that the quadratics Q(𝑡-1) in Algorithm 2 are good approximate John ellipsoid quadratics. We will use the following leverage score sampling theorem for this.Theorem 2.3 ([DMM06, SS11]). Let A ∈ R 𝑛×𝑑 . Let 𝜏 ′ 𝑖 ≥ 𝜏 𝑖 (A) and let 𝑝 𝑖 = min{1, 𝜏 ′ 𝑖 /𝛼} for 𝛼 = Θ(𝜀 2 )/ log(𝑑/𝛿). If S is a diagonal matrix with S 𝑖,𝑖 = 1/√ 𝑝 𝑖 with probability 𝑝 𝑖 and 0 otherwise for 𝑖 ∈ [𝑛], then with probability at least 1 -𝛿,for all x ∈ R 𝑑 , ‖SAx‖ 2 2 = (1 ± 𝜀)‖Ax‖ 2 2 .This theorem, combined with the bounds on 𝜒 2 products, gives the following guarantee for the approximate quadratics Q(𝑡) .Lemma 2.4. Fix a small constant 𝜃 > 0 and let 𝑘 = 𝑂(1/𝜃). Suppose that the leverage score sample in Line 6 oversamples by a factor of 𝑂(𝑛 2𝜃 ), that is, uses leverage score estimates 𝜏 ′ 𝑖 such that 𝜏 ′ 𝑖 ≥ 𝑂(𝑛 2𝜃 )𝜏 𝑖 ( √ U (𝑡) A). Then, with probability at least 1 -1/ poly(𝑛), we have for every𝑡 ∈ [𝑇 ] that Q(𝑡) = (1 ± 𝜀)A ⊤ V (𝑡) A (6)where v(𝑡) 𝑖 = ∏︀ 𝑡-1 𝑡 ′ =1 a ⊤ 𝑖 ( Q(𝑡 ′ ) ) -a 𝑖 for 𝑖 ∈ [𝑛]. Furthermore, Q(𝑡) can be computed in 𝑂(𝑛 2𝜃 )𝑇 𝜀 -2 𝑑 𝜔+1 log 𝑛 time in each iteration.Proof. We will first condition on the success of the event of Lemma 2.2, so that the 𝜒 2 products in (5) are bounded by 𝑛 𝜃 factors for every row 𝑖 ∈ [𝑛] and every iteration 𝑡 ∈ [𝑡]. We will also condition on the success of the leverage score sampling for all 𝑇 iterations. Note that by (5) and the bound the 𝜒 2 products, u(𝑡) 𝑖 is within a 𝑂(𝑛 2𝜃 ) factor of v (𝑡)𝑖 , and thus S (𝑡) is a correct leverage score sample for √ V (𝑡) A. We thus have thatQ(𝑡) = (S (𝑡) √︀ V (𝑡) A) ⊤ S (𝑡) √︀ V (𝑡) A = (1 ± 𝜀)A ⊤ V (𝑡) A.For the running time, note that S (𝑡) samples at most 𝑂(𝜀 -2 𝑛 2𝜃 𝑑 log 𝑛) rows in each iteration, and each sampled row 𝑖 requires 𝑂(𝑑 𝜔 𝑇 ) to compute v Let 𝐵 = 𝑂(𝑐 -1 𝜀 -1 ) and 𝑇 = 𝑂(𝑐 log(𝑛/𝑑)) for a sufficiently small constant 𝑐. for 𝑏 = 1 to 𝐵 do 5:Let { Q(𝑡) } 𝑇 𝑡=0 be given by APPROXQUADRATIC(A, w(0) ) (Algorithm 2).Let G (𝑡-1) for 𝑡 ∈ [𝑇 ] be a random 𝑑 × 𝑚 Gaussian for 𝑚 = 𝑂(𝜀 -2 (𝐵𝑇 ) 2 log 𝑛).7:Compute A • [( Q(0) ) -1/2 G (0) , ( Q(1) ) -1/2 G (1) , . . . , ( Q(𝑇 ) ) -1/2 G (𝑇 ) ]. 8: Let w(𝑏,𝑡) 𝑖 = ∏︀ 𝑡 𝑡 ′ =1 1 𝑚 ‖e ⊤ 𝑖 A( Q(𝑡 ′ -1) ) -1/2 G (𝑡 ′ -1) ‖ 2 2 for each 𝑖 ∈ [𝑛]. 9:Let w(0) = w(𝑏,𝑇 ) .10:Let w = 1 𝐵𝑇 ∑︀ 𝑏,𝑡w(𝑏,𝑡) .11:return A ⊤ WAWe now combine the subroutine for approximating quadratics from Section 2.1 with a resetting procedure using fast matrix multiplication to obtain our full algorithm for quickly computing John ellipsoids. The full algorithm is presented in Algorithm 3.We will need the following theorem, which summarizes results of [CCLY19] on guarantees of the fixed point iteration algorithm (2) under approximate leverage score computations. Theorem 2.5 (Lemma 2.3 and Lemma C.4 of [CCLY19] ). Let A ∈ R 𝑛×𝑑 and let 𝑃 = {x : ‖Ax‖ ∞ ≤ 1}. Let 𝑇 = 𝑂(𝜀 -1 log(𝑛/𝑑)). Suppose thatw (𝑡) 𝑖 = (1 ± 𝜀)𝜏 𝑖 ( √︀ W (𝑡-1) A) for all 𝑡 ∈ [𝑇 ] and 𝑖 ∈ [𝑛].Then, if 𝑄 is the ellipsoid given by the quadratic form A ⊤ W (𝑇 ) A, then1 √ 1+𝜀 • 𝑄 ⊆ 𝑃 ⊆ √ 𝑑 • 𝑄.We also use the Johnson-Lindenstrauss lemma, which is a standard tool for randomized numerical linear algebra, especially in the context of approximating leverage scores [SS11, DMMW12, LMP13 We will now give a proof of Theorem 1.6.Proof of Theorem 1.6. We first argue the correctness of this algorithm. By Theorem 2.5, it suffices to verify that w(𝑏,𝑡) 𝑖 computes (1 + 𝜀)-approximate leverage scores in order to prove the claimed guarantees of Theorem 1.6. For the updates within APPROXQUADRATIC (Algorithm 2), we have by Lemma 2.4 thatQ(𝑡) = (1 ± 𝜀)A ⊤ V (𝑡) A.Thus,v (𝑡-1) 𝑖 • ‖e ⊤ 𝑖 A( Q(𝑡-1) ) -1/2 ‖ 2 2 = (1 ± 𝜀)v (𝑡-1) 𝑖 • ‖e ⊤ 𝑖 A(A ⊤ V (𝑡-1) A) -1/2 ‖ 2 2 = (1 ± 𝜀)𝜏 𝑖 ( √︀ V (𝑡-1) A)and thus the weights v(𝑡)𝑖 output by APPROXQUADRATIC (Algorithm 2) indeed satisfy the requirements of Theorem 2.5.For the weights computed in Line 8 of Algorithm 3, note that we also compute these weights, but this time with approximation error from the application of the Gaussian matrix 1 𝑚 G (𝑡) to speed up the computation. Applying Theorem 2.6 with 𝜀 set to 𝜀/𝐵𝑇 and failure probability 𝛿 = 1/ poly(𝑛), we have that1 𝑚 ‖e ⊤ 𝑖 A( Q(𝑡 ′ -1) ) -1/2 G (𝑡 ′ -1) ‖ 2 2 = (1 ± 𝜀/𝐵𝑇 )‖e ⊤ 𝑖 A( Q(𝑡 ′ -1) ) -1/2 ‖ 2 2 .Note then that Line 8 takes a product of at most 𝐵𝑇 of these approximations, so the total error in the approximation is at most(1 ± 𝜀/𝐵𝑇 ) 𝐵𝑇 = (1 ± 𝑂(𝜀)).The running time is given by 𝐵𝑇 iterations of the inner loop of APPROXQUADRATIC and 𝐵 iterations of the fast matrix multiplication procedure in Line 7 of Algorithm 3. The inner loop of APPROXQUADRATIC requires 𝑂(𝑛𝑑) time to compute the product with the 𝑑 × 𝑘 Gaussian as well as the time to compute the approximate quadratic, which is bounded in Lemma 2.4. Altogether, this gives the claimed running time bound.Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes] Justification: We discuss limitations of our work in Section 3, and pose open questions towards closing these gaps. Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
future directions: 
 In this work, we developed fast algorithms and low-space streaming algorithms for the problem of computing John ellipsoids. Our fast algorithms use a combination of using lazy updates together with fast matrix multiplication to substantially improve the running time of John ellipsoids, and we apply similar ideas to obtain a low-space streaming implementation of the John ellipsoid algorithm.Our results have several limitations that we discuss here, which we leave for future work to resolve. First, our algorithm makes crucial use of fast matrix multiplication in order to get running time improvements. However, this makes it a less attractive option for practical implementations, and also makes the polynomial dependence on 𝜀 in the Õ(𝑛) poly(𝜀 -1 ) term rather large. Thus, it is an interesting question whether the running time that we obtain in Theorem 1.6 is possible without fast matrix multiplication. Question 3.1. Is there an algorithm with the guarantee of Theorem 1.6 that avoids fast matrix multiplication?More generally, it is an interesting question to design algorithms for approximating John ellipsoids with optimal running time. This is an old question which has been studied in a long line of work [Wol70, Atw73, KT93, NN94, Kha96, STT78, KY05, SF04, TY07, AST08, Yu11, HFR20], and we believe that the investigation of this question will lead to further interesting developments in algorithms research. Question 3.2. What is the optimal running time of approximating John ellipsoids?For instance, one interesting question is whether it is possible to obtain nearly linear time algorithms that run in time Õ(𝑛𝑑) + Õ(𝑛) poly(𝜀 -1 ), or even input sparsity time algorithms that run in time Õ(nnz(A))+ Õ(𝑛) poly(𝜀 -1 ). The resolution of such questions for the least squares linear regression problem has led to great progress in algorithms, and studying these questions for the John ellipsoid problem may have interesting consequences as well.John ellipsoids are closely related to ℓ 𝑝 Lewis weights, which give a natural ℓ 𝑝 generalization of John ellipsoids and leverage scores and have been a valuable tool in randomized numerical linear algebra. There has been a recent focus on fast algorithms for computing (1 + 𝜀)-approximate ℓ 𝑝 Lewis weights [FLPS22, AGS24] , and thus it is natural to ask whether developments in algorithms for John ellipsoids would carry over to algorithms for ℓ 𝑝 Lewis weights as well. Question 3.3. What is the optimal running time of approximating ℓ 𝑝 Lewis weights?Finally, we raise questions concerning streaming algorithms for approximating John ellipsoids. In Theorem 1.8, we gave a multi-pass streaming algorithm which obtains a (1+𝜀)-optimal John ellipsoid. A natural question is whether a similar algorithm can be achieved in fewer passes, or if there are pass lower bounds for computing (1 + 𝜀)-optimal John ellipsoids. Question 3.4. Are there small space streaming algorithms for approximating John ellipsoids up to a (1 + 𝜀) factor which make fewer than 𝑂(𝜀 -1 log(𝑛/𝑑)) passes, or do small space streaming algorithms necessarily require many passes?We note that there are one-pass small space streaming algorithms if we allow for approximation factors that scale as 𝑂( √ log 𝑛) rather than(1 + 𝜀) [WY22, MMO22, MMO23].1. Claims Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? Answer: [Yes] Justification: [TODO] Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes]Answer: [NA] Justification: [TODO] Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details. • While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark). Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: [TODO] Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [NA] Justification: [TODO] Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 13. New Assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA] Justification: [TODO] Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.For a weight vector w ∈ R 𝑛 , we will often write the corresponding 𝑛 × 𝑛 diagonal matrix diag(w) as the capitalized version W.We note that storing the pseudo-inverse may increase the space complexity by polynomial factors in the bit complexity model. 
 
 Appendix:  


Decision: This paper improves upon the computational state of the art for computing John's ellipsoids. The focus is on giving low-storage algorithms amenable to streaming data setting, and is based on a simple idea of delaying leverage scores via sampling, and then using fast matrix multiplication in batches to speed up the bottleneck with high approximation accuracy. The main weaknesses in my view are that the paper is presented in a way that feels narrow in scope, and as one reviewer points out, it can come across as a fairly incremental contribution. However, the results are correct and carefully executed, and the paper can be further improved incorporating some comments from useful reviews such as further acceleration using Woodbury approximations. I would encourage the authors to give more context on how these contributions can allow users to tackle problems that may have been out of reach due to the previous bottleneck, and just widen the scope of the paper even though the contribution is technical in nature. 

I recommend acceptance after studying the paper myself; due to the brevity and lack of substantive criticism in review hiPW, I chose to dismiss that score in the decision making process

Reviews: review: paper_summary: paper_summary: This paper considers the computing of an approximate John ellipsoid. They improve the algorithm by lazy update and fast matrix multiplication. They also give low-space streaming algorithms using similar ideas.
strength_weakness: strength_weakness: Major:
- This paper is badly written.
- The authors do not give justifications in the checklist.

Minor:
- Please check the capitalization in the references.
- can you explicitly give the space complexity in Theorem 1.8?
questions: questions: See weakness
limitations: limitations: NA
paper_summary: paper_summary: The authors give a leverage score approximation algorithm that runs in nearly linear time for dense matrices. Specifically, to find leverage scores for a matrix $A \in \mathbb{R}^{n \times d}$, they give an algorithm based on fast matrix multiplication that runs in time $\widetilde{O}(nd)$. This, combined with a "lazy update" trick, yield a new fast algorithm to approximate the John ellipsoid for a symmetric convex polyhedron represented as the set $\{x \in \mathbb{R}^d \colon \| Ax \|_{\infty} \le 1\}$. The same ideas also transfer to a $\sim \log n$-pass streaming setting to find the John ellipsoid, yielding a new low space complexity algorithm for finding the John ellipsoid that has a stronger approximation guarantee than existing one-pass solutions.
strength_weakness: strength_weakness: I wish more applications had been discussed. Leverage scores are a pretty fundamental primitive that get used in a lot of problems in optimization and numerical linear algebra, and it would have been nice to write down concrete runtime improvements (if any) that emerge from Theorem 1.3. I alluded to some in the previous section, and I would be happy to raise my score if the authors either confirm the above or present a few more settings in which concrete runtime improvements are realized. Alternately, I'd also raise the score if the authors discussed a bit about the barriers behind extending their approach to these settings.

Mild issues:

I think there are typos in Algorithm 1 and a suboptimal runtime guarantee. I think Line 8 should read:

$w_{i}^{(t)} = \prod_{t'=1}^{t} a_i^{\top}(A^{\top}W^{(t')}A)^{-1}a_i$. In particular, as written, Line 8 multiplies across the first $t$ rows for a fixed quadratic, whereas the loop counter should be running over the quadratics instead.

I also think Line 9 should read:

$Q^{(t)} \gets Q^{(t)} + w_i^{(t)}a_ia_i^{\top}$ (I think you need an outer product and not an inner product, same thing for Line 6)

I think you also need an averaging step at the very end, as the [CCLY19] algorithm returns the weights $w = 1/T \cdot \sum_{t'=1}^{T} w_{t'}$.

Finally, I think this algorithm actually can be implemented in $d^2T \le d^{\omega}T$ time per row. This is because of associativity:

$\prod_{t=1}^t a_i^{\top}Q_{t}^{-1}a_i = \prod_{t=1}^t a_i^{\top}(Q_{t}^{-1}a_i)$

Now, $(Q_{t}^{-1}a_i)$ is a vector that can be formed in time $d^2$ via naive multiplications, and then $a_i^{\top}(Q_{t}^{-1}a_i)$ is a dot product that can be found in $d$ time. We do this $t$ times for each term in the product and then multiply them all together for a total time of $d^2t$ for step $t$.
questions: questions: See all of the above.
limitations: limitations: Yes
paper_summary: paper_summary: The paper proposes provably faster algorithm for computing John ellipsoids. The first idea is to approximate Leverage score with an existing faster algorithm given by Theorem 1.3. The second idea is to use fast matrix multiplication and lazy updates to reduce polylog suboptimality.
strength_weakness: strength_weakness: - despite that the main result does not rely on Theorem 1.3, there is no proof or reference for it.
- contribution should clarified. The proof of Theorem 1.6 mainly relies on the existing results.
questions: questions: Could you please clarify on your particular contribution?

Could you please refer to the proof of Theorem 1.3?

I am not an expert in John Ellipsoids. For me it seemed that the contribution was incremental, because the theoretical result combines existing techniques. By the way, the new result is obtained.

I would also like you to comment on the most recent papers and preprints from 2023 and 2024 if any, to confirm your contribution.
limitations: limitations: NA

Comments: title: response to author response to reviewer response to author response to review
responding to author response to review
comment: OK, thanks a lot for checking the above! I have updated my confidence score accordingly.
Indeed, constant factor leverage score approximation takes time $\tilde O(\mathrm{nnz}(A) + d^\omega)$ time, as we cite in Theorem 1.2. We have informally referred to this as input sparsity time, as we assume a regime where $\mathrm{nnz}(A) \geq n \gg \mathrm{poly}(d)$. For dense matrices, our $\tilde O(nd)$ dominating running time would upper bound $\mathrm{nnz}(A)$. To re-iterate, our running time improvement is focused on $(1+\epsilon)$ approximation of leverage scores, where we improve from $\epsilon^{-2}\mathrm{nnz}(A) + \mathrm{poly}(d/\epsilon)$ to $\tilde O(nd) + \mathrm{poly}(d/\epsilon)$. However, this does not give substantial improvements when $\epsilon$ is constant.

Yes, we can use the Sherman--Morrison formula to improve the update time of maintaining the inverse of the quadratic in our streaming algorithm, thank you for pointing that out! It is unclear if the Sherman--Morrison would improve the running time of the offline fast algorithm.
Hi, thanks for the clarification! I have a couple of followup questions:

> if only constant factor approximations are required (as in the case of $\ell_p$ Lewis weight sampling), then leverage score approximation can be done in input sparsity time

I thought leverage scores for general inputs took $\widetilde{O}(nnz(A) + d^{\omega})$ time to compute (which looks like it could be worse than the runtime you get for dense matrices)? They only run in truly input sparsity time when $A$ has additional structure (e.g. if $A$ is a graph edge incidence matrix)? See, e.g., Lemmas 7-10 of https://arxiv.org/pdf/1408.5099. Of course, it's likely I misunderstood what you meant by input-sparsity time though.

> The [runtime] is due to computing the matrix inverse

Ah, I forgot to mention -- can you maintain $Q_{t}^{-1}$ using Sherman-Morrison (https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula)? In particular, each update to $Q$ looks like a rank-$1$ update, so the formula for updating $Q^{-1}$ follows from Sherman-Morrison (and it looks like to me that each update to $Q^{-1}$ runs in $d^2$ time). The formula as written in the wiki page only applies to invertible matrices but I think it holds (https://mathoverflow.net/questions/146831/sherman-morrison-type-formula-for-moore-penrose-pseudoinverse) for maintaining the pseudoinverse as well if the matrix you are keeping track of symmetric -- which it looks like it is.

Thanks again for following up!