Banded Square Root Matrix Factorization for Differentially Private Model Training

Authors: nikita  kalinin
christoph  lampert
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 13:58:29

Abstract: Current state-of-the-art methods for differentially private model training are based on matrix factorization techniques. However, these methods suffer from high computational overhead because they require numerically solving a demanding optimization problem to determine an approximately optimal factorization prior to the actual model training. In this work, we present a new matrix factorization approach, BSR, which overcomes this computational bottleneck. By exploiting properties of the standard matrix square root, BSR allows to efficiently handle also large-scale problems. For the key scenario of stochastic gradient descent with momentum and weight decay, we even derive analytical expressions for BSR that render the computational overhead negligible. We prove bounds on the approximation quality that hold both in the centralized and in the federated learning setting. Our numerical experiments demonstrate that models trained using BSR perform on par with the best existing methods, while completely avoiding their computational overhead.

Main paper: 
 introduction: 
 We study the problem of differentially private (DP) model training with stochastic gradient descent (SGD) in the setting of either federated or centralized learning. This task has recently emerged as one of the most promising ways to train powerful machine learning models but nevertheless guarantee the privacy of the used data, which led to a number of studies, both theoretical as well as application-driven [Abadi et al., 2016 , Yu et al., 2020 , Zhang et al., 2021 , Kairouz et al., 2021 , Denisov et al., 2022] . The state of the art in the field are approaches based on the matrix factorization (MF) mechanism [Li et al., 2015 , Henzinger et al., 2024] , which combines theoretical guarantees with practical applicability [Choquette-Choo et al., 2023a ,c,b, 2024] 1 It is based on the observation that all iterates of SGD are simply linear combinations of model gradients, which are computed at intermediate time steps. Consequently, the iterates can be written formally as the result of multiplying the matrix of coefficients, called workload matrix, with the row-stacked gradient vectors. To preserve the privacy of the training data in this process one adds suitably scaled Gaussian noise at intermediate steps of the computation. The MF mechanism provides a way to select the noise covariance structure based on a factorization of the workload matrix into two matrices.Identifying the minimal amount of noise necessary to achieve a desired privacy level requires solving an optimization problem over all possible factorizations, subject to data participation constraints. For some specific settings, the optimal solutions have been characterized: for streaming learning, when each data batch contributes at most once to the gradients, Li et al. [2015] presented a formulation of this problem as a semi-definite program. Henzinger et al. [2024] proved that a square root factorization of the workload matrix is asymptotically optimal for different linear workloads, including continual summation and decaying sums.In this work, our focus lies on the settings that are most relevant to machine learning tasks: workload matrices that reflect SGD-like optimization, and participation schemes in which each data batch can potentially contribute to more than one gradient vector, as it is the case for standard multi-epoch training. Assuming that this happens at most once every b steps, for some value b ≥ 1, leads to the problem of optimal matrix factorization in the context of b-min-separated participation sensitivity. Unfortunately, as shown in Choquette-Choo et al. [2023a] , finding the optimal matrix factorization in this setting is computationally intractable. Instead, the authors proposed an approximate solution by posing additional constraints on the solution set. The result is a semi-definite program that is tractable, but still has high computational cost, making it practical only for small to medium-sized problem settings.Subsequent work concentrated on improving or better understanding the factorizations for specific algorithms, such as plain SGD without gradient clipping, momentum, or weight decay [Koloskova et al., 2023] or specific, e.g. convex, objective functions [Choquette-Choo et al., 2024] . Often, streaming data was asssumed, i.e. each data item can contribute at most to one model update, which is easier to analyze theoretical, but further removed from real-world applications [Dvijotham et al., 2024] . A concurrent line of works has also focused on scaling matrix factorizations for large-scale training. McMahan et al. [2024] extended the Buffered Linear Toeplitz (BLT) mechanism to support multi-participation in federated learning, improving privacy-utility tradeoffs, while ensuring memory efficiency. McKenna [2024] improved DP Banded Matrix Factorization's scalability, enabling it to handle millions of iterations and large models efficiently. These advancements enhance the practicality of differentially private training in real-world applications.Our work aims at a general-purpose solution that covers as many realistic scenarios as possible. Our ultimate goal is to make general-purpose differentially private model training as simple and efficient to use as currently dominating non-private technique. Our main contributions are:1. We introduce a new factorization, the banded squared root (BSR), which is efficiently computable even for large workload matrices and agnostic to the underlying training objective. For matrices stemming from SGD optimization potentially with momentum and/or weight decay, we even provide closed form expressions. 2. We provide general lower bounds on the approximate error for any factorization, along with specific upper and lower bounds for the BSR, Square Root, and baseline factorizations, in the contexts of both single participation (streaming) and repeated participation (e.g., multi-epoch) training. 3. We demonstrate experimentally that BSR's approximation error is comparable to the state-of-the-art method, and that both methods also perform comparably in real-world training tasks.Overall, the proposed BSR factorization achieves training high-accuracy models with provable privacy guarantees while staying computationally efficient even for large-scale training tasks. 
background: 
 Our work falls into the areas of differentially private (stochastic) optimization, of which we remind the reader here, following mostly the description of Denisov et al. [2022] and Choquette-Choo et al. [2023a] . The goal is to estimate a sequence of parameter vectors, Θ = (θ 1 , . . . , θ n ) ∈ R d , where each θ i is a linear combination of update vectors, x 1 , . . . , x i ∈ R d , that were computed in previous steps, typically as gradients of a model with respect to training data that is meant to stay private. We assume that all x i have a bounded norm, x i ≤ ζ. Compactly, we write Θ = AX, where the lower triangular workload matrix A contains the coefficients and X ∈ R n×d is formed by stacking the update vectors as rows. With different choices of A, the setting then reflects many popular first-order optimization algorithms, in particular stochastic gradient descent (SGD), potentially with momentum and/or weight decay. Depending on how exactly x 1 , . . . , x n are obtained, the setting can express different centralized as well as federated training paradigms. To formalize this aspect, we adopt the concept of b-min-separated participation [Choquette-Choo et al., 2023a] . For some integer b ≥ 1, it states that if a data item (e. x i + ζ[C -1 Z] [i,•] θ i ← update(θ i-1 , xi ),// SGD model updates Output: Θ = (θ 1 , . . . , θ n ) federated learning) contributed to an update x i , the earliest it can contribute again is the update x i+b . Additionally, let 1 ≤ k ≤ n b be the maximal number any data point can contribute. In particular, this notion also allows us to treat in a unified way streaming data (b = n or k = 1), as well as unrestricted access patterns (k = n with b = 1), but also intermediate settings, such as multi-epoch training on a fixed-size dataset.The matrix factorization approach [Li et al., 2015] adopts a factorization A = BC of the workload matrix and computes Θ MF = B(CX + Z), where Z is Gaussian noise that is chosen appropriately to make the intermediate result CX + Z private to the desired level. Algorithm 1 shows the resulting algorithm in pseudocode. It exploits the fact that instead of explicit multiplication by C and B, standard optimization toolboxes can be employed with suitably modified update vectors, because also Θ MF = A(X + C -1 Z), and multiplication by A corresponds to performing the optimization.Different factorizations recover different algorithms from the literature. For example, B = A, C = Id recovers DP-SGD [Abadi et al., 2016] , where noise is added directly to the gradients. Conversely, B = Id, C = A simply adds noise to each iterate of the optimization [Dwork et al., 2006] . However, better choices than these baselines are possible, in the sense that they can guarantee the same levels of privacy with less added noise, and therefore potentially with higher retained accuracy. The reason lies in the fact that B and C play different roles: B acts as a post-processing operation of already private data. Hence, it has no further effect on privacy, but it influences to what amount the added noise affects the expected error in the approximation of Θ. Specifically, for Z ∼ N (0; s Id),EQUATIONIn contrast, CX is the quantity that is meant to be made private. Doing so requires noise of a strength proportional to C's sensitivity, sens(C) := sup X∼X CX -CX F , where the neighborhood relation, X ∼ X , indicates that the two sequences of update vectors differ only in those entries that correspond to a single data item2 As shown in Choquette-Choo et al. [2023a] , in the setting of b-min-separated repeated participation, it holds thatEQUATIONwhereΠ k,b = { π ⊂ {1, . . . , n} : |π| ≤ k ∧ ({i, j} ⊂ π ⇒ i = j ∨ |i -j| ≥ b) },is the set of possible b-min-separated index sets with at most k participation. Furthermore, (2) holds even with equality if all entries of C C are non-negative.Combining (1) with s = sens k,b (C) yields a quantitative measure for the quality of a factorization. Definition 1. For any factorization A = BC, its expected approximation error isEQUATIONwhere the 1/ √ n factor is meant to make the quantity comparable across different problem sizes.The optimal factorization by this reasoning would be the one of smallest expected approximation error. Unfortunately, minimizing (3) across all factorizations it is generally computationally intractable. Instead, Choquette-Choo et al. [2023a] propose an approximately optimal factorization. Definition 2. For a workload matrix A, let S be the solution to the optimization problem arg minS∈S n + trace[A A S -1 ] subject to diag(S) = 1 and S [i,j] = 0 for |i -j| ≥ b, (4)where S n + is the cone of positive definite n × n matrices. Then, A = BC is called the approximately optimal factorization (AOF), if C is lower triangular and fulfills C C = S.The optimization problem ( 4) is a semi-definite program (SDP), and can therefore be solved numerically using standard packages. However, this is computationally costly, and for large problems (e.g. n > 5000) computing the AOF solution is impractical. This poses a problem for real-world training tasks, where the number of update steps are commonly thousands or tens of thousands.Solving (4) itself only approximately can mitigate this problem to some extent, but as we will discuss in Section 4, this can lead to robustness problems, especially because the recovery of C from S in Definition 2, e.g. by a Cholesky decomposition, tends to be sensitive to numerical errors. 
banded square root factorization: 
 In the following section, we introduce our main contribution: a general-purpose factorization for the task of differentially private stochastic optimization that can be computed efficiently even for large problem sizes. Definition 3 (Banded Square Root Factorization). Let A ∈ R n×n be a lower triangular workload matrix with strictly positive diagonal entries. Then, we call A = C 2 the square root factorization (SR), when C denotes the unique matrix square root that also has strictly positive diagonal entries. Furthermore, for any bandwidth p ∈ {1, . . . , n}, we define the banded square root factorization of bandwidth p (p-BSR) of A, asA = B |p| C |p| (5)where C |p| is created from C by setting all entries below the p-th diagonal to 0,EQUATION)Note that determining the SR, and therefore any p-BSR, is generally efficient even for large workload matrices, because explicit recursive expressions exist for computing the square root of a lower triangular matrix [Björck and Hammarling, 1983, Deadman et al., 2012] .In the rest of this work, we focus on the case where the workload matrix stems from SGD with momentum and/or weight decay, and we show that then even closed form expressions for the entries of C |p| exist that renders the computational cost negligible.We recall the update steps of SGD with momentum and weight decay:θ i = αθ i-1 -ηm i for m i = βm i-1 + x i (7)where x 1 , . . . , x n are the update vectors, η > 0 is the learning rate and 0 ≤ β < 1 is the momentum strength and 0 < α ≤ 1 is the weight decay parameter. Note that our results also hold for β = 0, i.e. without momentum, and for α = 1, i.e., without weight decay. In line with real algorithms and to avoid degenerate cases, we assume β < α throughout this work. The update vectors are typically gradients of the model with respect to its parameters, but additional operations such as normalization or clipping might have been applied.Unrolling the recursion, we obtain an expression for θ i as a linear combination of update vectors asEQUATIONConsequently, the workload matrix has the explicit form A = ηA α,β forEQUATION)As one can see, A α,β is a lower triangular Toeplitz-matrix, so it is completely determined by the entries of its first column. In the following, we use the notation LDToep(m 1 , . . . , m n ) to denote a lower triangular Toeplitz matrix with first column m 1 , . . . , m n , i.e. A α,β = LDToep(a 0 , . . . , a n-1 ).Our first result is an explicit expression for the positive square root of A α,β (and thereby its p-BSR).Theorem 1 (Square-Root of SGD Workload Matrix). Let A α,β be the workload matrix (9). ThenA α,β = C 2 α,β for C α,β = LDToep(c 0 , . . . , c n-1 ), with c 0 = 1 and c j = j i=0 α j-i r j-i r i β i for j = 1, . . . , n -1 with coefficients r i = -1/2 i .For any p ∈ {1, . . . , n}, the p-banded BSR matrix C |p| α,β is obtained from this by setting all coefficients c j = 0 for j ≥ p.Proof sketch. The proof be found in Appendix F.1. Its main idea is to factorize A α,β into a product of two simpler lower triangular matrices, each of which has a closed-form square root. We show that the two roots commute and that the matrix C α,β is their product, which implies the theorem.We first establish that the p-BSR for SGD can be computed efficiently even for large problem sizes.Lemma 1 (Efficiency of BSR). The entries of C |p| α,β can be determined in runtime O(p log p), i.e., in particular independent of n.Proof sketch. As a lower triangular Toeplitz matrix, C |p| α,β is fully determined by the values of its first column. By construction c p+1 , . . . , c n = 0, so only the complexity of computing c 1 , . . . , c p-1 matters. These can be computed efficiently by writing them as the convolution of vectors (α i r i ) i=0,...,p-1 and (β i r i ) i=0,...,p-1 and, e.g., employing the fast Fourier transform.Note that for running Algorithm 1, the matrix B of the factorization is not actually required. However, one needs to know the sensitivity of C |p| α,β , as this determines the necessary amount of noise. The following theorem establishes that for a large class of matrices, including the BSR in the SGD setting, this is possible exactly and in closed form. Theorem 2 (Sensitivity for decreasing non-negative Toeplitz matrices). Let M = LDToep(m 0 , . . . , m n-1 ) be a lower triangular Toeplitz matrix with decreasing non-negative entries, i.e.EQUATION)where M [•,1+jb] denotes the (1 + jb)-th column of M .Proof sketch. The proof can be found in Appendix F.2. It builds on the identity (2), which holds with equality because of the non-negative entries of M . Using the fact that the entries of M are non-increasing one establishes that an optimal b-separated index set is {1,1 + b, • • • , 1 + (k -1)b}.From this, the identity (10) follows.Corollary 1. The sensitivity of the p-BSR for SGD can be computed using formula (10).Proof sketch. It suffices to show that the coefficients c 0 , . . . , c n-1 of Theorem 1 are monotonically decreasing. We do so by an explicit computation, see Appendix F.3.Having established the efficiency of BSR, we now demonstrate its suitability for high-quality model training. To avoid corner cases, we assume that n b is an integer, which does not affect the asymptotic behavior. We also discuss only the case in which the update vectors have bounded norm ζ = 1. Results for general ζ can readily be derived using the linearity of the sensitivity with respect to ζ.We first discuss the case of model training with single participation (k = 1), where more precise results are possible than the general case. Our main result are bounds on the expected approximation error of the square root factorization that, in particular, prove its asymptotic optimality. Theorem 3 (Expected approximation error with single participation). Let A α,β ∈ R n×n be the workload matrix (9) of SGD with momentum 0 ≤ β < 1 and weight decay parameter 0 < α ≤ 1, where α > β. Assume that each data item can contribute at most once to an update vector (e.g. single participation, k = 1). Then, the expected approximation error of the square root factorization,A α,β = C 2 α,β , fulfills 1 ≤ E(C α,β , C α,β ) ≤ 1 (α -β) 2 log 1 1 -α 2 (11)for α < 1, andEQUATION)Proof sketch. For the proof, we establish a relations between sens 1,n (C) and C α,β F , and then we bound the resulting expressions by an explicit analysis of the norm. For details, see Appendix F.5.The following two results provide context for the interpretation of Theorem 3. Theorem 4. Assume the setting of Theorem 3. Then, for any factorization A α,β = BC with C C ≥ 0, the expected approximation error fulfillsEQUATION)Proof sketch. The theorem is the special case k = 1 of Theorem 8, which we state in the next section and prove in Section F.9.Theorem 5. Assume the setting of Theorem 3. Then, the baseline factorizationsEQUATION)Proof sketch. The result follows from an explicit analysis of the coefficients, see Appendix F.6.Discussion. Theorems 3 to 5 provide a full characterization of the approximation quality of the square root factorization as well as its alternatives: 1) the square root factorization has asymptotically optimal approximation quality, because the upper bounds in Equation ( 12) match the lower bounds in Equation (13); 2) the AOF from Definition 2 also fulfills the conditions of Theorem 4. Therefore, it must also adhere to the lower bound (13) and cannot be asymptotically better than the square root factorization; 3) the approximation qualities of the baseline factorizations in Equation ( 14) and ( 15) are asymptotically worse than optimal in the α = 1 setting, and worse by a constant factor for α < 1.The BSR factorization can be applied even in more general scenarios, such as with varying learning rates. However, in this case, the workload matrix will no longer be Toeplitz. This makes it difficult to provide analytical guarantees for the matrix, but it can still be applied numerically.3.4 Approximation Quality -Repeated Participation.We now provide mostly asymptotic statements about the approximation quality of BSR and baselines in the setting where data items can contribute more than once to the update vectors.Theorem 6 (Approximation error of BSR). Let A α,β ∈ R n×n be the workload matrix (9) of SGD with momentum 0 ≤ β < 1 and weight decay 0< α ≤ 1, with α > β. Let A α,β = B |p| α,β C |p| α,β, be its banded square root factorization as in Definition 3. Then, for any b ∈ {1, . . . , n}, p ≤ b, and k ∈ {1, . . . , n b } it holds:EQUATION)Proof sketch. For the proof, we separately bound the sensitivity of C |p| α,β and the Frobenius norm of B |p| α,β . The former is straightforward because of the matrix's band structure. The latter requires an in-depth analysis of the inverse matrix' coefficient. Both steps are detailed in Appendix F.7.The following results provide context for the interpretation of Theorem 6.Theorem 7 (Approximation error of Square Root Factorization). Let A α,β ∈ R n×n be the workload matrix (9) of SGD with momentum 0 ≤ β < 1 and weight decay 0 < α ≤ 1, with α > β. Let A α,β = C 2 α,β be its square root factorization. Then, for any b ∈ {1, . . . , n} and k = n b it holds:EQUATIONProof sketch. We bound sens k,b (C α,β ) and C α,β F using the explicit entries for C α,β from Theorem 1. Details are provided in Appendix F.8.Theorem 8. Assume the setting of Theorem 6. Then, for any factorization A α,β = BC with C C ≥ 0, the approximation error fulfillsEQUATIONProof sketch. The proof is based on the observation that X F Y F ≥ XY * for any matrices X, Y , where Theorem 9. Assume the setting of Theorem 6. Then, the baseline factorizationsEQUATION)Proof sketch. The proof relies on the fact that the workload matrices can be lower bounded componentwise by simpler matrices: A α,β ≥ A α,0 and A α,0 ≥ Id. For the simpler matrices, the bounds (19) can then be derived analytically, and the general case follows by monotonicity.Discussion. Analogously to the case of single participation, Theorems 6 to 9 again establish that the proposed BSR is asymptotically superior to the baseline factorizations if α = 1. A comparison of Theorems 6 and 7 suggests that, at least for maximal participation, k = n b and p = b, the bandedness of the p-BSR improves the approximation quality, specifically in the practically relevant regime where b n. While none of the methods match the lower bound of Theorem 6, we conjecture that this is not because any asymptotically better methods would exist, but rather a sign of Equation ( 18) is not tight. Both theoretical consideration and experiments suggest that a term linear in k should appear there. For α < 1, all studied methods are asymptotically identical and, in fact, optimal. 
experiments: 
 To demonstrate that BSR can achieve high accuracy not only in theory but also in practice, we compare it to AOF and baselines in numerical experiments. Our results show that BSR achieves quality comparable to the AOF, but without the computational overhead, and it clearly outperforms the baseline factorizations. The privacy guarantees are identical for all methods, so we do not discuss them explicitly.Implementation and computational cost. We implement BSR by the closed-form expressions of Theorem (1). For single data participation, we use the square root decomposition directly. For repeated data participation we use p-BSR with p = b. Using standard python/numpy code, computing the BSR as dense matrices are memory-bound rather than compute-bound. Even sizes of n = 10, 000 or more take at most a few seconds. Computing only the Toeplitz coefficients is even faster, of course.To compute AOF, we solve the optimization problem (4) using the cvxpy package with SCS backend, see Algorithm B for the source code3 . With the default numerical tolerance, 10 -4 , each factorization took a few minutes (n ≤ 100) to hours (n ≤ 500) to several days (n ≥ 700) of CPU time. Note that this overhead reappears for any change in the number of update steps, n, weight decay, α, or momentum, β, as these induce different workload matrices. In our experiments, when the optimization for AOF did not terminate within 10 days, we reran the optimization problem with the tolerance increased by a factor of 10. The runtime depends not only on the matrix size but also on the entries. In particular, we observe matrices with momentum to be harder to factorize than without. For large matrix sizes we frequently encountered numerical problems: the intermediate matrices, S, in (4), often did not fulfill the positive definiteness condition required to solve the subsequent Cholesky decomposition for C. Unfortunately, simply projecting the intermediates back to the cone of positive semi-definite matrices is not enough, because the resulting C matrices also have to be invertible and not too badly conditioned. Ultimately, we adopted a postprocessing step for S that ensures that all its eigenvalues were at least of value 1/n, which we find to be a reasonable modification to ensure the stability of the convergence. Enforcing this empirically found value leads to generally good results, as our experiments below show, but it does add an undesirable extra level of complexity to the process. In contrast, due to its analytic expressions, BSR does not suffer from numerical problems. It also does not possess additional hyperparameters, such as a numeric tolerance or the number of optimization steps.Apart from the factorization itself, the computational cost of BSR and AOF are nearly identical. Both methods produce (banded) lower triangular matrices, so computing the inverse matrices or solving linear systems can be done within milliseconds to seconds using forward substitution. Note that, in principle, one could even exploit the Toeplitz structure of p-BSR, but we found this not to yield any practical benefit in our experiments. Computing the sensitivity is trivial for p-BSR using Corollary 1, and it is still efficient for AOF by the dynamic program proposed in Choquette-Choo et al. [2023a] .As a first numeric experiment, we evaluate the expected approximation error for workload matrices that reflect different SGD settings. Specifically, we use workload matrices (9) for n ∈ {100, 200, . . . , 1000, 1500, 2000}, with α = {0.99, 0.999, 0.9999, 1}, and β ∈ {0, 0.9}, either with single participation, k = 1, or repeated participation, b = 100, k = n/100. The results confirm our expectations from the theoretical analysis: in particular, BSR's expected approximation error is quite close to AOF's, typically within a few percent (left plot). Both methods are clearly superior to the naive factorizations. For large matrix sizes, BSR sometimes even yields slightly better values than AOF (right plot). However, we believe this to be a numeric artifact of us having to solve AOF with less-than-perfect precision.Private Model Training on CIFAR-10. To demonstrate the usefulness of BSR in practical settings, we follow the setup of Kairouz et al. [2021] and report results for training a simple ConvNet on the CIFAR-10 dataset (see Table 1 in Appendix C for the architecture). Specifically, we adapt Google's reference implementation of DP-SGD in jax Bradbury et al. [2018] In both cases, 20% of the training examples are used as validation sets to determine the learning rate η ∈ {0.01, 0.05, 0.1, 0.5, 1}, weight decay parameters α ∈ {0.99, 0.999, 0.9999, 1}, and momentum β ∈ {0, 0.9}. Figure 2 shows the test set accuracy of the model trained with hyperparameters that achieved the highest validation accuracy. 4 One can see the expected effect that in DP model training, more update steps/epochs do not necessary lead to higher accuracy due to the need to add more noise. The quality of models trained with BSR is mostly identical to AOF. When training for a large number of epochs it achieves even better slightly results, but this could also be an artifact of us having to solve AOF with reduced precision in this regime. Both methods are clearly superior to the baselines. 
conclusion and discussion: 
 We introduce an efficient and effective approach to the matrix factorization mechanism for SGDbased model training with differential privacy. The proposed banded square root factorization (BSR) factorization achieves results on par with the previous state-of-the-art, and clearly superior to baseline methods. At the same time, it does not suffer from the previous method's computational overhead, thereby making differentially private model training practical even for large scale problems.Despite the promising results, some open questions remain. On the theoretical side, the asymptotic optimality of BSR without weight decay is still unresolved because the current upper bounds on the expected approximation error do not match the provided lower bounds. Based on the experimental results, we believe this discrepancy lies with the lower bounds, which we suspect should be linear in the number of participations. We observe that BSR achieves results comparable to AOF, although we cannot currently prove this due to the insufficient understanding of AOF's theoretical properties; nonetheless, we consider it a promising research direction. On the practical side, it would be interesting to extend the guarantees to even more learning scenarios, such as variable learning rates. 
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] 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: [Yes] Justification: We have cited all the papers which have results and code used in our paper.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] 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. Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • 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. 
 
 Appendix: A General introduction to differential privacy.Differential privacy [Dwork, 2006] is a robust framework designed to provide strong privacy guarantees for statistical analyses and data sharing. It aims to protect individual data points in a dataset while still allowing meaningful aggregate information to be extracted. Unlike traditional data anonymization techniques, which might involve removing identifiers or aggregating data, differential privacy offers a mathematical definition of privacy that quantifies the amount of privacy loss and ensures that the risk of identifying any individual's data remains low, even when combined with other data sources.To formalize this concept, a randomized mechanism M is said to provide (ε, δ)-differential privacy if, for all data sets D and D that differ in one element, and for all subsets of the mechanism's output space S:Pr[M (D) ∈ S] ≤ e ε • Pr[M (D ) ∈ S] + δFor a detailed introduction to differential privacy, we recommend the books "The Algorithmic Foundations of Differential Privacy" by Dwork and Roth [2014] and "The Complexity of Differential Privacy" by Vadhan [2017] .Algorithm 2 Source code for computing AOF using cvxpy.import cvxpy as cp import numpy as np def banded_factorization (A , b ) :n = len ( A ) X = cp . Variable (( n , n ) , PSD = True ) # cp . matrix_frac (A , X ) = tr ( A . T @ X ^-1 @ A ) objective = cp . Minimize ( cp . matrix_frac ( A .T , X ) * np . ceil ( n / b ) ) constraints = [ cp . diag ( X ) == 1] for i in range (b , n ) : constraints += [ cp . diag (X , i ) == 0 , cp . diag (X , -i ) == 0] prob = cp . Problem ( objective , constraints )return prob . solve ( solver = ' SCS ') , X . value Conv2D(channels=32, kernel=(3, 3), strides=(1, 1), padding='SAME', activation='relu') Conv2D(channels=32, kernel=(3, 3), strides=(1, 1), padding='SAME', activation='relu') MaxPool(kernel=(2, 2), strides=(2, 2)) Conv2D(channels=64, kernel=(3, 3), strides=(1, 1), padding='SAME'), activation='relu') Conv2D(channels=64, kernel=(3, 3), strides=(1, 1), padding='SAME'), activation='relu') MaxPool(kernel=(2, 2), strides=(2, 2)) Conv2D(channels=128, kernel=(3, 3), strides=(1, 1), padding='SAME'), activation='relu') Conv2D(channels=128, kernel=(3, 3), strides=(1, 1), padding='SAME'), activation='relu') MaxPool(kernel=(2, 2), strides=(2, 2)) Flatten() Dense(outputs=10)In this section we provide additional experiments comparing BSR, AOF and baselines: Figures 3 and 4 and following tables show their expected approximation error (lower is better) for workload matrices stemming from SGD with different hyperparameter settings. Figure 5 ) and following tables show the accuracy of resulting classifiers on CIFAR-10 (higher is better) for different privacy levels.The results show the same trends as the one in Section 4. BSR achieves almost identical expected approximation error as AOF, and results in equally good classifiers. In some cases, results for BSR even improve over AOF's. Presumably this is because of numerical issues in solving the optimization problem for AOF. 50 29.8 ± 2.7 31.5 ± 1.4 30.2 ± 1.5 20.8 ± 0.8 100 31.2 ± 0.9 34.0 ± 1.5 29.8 ± 1.2 16.5 ± 2.9 200 30.9 ± 1.4 31.5 ± 1.9 26.0 ± 1.6 15.9 ± 2.3 400 32.1 ± 0.9 33.3 ± 0.8 24.8 ± 2.4 14.7 ± 2.0 500 31.6 ± 2.3 32.8 ± 1.8 23.8 ± 1.4 14.9 ± 2.0 1000 30.2 ± 1.6 30.1 ± 1.1 17.5 ± 0.7 12.5 ± 1.5 2000 29.9 ± 2.0 30.5 ± 1.2 16.1 ± 1.0 11.3 ± 1.1 Table 20 : Numeric values for results in Table 5 top left plot (CIFAR-10, single participation, ( , δ) = (2, 10 -5 ).Id A 1 31.2 ± 0.9 34.0 ± 1.5 29.8 ± 1.2 16.5 ± 2.9 2 34.5 ± 1.0 36.5 ± 0.6 30.1 ± 0.9 18.6 ± 2.1 3 36.8 ± 2.3 38.0 ± 1.2 31.3 ± 1.4 17.6 ± 1.6 4 39.2 ± 2.1 40.1 ± 1.6 28.8 ± 0.9 17.4 ± 2.1 5 40.4 ± 1.9 41.0 ± 0.9 28.4 ± 1.0 16.1 ± 2.1 6 41.3 ± 0.7 41.8 ± 1.3 28.3 ± 0.9 15.2 ± 1.3 7 42.2 ± 1.6 44.1 ± 0.9 28.4 ± 1.0 15.4 ± 1.9 8 43.3 ± 1.2 44.4 ± 1.1 28.9 ± 1.4 16.1 ± 1.7 9 44.1 ± 0.4 45.2 ± 0.5 28.6 ± 0.5 15.6 ± 1.6 10 44.3 ± 1.8 43.8 ± 0.3 27.7 ± 1.5 15.8 ± 0.9 15 44.9 ± 1.1 43.9 ± 0.8 27.3 ± 1.1 14.9 ± 1.0 20 45.7 ± 1.2 43.0 ± 1.3 26.9 ± 1.5 14.2 ± 0.9 In this section we report on experimental results when different optimizers are used to (approximately) solve the AOF optimization problem (4). Besides cvxpy (CVX) these are standard gradient descent (GD) and the Limited-Memory Broyden-Fletcher-Goldfarb-Shanno algorithm (LBFGS). The latter two we implement in jax using the optax toolbox. Similar to [Granqvist et al., 2024, ftrl_mechanism .py], we use an adaptive line-search for the step size of the gradient-based methods, which at the same time ensures the positive definiteness constraints of the optimization problem. Our implementation differs from theirs, however, in that our learning rate is not restricted to shrink monotonically, thereby avoiding premature termination.We report the runtimes for the different methods in Tables 24 to 39 . For comparison, we also include results for BSR and the CVX optimizers with three tolerance levels in the same settings, where practically feasible. Note that while the experiments for BSR and CVX used a single-core CPU-only environment, the experiments for GD and LBFGS were run on an NVIDIA H100 GPU with 16 available CPU cores. As a consequence, the absolute runtimes are not directly comparable between the methods, but they should rather be seen as illustrations of the scaling behavior of the method for different workload types and problem sizes. Indeed, the results show a clear trend: BSR is the fastest, with almost no overhead. Even for the largest problem sizes of n = 10 000, BSR never took more than 2.5s to despite running in the single-core CPU-only setup. GD and LBFGS benefit strongly from the GPU hardware. In the multiple participation setting (p = 100, k = n/p), they solve most workload sizes within a few minutes, except the largest ones, which for GD can take a few hours. In the single participation setting (k = 1), LBFGS also occasionally need several hours to converge. In general, stronger weight decay (smaller α) tends to lead to lower runtimes, while the use of momentum (β = 0.9) to higher times until convergence. CVX (on weak hardware) is orders of magnitude slower than the other methods. Furthermore, its runtime grow approximately cubic with the problem size, whereas for GD and LBFGS the relation is not too far from linear. Note that despite the stable patterns described above, all runtime results should be taken with caution, because internal parameters of the optimization, such as the convergence criterion and the specific implementation of the line search can substantially influence the overall runtime as well.Figures 6 and 7 report the expected approximation errors achieved by the different optimizers of AOF (4) and by BSR. For CVX, we report the smallest error value across all tolerance levels for which the optimization converged.The curves show several clear trends. GD and LBFGS generally perform similarly, and achieve expected approximation errors slightly (at most a few percent) lower than BSR. An exception are the problems with momentum (β = 0.9) in the single participation setting, where it appears that SGD occasionally fails to find the optimum for large problem sizes (n ≥ 1500). CVX performs comparably to the other methods for problems without momentum (β = 0). With momentum, however, the solutions it find are often worse than the other methods, especially in the single-participation setting and for medium to large problem sizes (n ≥ 500). Presumably, even smaller tolerance values would be require here, which, however, would result in even longer runtimes. In this section, we provide the complete proofs for our results from the main manuscript. For the convenience of the reader, we also restate the statements themselves.Theorem 1 (Square-Root of SGD Workload Matrix). Let A α,β be the workload matrix (9). Then A α,β = C 2 α,β for C α,β = LDToep(c 0 , . . . , c n-1 ), with c 0 = 1 and c j = j i=0 α j-i r j-i r i β i for j = 1, . . . , n -1 with coefficients r i = -1/2 i . For any p ∈ {1, . . . , n}, the p-banded BSR matrix C |p| α,β is obtained from this by setting all coefficients c j = 0 for j ≥ p.Proof. We observe that A α,β can be written asA α,β =     1 0 . . . 0 α 1 . . . 0 . . . . . . . . . . . . α n-1 α n-2 . . . 1     ×     1 0 . . . 0 β 1 . . . 0 . . . . . . . . . . . . β n-1 β n-2 . . . 1     =: E α × E β (20)Relying on the result from Henzinger et al. [2024] , that E1/2 1 =     1 0 . . . 0 r 1 1 . . . 0 . . . . . . . . . . . . r n-1 r n-2 . . . 1     with r k = -1/2 k, one can check that the square roots of the matrices E α , E β are:EQUATION)An explicit check yields that these matrices commute, i.e.E 1/2 α E 1/2 β = E 1/2 β E 1/2 α . Therefore C α,β = A 1/2 α,β = E 1/2 α × E 1/2 β =     1 0 . . . 0 c 1 1 . . . 0 . . . . . . . . . . . . c n-1 c n-2 . . . 1     , with c k = k i=0 α i r i r k-i β k-i . (22) F.2 Proof of Theorem 2Theorem 2 (Sensitivity for decreasing non-negative Toeplitz matrices). Let M = LDToep(m 0 , . . . , m n-1 ) be a lower triangular Toeplitz matrix with decreasing non-negative entries, i.e.EQUATION)where M [•,1+jb] denotes the (1 + jb)-th column of M .Proof. Because all entries of M are positive, so are the entries of M M . Therefore, the condition is fulfilled such that (2) holds with equality, andEQUATIONwhereΠ k,b = { π ⊂ {1, . . . , n} : |π| ≤ k ∧ ({i, j} ⊂ π ⇒ i = j ∨ |i -j| ≥ b) }.We now establish that {1, 1 + b, . . . , 1 + (k -1)b} is an optimal index set, which implies the statement of the theorem.To see this, let π * be any optimal solution and let i * ∈ π * be a column index that is not as far left as possible, i.e., π * \ {i * } ∪ {i * -1} would be a valid index set in Π k,b . If such an index i * exists, we split π * into the left indices, which are smaller than i * , and the remaining, right, ones:π * = π * L ∪π * R with π * L = {i | i ∈ π ∧ i < i * }, π * R = {i | i ∈ π ∧ i ≥ i * }.Then, we construct a new index set in which the left indices are kept but all right ones are shifted by one position to the left:π = π * L ∪ ← - π * R with ← - π * R = {i -1 | i ∈ π * R }.By the condition on i * , we know π ∈ Π k,b . We now prove that f (π , π ) ≥ f (π, π), so π must also be optimal. First, we observe two inequalities: for any i, j > 1:EQUATIONand for i ≥ 1, j > 1:EQUATIONwhere the last inequality holds because by assumption m l-i-1 ≥ m l-i for l ≥ i + 1. Now, we split f (π , π ) and f (π , π ) into three terms: the inner products of indices below i * , the ones of terms above i * and the ones between both,EQUATION)The first term appears identically in both expressions. The second term fulfillsEQUATIONby Equation ( 24). The third term fulfillsEQUATIONby Equation ( 28). In combination, this establishes f (π , π ) ≥ f (π * , π * ), and since π * was already optimal, the same must hold for π .Using the above construction, we can create a new optimal index sets, π * , until reaching one that does not contain any index i * as described anymore. Then π * = {1, 1+b, . . . , 1+(l -1)b} for some l ∈ N must hold. If l = k, the statement of Theorem 2 is confirmed. Otherwise, π = {1, . . . , 1 + (k -1)b} is superset of π * , so because of the positivity of entries, f (π , π ) ≥ f (π * , π * ) must hold. Once again, because π * was optimal, the same must hold for π , which concludes the proof.Corollary 1. The sensitivity of the p-BSR for SGD can be computed using formula (10).Proof. From (1) we know that C α,β is a Toeplitz matrix with coefficients (1, c 1 , . . . , c n-1 ), wherec j = j i=0 α i r i r j-i β j-i for 0 ≤ β < α ≤ 1, with r i = | -1/2 i | = Bi 4 i , where B i = 2i i is the i-central binomial coefficient. It suffices to show that c j ≥ c j+1 for any j ∈ {1, . . . , n -1}.First, we show for the r i coefficients:EQUATIONwhereEQUATIONUsing the classic identity between Catalan numbers, 2 [Batir et al., 2021, Identity 4 .2] we obtaink i=0 B k-i C i = B k+1 , e.g.EQUATIONwhere the last inequality follow from the fact that β -i -β > 0 for each i = 0, . . . , k and any β < 1. This proves the monotonicity of c k .For α < 1, we observe that c j = α j j i=1 r i r k-i γ j-i for γ = α β < 1. Clearly, the sequence α j is decreasing, and by the above argument, the sum is decreasing, too. Consequently, c j is the product of two decreasing sequences, so it is also decreasing, which concludes the proof.Before the remaining proofs, we establish a number of useful lemmas. Lemma 2. For any C ∈ R n×n with C C ≥ 0 it holds for any b ∈ {1, . . . , n} thatEQUATIONwhereC 2 2,∞ = max i=1,...,n C [•,i] 2 .Proof. This follows directly from Theorem 2:EQUATIONProof. We first show the upper bound. Observe that for any π ⊂ [n]:i,j∈π[C C] i,j = i,j∈π C [•,i] , C [•,j] ≤ i,j∈π C [•,i] C [•,j] (45) = ( i∈π C [•,i] ) 2 ≤ |π| i∈π C [•,i] 2 ≤ |π| C 2 F (46)Therefore, using Theorem 2:EQUATION)For the lower bound, we introduce some additional notation. Let Πk be the set of b-separated index sets with exactly k elements. Then, from Theorem 2, we obtain EQUATIONEQUATIONThis proves the statement (49), because 1 b = k n in this case. As an induction step, we prove that if (49) holds for some value k ≤ n b , then it also holds for k -1 ≥ 1.Let π * ∈ argmax π∈ Πk S(π) and j * = argmin j∈π * S({j}), such that we know that S({j}) ≤EQUATIONwhere in the last step we used the induction hypothesis. This concludes the proof.Lemma 4. For C α,β as in (1), and k = 1, it holds that1 n n j=1 n-j i=0 c 2 i ≤ E(C α,β , C α,β ) ≤ n-1 i=0 c 2 i (52)Proof. From Lemmas 3 and 2 we obtainEQUATIONwhere the last identify follows from the explicit form of C α,β . The lower bound follows fromEQUATION)and again the explicit form of C 2 F .Lemma 5. Forr j = | -1/2 j | = 1 4 j 2jj it holds that:EQUATIONProof. The double inequality is a particular case of a more general pair of binomial inequalities when k = j and m = 2j:EQUATIONwhere H(k/m) is the binary entropy function, with H(1/2) = 1. The proof of the general result (56), can be found in MacWilliams and Sloane [1977, Chapter 10, Lemma 7, p309] .Lemma 6. Let c k = k j=0 α j r j r k-j β k-j as in (1). Then c 0 = 1, and for j ≥ 1:EQUATIONProof. We exploit the upper and lower bounds from Lemma 5. First, we write c k = α k k j=0 r j r k-j γ j with γ := β α . Then we check immediately that c 0 = 1 andc 1 = 1 2 (α + β) = α 2 (1 + γ) ≤ α 2 11-γ . For j ≥ 2 we derive the upper bound byEQUATIONEQUATIONwhich proves the upper bound on a j . The lower bound for j ≥ 1 follows trivially fromEQUATIONLemma 7. For j ∈ {1, . . . , n} it holdsEQUATIONfor α = 1, and otherwiseEQUATIONProof. We first prove the result for α = 1. Combining Lemmas 4 and 6 we obtainEQUATIONEQUATIONFor α < 1, if follows analogously:EQUATIONEQUATIONwhere the last term emerges from∞ i=j+1 α 2i i ≥ α 2(j+1) j+1 ∞ i=0 α 2i = α 2(j+1) j+1 1 1-α 2 . Lemma 8. Let 0 ≤ β < α ≤ 1. Let σ 1 ≥ • • • ≥ σ n bethe sorted list of singular values of A α,β . If α < 1, then for j = 1, . . . , n:EQUATIONandEQUATIONIf α = 1, then for j = 1, . . . , n, 2 πEQUATIONand consequentlyEQUATIONProof. The statements on the singular values follow from the following Lemma 9, because A α,β = E α E β . Because E α and E β are diagonalizable and they commute, we haveσ n (E β )σ j (E α ) ≤ σ j (E α E β ) ≤ σ 1 (E β )σ j (E α ).For α < 1 the lower bound follows from A α,β * ≥ trace A α,β , and the upper bound follows from the identity A α,β * = n j=1 σ j . For α = 1, the bounds follow from the same identity together with the fact thatEQUATIONEQUATION)Proof. We follow the steps of SebastienB [2017] , and use that the singular values of E t are the reciprocals of the singular values of E -1 t , which themselves are the eigenvalues of (E t ) -1 ((E t ) -1 ) =: T , i.e., for i = 1, . . . , n:EQUATIONThe E -1 t and T can be computed explicitly asEQUATIONLemma 10. All eigenvalues, µ, of T fulfillEQUATIONProof. By Gershgorin's circle theorem [Gershgorin, 1931] , we know thatµ fulfills i) |1 -µ| ≤ t, i.e. 1 -t ≤ t ≤ µ ≤ 1 + t or ii) |1 + t 2 -µ| ≤ 2t, i.e. 1 -2t + t 2 ≤ µ ≤ 1 + 2t + t 2 . For t ∈ [0, 1]the first condition implies the second, so (79) must hold.Case I: For t < 1, the statement (76) follows from Lemma 10 in combination with (77).Case II: For t = 1 the matrix simplifies toT =       1 -1 0 . . . 0 -1 2 -1 . . . 0 0 . . . . . . . . . . . . 0 . . . -1 2 -1 0 . . . 0 -1 2      . Note that T is not exactly Toeplitz, because of the top left entry, so closed-form expressions for the eigenvalues of tridiagonal Toeplitz matrices do not apply to it. Instead, we can compute its eigenvalues explicitly. Matrices of this form have been studied by Elliott [1953] ; for completeness, we provide a full proof here.Let µ be an eigenvalue of T with eigenvector Ψ = (Ψ 0 , . . . , Ψ n-1 ). From the eigenvector equation T Ψ = µΨ we obtainEQUATIONwhich yields a linear recurrence relationEQUATIONwith two boundary conditionsEQUATION)We solve the recurrence relation using the polynomial method [Greene and Knuth, 1990] . The characteristic polynomial of ( 83) is P (z) = z 2 + (µ -2)z + 1. Its roots areEQUATIONfor some value θ ∈ [0, 2π). Note that the expression under the second square root is positive, because of Lemma 10. The last equation is a consequence of,EQUATIONFrom standard results on linear recurrence, it follows that any solution to (83) has the form Ψ j = c 1 (r + ) j + c 2 (r -) j for some constants c 1 , c 2 ∈ C. The fact that Ψ j must be real-valued implies that c 1 = c 2 =: αe iφ for some values α ∈ R, φ ∈ [0, 2π). Dropping the normalization constant (which we could recover later if needed), we obtainEQUATIONNext, we use the boundary conditions to establish values for φ and θ.Equation ( 85) can be rewritten asEQUATIONwhich, using 2 cos(α + β) = cos(a + b) + cos(a -b), simplifies toEQUATIONConsequently, φ + nθ = 1 2 π + kπ must hold for some k ∈ N. Equation ( 84) can be rewritten asEQUATIONwhich simplifies toEQUATIONOne solution to this would be θ = 0, but that would implies µ = 0, which is inconsistent with T being an invertible matrix. So instead, it must hold that φ = θ 2 + kπ for some k ∈ N. Combining both conditions and solving for θ we obtainEQUATION)Each such value θ k for k ∈ {0, . . . , n -1} yields an eigenvector with associated eigenvalue µ = 2 -2 cos θ k = 4 sin 2 (θ k /2). Now, (76) follows from this in combination with (77).Theorem 3 (Expected approximation error with single participation). Let A α,β ∈ R n×n be the workload matrix (9) of SGD with momentum 0 ≤ β < 1 and weight decay parameter 0 < α ≤ 1, where α > β. Assume that each data item can contribute at most once to an update vector (e.g. single participation, k = 1). Then, the expected approximation error of the square root factorization,A α,β = C 2 α,β , fulfills 1 ≤ E(C α,β , C α,β ) ≤ 1 (α -β) 2 log 1 1 -α 2 (11)for α < 1, andEQUATION)Proof. The proof consists of a combination of Lemmas 4 and 7. Because in the single participation k = 1, so we need just the first column of matrix C α,β :EQUATION)which proves the upper bounds. For the lower bounds, for any α ≤ 1:EQUATIONAlso, for α = 1:EQUATIONbecause log( (n + 1)! ) ≥ (n + 1) log(n + 1) -n.F.6 Proof of Theorem 5Theorem 5. Assume the setting of Theorem 3. Then, the baseline factorizationsEQUATIONProof. For α = 1, by Lemma 2 we have:EQUATION)For α < 1:EQUATIONTogether withEQUATIONas β < 1 and the sumn-1 j=0 jβ j+1 is uniformly bounded by β ∞ j=0 jβ j = β 2 (1-β) 2For α < 1:EQUATIONwhere the second equality due to the fact that we average over the sequence jx j+1 which converges to 0 for |x| < 1. EQUATIONwhere c j = j i=1 a i b n+1-i so the LTT structure is preserved with multiplication that allows us to work with sequences and their convolutions rather than matrix multiplication. For instance, we would write the previous product in the form:EQUATION)The inverse of the Lower Triangular Toeplitz matrix remains a Lower Triangular Toeplitz (LTT) matrix because we can find a unique sequence (c 1 , . . . , c n ) such that:EQUATIONEQUATIONwith the restriction that a 1 = 0; otherwise, the original matrix was not invertible. We consider the banded square root factorization A α,β = B EQUATIONBy induction argument, we can see that b p+j = 2c p+j for 0 ≤ j ≤ p -1. For the remaining n -2p values we will prove convergence to a constant.EQUATION)We make an ansatz for the solution of the linear recurrence in the form:EQUATIONwhere y j represents the error terms, which will be proven to converge to 0. The sequence y j satisfies the following recurrence formula:EQUATION)We denote w j = c j α -j which is a decreasing sequence because the values correspond to the C 1,β/α matrix. We rewrite the recurrence in matrix notation: EQUATIONEQUATION)To show that the error terms y j goes to 0 as j goes to infinity, we first study the characteristic polynomial of the associate homogeneous relations: as:g(λ) = λ p-1 + w 1 λ p-2 + • • • + w p-2 λ + w p-1 . (119) Because 1 > w 1 > w 2 > • • • > w b-1 > 0, itEQUATIONWe use the following lower bound for the sum of w j :EQUATIONCombining these bounds we can upper bound the Frobenious norm of the matrix B |p| α,β the following way:EQUATIONSimplifying for the leading terms in asymptotics, we have: EQUATIONEQUATIONEQUATIONCombining sensitivity with the upper bound for the Frobenious norm we obtain:EQUATIONF.8 Proof of Theorem 7 for Square Root Factorization Theorem 7 (Approximation error of Square Root Factorization). Let A α,β ∈ R n×n be the workload matrix (9) of SGD with momentum 0 ≤ β < 1 and weight decay 0 < α ≤ 1, with α > β. Let A α,β = C 2 α,β be its square root factorization. Then, for any b ∈ {1, . . . , n} and k = n b it holds:EQUATIONProof. We prove the case without weight decay (α = 1) and with weight decay (α < 1) separately.Case 1) no weight decay (α = 1).We start by bounding the b-min-separation sensitivity:EQUATIONConsider a scalar product for a general pair of indices, j > i:EQUATIONUsing the bounds on c k (6) for α = 1 we can lower and upper bound this sum by:EQUATION)We can compute the indefinite integral explicitly:dx x(j -i + x) = F j -i x + C (130) for F (a) = 2 log 1 a + 1 + 1 a .In combination, we obtain the upper and lower bound for (126):EQUATIONNow we are ready to bound the sensitivity of the matrix C 1,β :EQUATIONand, analogouslyEQUATIONFirstly, using ( k e ) k ≤ k! ≤ k k , we bound the sum of the logarithms:EQUATION)To upper bound the last term in sensitivity lower bound (134), we use the auxiliary inequalityf (a) = 2 log 1 a + 1 + 1 a ≤ 4 √ a to derive: 1 2 0≤i<j≤k-1 f (b(j -i)) ≤ 2 √ b 0≤i<j≤k-1 1 √ j -i = 2 √ b k-1 j=1 j t=1 1 √ t ≤ 4 √ b k-1 j=1 j ≤ 8 3 √ b k 3/2 (137)To bound the final term we establish the following inequalities for f (a):EQUATIONThen we can bound the first double sum in (134) asEQUATION)To bound the term 0≤i<j≤k-1 log k-i j-i we use the following identities:EQUATIONNow, using that x log x is a monotonically increasing function,EQUATIONAs a lower bound, we obtain0≤i<j≤k-1 log k -i j -i ≥ 2 k-1 1 x log xdx -k log k + k log k -k(k -1) log(k -1) (149) = (k -1) 2 log(k -1) - (k -1) 2 2 + k log k -k(k -1) log(k -1) (150) = -(k -1) log(k -1) - (k -1) 2 2 + k log k (151) ≥ - k 2 2 (152)Therefore, combining the upper bound (146) and the lower bound (149) yieldsEQUATIONEQUATIONCombining all three terms together we obtain the following bounds for the squared sensitivity (134):sens 2 k,b (C 1,β ) ≤ k (1 -β) 2 (log n + 1) + 6 (1 -β) 2 k 2 (155) sens 2 k,b (C 1,β ) ≥ k 4 (log n -1) - 8 3 √ b k 3/2 + 2 5 k 2 (156)Now, we recall the bounds for the Frobenius norm of the matrix C 1,β 7 and (96):EQUATION)With the auxiliary inequality aEQUATIONcombining (155) and the bounds on Frobenius norm ( 7) and ( 96) we get that:EQUATIONMaking the lower bound well-defined requires n ≥ 45, otherwise one can simply takeE(C 1,β , C 1,β ) ≥ 1.As a final step, we combine both inequalities in the following asymptotic statement:EQUATIONwhich concludes the proof of the case without weight decay.Case 2) with weight decay (α < 1).As above, we first express the b-min-separation sensitivity of the matrix C α,β in terms of inner products,EQUATIONand then consider a scalar product for a general pair of indexes j > i:EQUATIONNow, we use the bounds on c t from Lemma 6 for α < 1, to upper and lower bound this sum with the following expression, where γ = β α :EQUATION)We substitute these bounds into Equation ( 161) to obtain the following upper bound for sensitivity of matrix C α,β :EQUATIONwhere the second inequality is due to Equation (95).A lower bound for the sensitivity follows directly from Lemma 3:sens 2 k,b (C α,β ) ≥ k C α,β F ≥ k.(168) The Frobenius norm of the matrix C α,β is the same as that for one round participation; thus, we could reuse Inequalities (95):EQUATION)By merging the bounds for sensitivity and Frobenius norm, we derive the following bounds for error:EQUATIONThe combination of these results yields the following asymptotic statement: EQUATIONProof. Let A α,β = BC be any factorization with CC ≥ 0. From Lemma 3 it follows that EQUATIONAgain, we use the fact that A α,0 ≤ A α,β ≤ 1 α-β A α,0 . Now A α,0 fulfills the conditions of Theorem 2, so from Equation (10) we knowEQUATIONWe first study (176) for α = 1. Then, from the explicit structure of A 1,0 = LDToep(1, 1, . . . , 1) one sees that the vectors inside the norm have a block structure EQUATIONEQUATIONbecause bk ≥ n. ConsequentlyEQUATIONwhich concludes the proof of this case. For α < 1, A α,β ≥ Id componentwise readily impliesEQUATIONwhich implies the statement of the theorem.• 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).Note that this specific type of MF should not be confused with other occurrences of matrix factorization in, potentially private, machine learning, such as in recommender systems[Shin et al.,2018 , Li et al., 2021]].38th Conference on Neural Information Processing Systems (NeurIPS 2024).As proved inDenisov et al. [2022, Theorem 2.1], establishing privacy in this non-adaptive setting suffices to guarantee also privacy in the adaptive setting, where the update vectors depend not only on the data items but also the intermediate estimates of the model parameters, as it is the case for private model training.Additional experiments with gradient-based optimizers can be found in Appendix E.Such a setting would not optimal for real-world private training, because the many repeated experiments reduce the privacy guarantees[Papernot and Steinke, 2021, Kurakin et al., 2022, Ponomareva et al., 2023]. We nevertheless adopt it here to allow for a simpler and fair comparison between methods. 


Decision: This paper studies the problem of optimal matrix factorization mechanism for differentially private training. Specifically, the authors develop a banded square root factorization approach and provide formal lower and upper bounds on the approximation error.

All the reviewers agree that the paper presents solid results and makes significant contributions to the field of differentially private training. In light of this, I recommend acceptance and encourage the authors to incorporate the reviewers' suggestions into the revision.

Reviews: review: paper_summary: paper_summary: The paper considers the problem of adding correlated noise $C^{-1}z$ instead of independent noise $z$ across iterations in continual counting (equivalently, DP-SGD). Past work gave an algorithm to choose $C^{-1}$ that optimizes some objective on the noise under b-min-separated participation, but this algorithm requires solving an SDP which is infeasible when using a large number of iterations, i.e. optimizing over $C$ that has a large number of rows / columns $n$. The authors propose a choice of $C$, which they call banded square root, which has a succinct representation that is efficiently computable in time independent of $n$. Furthermore, their choice of $C$ is banded, i.e. is non-zero below the $b$-th diagonal for $b \geq 1$, which gives it some nice properties such as being able to compute $C^{-1}z$ efficiently in a streaming manner. To define the banded square root $C$, they authors take $A$ which represents the workload, effectively a representation of the updates in SGD which can include momentum and weight decay, and take its matrix square root. They then truncate the matrix square root to its first $b$ diagonals. The authors give an explicitly and efficiently calculable formula for the computing the first column of this matrix, which specifies the whole matrix as it is Toeplitz. The authors analyze the asymptotic error guarantees of different choices of $C$, including the banded square root, in the single- and multiple-participation settings. Specifically, their error guarantees show that banded square root improves on $C = I$ or $C = A$, and nearly matches a lower bound on any factorization if using enough bands. The authors conduct experiments training models for classification on CIFAR10 and show that banded square root is competitive with the choice of $C$ given by solving a more expensive SDP.
strength_weakness: strength_weakness: * While the previous work did require solving an SDP, to my understanding this is for the case where $C$ can be is only required to be a banded lower-triangular matrix. If $C$ is required to be Toeplitz as well, as the banded square root factorization is, then one only needs to optimize over $b$ variables instead of $\approx nb$ and it is not clear that the computations in the previous work are expensive, which mitigates the improvements in this paper.
* It is worth pointing out there is a work of https://arxiv.org/abs/2404.16706, that also gives an efficient-to-compute way to choose $C$, and leads to faster streaming noise generation for than banded matrices. This work I would consider concurrent, so I did not account for any possible overlap between the two papers' results and impact when assigning my score.
questions: questions: * wrt the first weakness, did the authors consider this alternate approach / if so, do you believe banded square root still offers speedups in this setting?
* The previous work of Choquette-Choo et al. shows that banded $C$ are compatible with privacy amplification with sampling. Have you considered combining your error analysis with their privacy analysis (e.g., maybe for RDP to simplify) in the setting where batches are sampled in DP-SGD? It might be an interesting question to see how the number of bands affects error in the presence of amplification, since they observed more bands reduce the benefits of amplification.
limitations: limitations: Yes
paper_summary: paper_summary: This paper addresses optimal matrix factorization mechanism for differential privacy, focussing on stochastic gradient descent (SGD) optimization.
It introduces the banded squared root (BSR) factorization and provide formal lower and upper bounds to the approximation error.
The BSR achieves competitive formal approximation error compared to state-of-the-art methods and demonstrates practical utility in real-world training tasks. 
The authors provide theoretical bounds for single and repeated participation training, showing that BSR maintains high training accuracy with provable privacy guarantees while being computationally efficient.
strength_weakness: strength_weakness: - The experimental evaluation is limited to a small set of real-world and synthetic experiments.
- The range of synthetic datasets used in the experiments is restricted, which may not be representative of diverse scenarios.
- The data generation process is unclear, and only homogeneous partitions are considered, with no variation in data distributions.
- The real-world experiments are limited to a single dataset, CIFAR-10, which is insufficient to draw general conclusions.
- The post-processing step (line 275) lacks clear explanation, leaving it unclear whether it addresses a principled problem or a numerical issue.
- The experiments do not provide a comparison with non-private approaches, which makes it difficult to understand the practical limitations of the method.
- The evaluation is limited to a single real-world dataset, CIFAR-10, and a simple ConvNet model.
- The experiments do not discuss potential practical limitations of BSR, which may impact its applicability.
- The comparative analysis is limited and does not include a wide range of existing differential privacy methods.
- The formal proof of Theorem 2 (Appendix D.2) employs unclear and undefined notation (Π), which makes it difficult to follow.
- Neither the proof sketch nor the formal proof are straightforwardly understandable, which may hinder comprehension.

Minor:

- The tone of "supposedly optimal" (line 254) could be perceived as inappropriate and may benefit from rephrasing.
- The proof sketch of Theorem 2 refers to the wrong appendix (Appendix D.3 instead of Appendix D.2, presumably).
questions: questions: ./.
limitations: limitations: ./.
paper_summary: paper_summary: This paper proposes the Banded Square Root Matrix Factorization, an efficient approximation of the optimal banded matrix factorization for differentially private machine learning applications.  The authors give a closed form expression for the BSR C matrix for any SGD + Momentum + Weight Decay workload by exploiting the Toepltiz structure of the workload and factors.  These factorizations can be used in place of the prior work and enjoy some efficiency advantages.
strength_weakness: strength_weakness: * The implementation of AOF seems suboptimal and hence the comparison appears to be biased.  Multiple-day runtimes for n ~ 700 seems high when prior work conducts experiments for n >= 2000 (https://arxiv.org/pdf/2306.08153).  Doing a little digging, I found https://www.kdd.org/kdd2016/papers/files/rpp0224-yuanA.pdf, which ran experiments up to n ~ 8192.  The limits of the BSR approach are not really demonstrated besides the statement "Even sizes of n = 10,000 or more take at most a few seconds."  How much beyond this can you scale?  Would be good to add a scalability experiment.
questions: questions: * Assuming your method scales well beyond the prior work of n ~ 8192, can you make a convincing case that in modern DP + ML applications we need to scale beyond that point?  I'm not super familiar with how many iterations DP-SGD is usually run for.

* Thm 2 is an interesting result that allows you to compute sensitivity efficiently under this special Toepltiz structure without requiring bandedness, but it's not clear where you are using this Thm.  Is it better to find a banded factorization with bands = minsep or use a non-banded factorization with Thm 2?
>* Also minor point: might be better to write Eq 10 in terms of m_0, ... m_{n-1} instead of M (or maybe include both expressions).
>* There is an incorrect reference to Eq 10 in the experiments
>* Typo: n_{n-1}

* Regarding the statement: "Apart from the factorization itself, the computational cost of BSR and AOF are nearly identical."  Is the factorization time usually the dominating term in the cost for AOF?  Is that still the case for BSR, even when scaling up well beyond n ~ 8192?  What factor(s) other than n (if any) impact the complexity of these mechanisms?

* It is remarkable that the approximation quality of BSR and AOF is so close.  It looks like you evaluated a few settings, one of them being b=100, k=n/100 in Fig 1.  Does this hold more generally across other settings, and can you demonstrate that with an experiment?  There are two extremes where it I might expect *some* gap: (a) 2 bands (b) n bands.  
>* In Fig 1, the most interesting comparison to me is BSR vs. AOF, the baselines have already been established to be far from optimal in prior work and their presence makes it difficult to perceive difference between the two most interesting lines, would be good to update the plot accordingly.
limitations: limitations: The authors have discussed some limitations, but perhaps adding some of the things I mentioned in this review (or addressing them directly) would be good.
paper_summary: paper_summary: The paper proposes the Banded Square Root (BSR) matrix factorization to speed up banded matrix factorization. The authors demonstrate that the workload matrix for SGD with momentum and weight decay can be expressed as a lower triangular Toeplitz matrix. They utilize explicit recursive expressions to compute the square root of this workload matrix, and Theorem 1 provides a closed-form expression for the square root matrix. Theoretical analysis shows that the calculation of the square root matrix is efficient, and the sensitivity of the mechanism is discussed. Detailed analysis of the approximation error is also provided, demonstrating that the square root factorization has asymptotically optimal approximation quality. Experiments show that BSR can achieve performance comparable to previous state-of-the-art methods.
strength_weakness: strength_weakness: 1. Using CVXPY with SCS to compute AOF is not efficient; L-BFGS can solve this problem more efficiently. For instance, as shown in [1], an implementation using GPUs can solve a problem size of 10,000 in less than an hour. Even with CPUs, the implementation in [2] can solve a problem size of 2,000 in less than an hour, with slight modifications to add the b-min-step constraint.

2. The current experiment only demonstrates results for matrix sizes up to 2000, which is insufficient. A comparison involving matrix sizes of 10,000 or more would be more convincing. How long would it take to obtain a solution for problem sizes of 10,000 and 100,000 using BSR?

[1] https://github.com/apple/pfl-research/blob/2a9d65bd66dc89ef80a24e1be0d28446a0422133/pfl/privacy/ftrl_mechanism.py#L140 

[2] https://github.com/dpcomp-org/hdmm/blob/7a5079a7d4f1a06b0be78019adadf83c538d0514/src/hdmm/templates.py#L391
questions: questions: 1. Enforcing positive semi-definiteness (PSD) for the matrix is crucial for convergence in AOF. In line 274, you mentioned a post-processing step for S to ensure that all its eigenvalues are at least $\sqrt{1/n}$. Could you provide a more detailed description of this step? I couldn't find it in the source code provided in Appendix A.

2. In equation (6), b should be p, right?
limitations: limitations: 1. If the learning rate is not constant, then the workload matrix will not be a Toeplitz matrix, right? This would significantly impact the efficiency of BSR. During the training process, a proper learning rate schedule can enhance convergence. Therefore, a discussion about varying learning rates should be included.

2. The paper does not provide a general introduction to differential privacy, and Figure 2 presents (epsilon, delta) without sufficient background information. This lack of context makes the paper not self-contained.

Comments: comment: Thank you for addressing my concerns. I have raised my rating.

The runtime was observed when I did several experiments. I really encourage to have a comparison experiment on the latest algorithm. Although I agree that the main purpose of the paper is not for runtime comparison.
Thanks for your response. Of course, feel free to pursue any of the directions discussed here, if they interest you. I remain in support of accepting the paper.
Thank you for your feedback and for increasing your rating. We appreciate your insights and will address your recommendations in the revision. Below are responses to your follow-up questions:

 **\> Memory Usage**:  
Yes, running Algorithm 3 with $n=100{,}000$ does require around 40 GB of RAM and $n=500{,}000$ approximately 1 TB. We used a machine with 1.5 TB of RAM for these experiments. However, setting `full='true'` is not necessary for using this mechanism in an ML training pipeline.

**\> Matrix Structures and Efficient Representation**:
The matrices $A$, $B$, and $C^{-1}$ are indeed all Toeplitz. However, $A$ and $B$ are not explicitly needed for MF-SGD, only implicitly (see Alg.1). In our demo code, we instantiate $C$ and $C^{-1}$ explicitly for simplicity, but for large-scale applications, this would not be necessary. The Toeplitz structure allows for efficient representation, avoiding the need to materialize these matrices fully.
Thank you for the response, I went ahead and bumped up my rating -- I think this is a solid paper with nice impact that should be of interest to people working in the space DP ML.   Three things I'd like to see addressed in the revision:

1. I would like to see in the final version the results from Appendix C expanded and included in the main text.  From what I can tell, there is a 10%+ degredation in expected error for BSR in Table 2, which could be meaningful in practice.  Including plots that show how this degredation changes with n and b would be informative and strengthen the experiments.  

2. The results on AOF should also be cleaned up -- make sure you are using an implementation that converges to the optimal solution so you don't end up with strange artifacts where BSR outperforms AOF.  If your custom implementation of AOF is not converging, try modifying an existing implementation, such as one of the links from Reviewer VxsQ below, or finding the source from the kdd paper.  Alternatively asking the authors of arXiv:2306.08153 for their source code might be worth a shot (if you haven't done so already).

3. Add limitations section, discussing challenges to scaling to large n beyond solving Eq (4).

Also, please clarify two follow-up questions:

I'm not sure how you were able to run Algorithm 3 BSR(..., full='true') for n >= 100000 -- wouldn't that require 40 GB of RAM?  Do you need to set full='true' to use this mechanism in an ML training pipeline?

One final comment: It is clear that C has a closed form expression but from Eq (6) it seems like you compute B by materializing A and C, and computing A C^{-1} using something like numpy, which I wouldn't expect to scale well due to the size of the matrices.  Does B also have a Toeplitz structure that you exploit to represent it efficiently?