Multi-Epoch Matrix Factorization Mechanisms for Private Machine Learning

Introduction

Differentially private stochastic gradient descent (DP-SGD) is the de facto standard algorithm for DP machine learning (ML) (Song et al. 2013; Bassily et al. 2014; Abadi et al. 2016a). However, obtaining state-of-the-art privacy-utility tradeoffs critically requires use of privacy amplification techniques like shuffling (Erlingsson et al. 2019; Feldman et al. 2022) or (Poisson) subsampling (Bassily et al. 2014; Zhu & Wang 2019; Wang et al. 2019). These in turn require strong assumptions on the manner in which data is processed that often do not hold under the processing performed by centralized ML pipelines, and are particularly challenging in cross-device federated learning (Kairouz et al. 2021).

Kairouz et al. 2021 recently proposed the DP-FTRL framework that avoids reliance on amplification by sampling, instead leveraging DP streaming of prefix sums (Dwork et al. 2010; Chan et al. 2011; Honaker 2015). DP-FTRL can match (or outperform) DP-SGD in privacy-utility tradeoffs. This algorithm enabled McMahan & Thakurta 2022 to train the first known provably DP ML model on user data in a production setting.

Several works have since focused on this primitive as an instantiation of the streaming matrix mechanism (see Eq. 1) (Henzinger et al. 2022; Fichtenberger et al. 2022; Denisov et al. 2022); in particular, Denisov et al. 2022 showed that leveraging the flexibility inherent in this formulation to design optimal matrices led to significant empirical improvements, though their work was restricted to the single-epoch setting.

This single-epoch restriction is unnatural from the perspective of modern ML, where many passes over the training data are common. We extend matrix factorization mechanisms for ML to the multi-epoch setting by tackling several intertwined problems. This enables state-of-the-art mechanisms for DP-ML, with potentially broader applicability to, e.g., online PCA, marginal estimation, and top-k selection, as discussed in Denisov et al. 2022.

1) We provide a framework for computing the sensitivity of matrix mechanisms under general participation schemas with vector contributions: these are essential to ML applications where we wish to privatize high-dimensional models. However, the efficient computation of optimal matrix factorizations only allows for sensitivity constraints on scalar contributions. In the single participation case a trivial argument equates sensitivity under scalar and vector contributions and the computation of sensitivity is $\mathcal{O}(n)$ by inspection. The situation becomes dramatically more complex under multiple participations. To our surprise, computing even scalar sensitivity can be NP-hard, and the vector-to-scalar reduction does not hold in general (see Section H.2). Nevertheless, we establish sufficient conditions for both computational tractability and the necessary reduction (2.1, H.1), enabling our next contribution:

2) We obtain a closed-form representation of the Lagrange dual function for the loss-minimization problem of computing an optimal matrix factorization subject to multiple-participation sensitivity constraints, in Eq. 6, substantially generalizing the approach of Denisov et al. 2022. This dual formulation is used to efficiently compute optimal factorizations for our experiments.

3) We explore the computational tradeoffs of our approaches. Computing optimal matrix factorizations may become relatively expensive when more than $\approx 10,000$ steps are required. While this is uncommon in federated algorithms for user-level DP, it can be a limitation when training with SGD for example-level privacy. To reduce this cost, we propose and investigate an approach based on the Fast Fourier Transform (FFT) (Cooley & Tukey 1965). Careful analysis of the DFT, shows it is near-optimal for the single-epoch setting and efficiently computable for most, if all, ML settings. Indeed, we find this approach still outperforms the mechanisms from the extant literature, even under multiple participations.

4) We perform detailed empirical comparisons of our mechanisms, e.g., above in Fig. 1. We compare with both the prior matrix mechanism approaches and DP-SGD. We show that the methods proposed here outperform all others (in particular, DP-SGD with amplification), to privacy budgets as low as $\varepsilon\approx 2$, and without any need for privacy amplification. We also find in Fig. 3 that our methods can achieve near the folklore upper bound of full-batch DPGD, but at $340$x less compute. Our code is at: https://github.com/google-research/federated/tree/master/multi_epoch_dp_matrix_factorization.

The core privacy primitive here is the matrix mechanism (Li et al. 2015). Its long history of study and application was, until recently, primarily oriented towards offline, statistical queries (McKenna et al. 2018; Edmonds et al. 2020; Yuan et al. 2016; Hardt & Talwar 2010). Fichtenberger et al. 2022; Denisov et al. 2022 independently applied it to the adaptive streaming setting, where outputs are released one-by-one and privacy analysis must account for an adversary adaptively defining the inputs. Denisov et al. 2022 connected the matrix mechanism to DP ML, via the DP-FTRL algorithm of Kairouz et al. 2021, and showed that computing optimal factorizations significantly improves the privacy-utility-computation tradeoffs when making only a single pass (epoch) over the training data.

Differential Privacy for Adaptive Streams with Multiple Participations

We focus on a generalization of the above two, $(k,b)$-participation, where each example participates at most $k$ times, with any adjacent participations exactly $b$ steps apart: formally, $\Pi$ is the set of all $\pi$ such that $|\pi|\leq k$, and if $\pi=\{i_{1},\dots,i_{k}\}$ indexed in increasing order, we have $\forall j\in\{2,\dots,k\},i_{j}-i_{j-1}=b$. Note $(k{=}1,b{=}n)$-participation recovers the single-epoch setting, and $(k{=}n,b{=}1)$-participation recovers every-step participation, and for example $(k{=}3,b{=}2)$-participation has $\Pi=\{\{1,3,5\},\{2,4,6\}\}$. We focus on this participation schema because of the following three reasons. 1) It encompasses multi-epoch SGD training using a data processing pattern well-supported by modern ML infrastracture. This is in contrast to Poisson or independent fixed-sized batch sampling, with replacement across steps, as is assumed by many works (Abadi et al. 2016a; Bassily et al. 2014; Zhu & Wang 2019; Wang et al. 2019). Many works in fact process batches in a shuffled order without replacement and then incorrectly apply DP analysis for, e.g., Poisson sampling. Indeed, we use this same—incorrect—analysis for our DP-SGD baseline because it reproduces the previous state-of-the-art results for DP-SGD. The only requirement is that rather than shuffling the dataset for each epoch, the dataset is shuffled once and the same order of minibatches is used for each epoch. With this setup, $k$ epochs of training on a dataset of size $m$ with a batch size $B$ gives $n=mk/B$ total training steps, and satisfies $(k,m/B)$-participation. 2) We show that in important cases, e.g., Eq. 3, this participation schema allows for the efficient computation of sensitivity. 3) We will see in Section 3 that the more possible participation patterns $\pi$, the more constrained the problem of finding optimal mechanisms becomes (that is, fewer matrix factorizations satisfy sensitivty $\leq 1$). Hence, a relatively limited (but practical) schema like $(k,b)$-participation yields more favorable privacy-utility tradeoffs when we directly optimize matrices for this schema.

Gradient descent with fixed learning rate represents a canonical example of this setup: letting $\mathbf{x}$ represent the stream of unnoised gradients generated by a training procedure, the partially trained model at every step of the training procedure is simply a scalar multiple (with scalar being the negative learning rate) of the partial sums of this gradient stream. This partial-sum operation is represented by the matrix $\mathbf{A}$ of all 1s on the lower triangle. The ‘adaptivity’ of the stream $\mathbf{x}$ reflects the fact that the point at which we compute gradients during the training procedure depends on the previously released models, and is therefore required to capture privacy guarantees for ML model training.

We utilize the matrix mechanism (Li et al. 2015), which, provided a factorization $\mathbf{A}=\mathbf{B}\mathbf{C}$, computes the estimate

The sensitivity of the matrix factorization mechanism Eq. 1 is defined as

In the $d=1$ case, we have a much simpler polytope, $\mathfrak{D}^{1}_{\Pi}=\conv(\mathcal{D}^{1}_{\Pi})$ where

One might hope to show $\sens_{\mathfrak{D}^{d}_{\Pi}}(\mathbf{C})\leq\sens_{\mathfrak{D}^{1}_{\Pi}}(\mathbf{C})$, and the authors in fact initially conjectured this to be true. To our surprise, while this inequality holds under a variety of assumptions, it does not hold in general (Section H.2 gives a counterexample). We conjecture it is “almost” true; tightly bounding the necessary error term is an interesting open question. Empirically we have observed that for various query matrices $\mathbf{A}$ and $(k,b)$-participation with $d=1$, the optimal $\mathbf{C}$ satisfy (or almost satisfy) the condition $\mathbf{C}^{\top}\mathbf{C}\geq 0$ (element-wise non-negativity). In this case, we can show:

When per-step contributions bounded by $\zeta=1$, for any participation schema $\Pi$ and dimensionality $d\geq 1$, when $\mathbf{C}^{\top}\mathbf{C}\geq 0$ elementwise, we have $\sens_{\mathfrak{D}^{d}_{\Pi}}(\mathbf{C})=\sens_{\mathfrak{D}^{1}_{\Pi}}(\mathbf{C}).$

In particular, this implies that if $\mathbf{C}$ is optimal in the $d=1$ case and satisfies $\mathbf{C}^{\top}\mathbf{C}\geq 0$, it is also optimal in the $d>1$ case. This result is a corollary of H.1, which establishes additional conditions under which $\sens_{\mathfrak{D}^{d}_{\Pi}}(\mathbf{C})\leq\sens_{\mathfrak{D}^{1}_{\Pi}}(\mathbf{C})$ holds. All proofs are in Appendix H onwards.

Let $\mathbf{X}=\mathbf{C}^{\top}\mathbf{C}.$ When $\mathbf{X}$ has only nonnegative elements, one may reduce the problem of computing $\sens_{\mathfrak{D}^{d}_{\Pi}}(\mathbf{C})$ to

As an alternative to structural conditions on $\mathbf{C}$ or $\mathbf{X}$ allowing efficient exact computation of sensitivity for $d>1$, we can look to (reasonably tight) upper bounds on the sensitivity of $\mathbf{C}$. In the case of $(k,b)$-participation, one efficient method of computing upper bounds for the multiple-participation sensitivity of $\mathbf{C}$ has shown itself to be particularly useful:

In the $(k,b)$-participation case, $|\Pi|=b$. The complexity of computing the largest eigenvalue of the subselected $\mathbf{C}$ matrix is cubic in $k$. Thus, computing this upper bound is $\mathcal{O}(bk^{3})$, easily computable for the range of $k,b$ considered here ($k\leq 100,b\leq 500$).

Using our generalization of adaptive streams to multiple participations we obtain the following result (a straightforward generalization of Denisov et al. 2022). The proof is identical to (Denisov et al. 2022), except we replace the sensitivity bound with that for multiple participations obtained via 2.1.

Optimal Matrices for Multiple Epochs

We now present methods for computing optimal matrix mechanisms that are specialized to (optimized for) a specific participation schema $\Pi$ and query matrix $\mathbf{A}$. For example, Fig. 2 shows the optimal factorization for $\mathbf{A}$ representing SGD with momentum and learning-rate cooldown under $(k{=}6,b{=}342)$-participation, used in Section 5.3. Specializing the mechanism to both the participation pattern and specific query workload enables us to obtain state-of-the-art results in ML (Section 5).

We build on the approach of Denisov et al. 2022. We begin by defining the loss of interest, i.e., the total variance of noise added, for the mechanism defined in Eq. 1. Note that this loss characterizes other downstream tasks like DP mean estimation. Given $\mathfrak{D}$, assume that we may represent $\mathfrak{D}=\conv(\mathcal{D})$ for some finite set $\mathcal{D}$—as we have seen, this is the case, e.g., in $(k,b)$-participation. Then the loss which corresponds to total squared error of a factorization, at a fixed privacy level, may be expressed as:

Observing that $\forall\alpha$ the mechanism $\mathbf{A}=\big(\alpha\mathbf{B}\big)\big(\frac{1}{\alpha}\mathbf{C}\big)$ has identical loss, we conclude that we may consider the constrained version of the problem of minimizing this loss where $\sens_{\mathfrak{D}}(\mathbf{C})\leq 1$. Since for any $\mathbf{C}$, $\mathbf{B}=\mathbf{A}\mathbf{C}^{\dagger}$ produces the minimum-Frobenius norm $\mathbf{B}$-matrix, it is sufficient to solve:

With the change of variables $\mathbf{X}=\mathbf{C}^{\top}\mathbf{C}$, equivalently:

One of our main contributions is the following theorem which leads directly to efficient algorithms with provable optimality gaps for the mathematical program Eq. 6:

In the same setup as 3.1, the gradient of the dual function $g$ is: $\frac{\partial g}{\partial\mathbf{v}_{i}}=\mathbf{u}_{i}^{\top}\mathbf{U}^{-\frac{1}{2}}\big(\mathbf{U}^{\frac{1}{2}}\mathbf{A}^{\top}\mathbf{A}\mathbf{U}^{\frac{1}{2}}\big)^{\frac{1}{2}}\mathbf{U}^{-\frac{1}{2}}\mathbf{u}_{i}-1.$ Moreover, a maximizer of the dual ${\mathbf{v}^{\star}}$ must satisfy:

In the single-participation case of Denisov et al. 2022, $\mathbf{H}_{\mathbf{v}}=\diag(\mathbf{v})^{\frac{1}{2}}$, and Eq. 8 recovers the fixed point expression of that paper’s Theorem 3.2. Our 3.1 implies that the optimization methods presented in Denisov et al. 2022 may be applied, with suitable translation, to our setting; we use these methods to generate the optimal matrices studied empirically in Section 5.

While Eq. 6 gives a mathematical program for finding optimal matrices, two challenges arise: 1) for $(k,b)$-participation, we have $|\mathcal{D}^{1}_{\Pi}|=b2^{k}$, equal to the number of Lagrange multipliers that must be used to represent the constraint in Eq. 6. Hence, solving the dual problem will be intractible for moderately large $k$. 2) Even if we could surmount this issue, we cannot in general use 2.1 to bound sensitivity under vector contributions. In order to resolve these issues, in our experiments we introduce an additional constraint $\mathbf{X}\geq 0$ element-wise, requiring only $n^{2}$ additional Lagrange multipliers. This lets us only consider the positive vectors $\mathbf{u}\in\mathcal{D}^{1}_{\Pi}$ in the sensitivity constraint (following Eq. 3), and hence we need only a total of $n^{2}+b$ Lagrange multipliers in the dual optimization. Our empirical results indicate that $\mathbf{X}\geq 0$ holds for the optimal solution for the prefix-sum workload $\mathbf{A}$ (proving this conjecture is an interesting open problem). However, this conjecture does not hold in general, and in particular it fails for momentum workloads; see Section I.3. Nevertheless, enforcing $\mathbf{X}\geq 0$ produces mechanisms that lead to state-of-the-art learning performance in all cases. Section I.3 demonstrates these gaps on small examples and provides additional discussion.

FFT-based Matrix Factorization

Our work has two types of computation costs: optimization costs are those associated with optimizing and generating (or, computing) a mechanism whereas noise generation costs are those associated with using the mechanism to sample noise as part of a ML training algorithm. Once optimized, a single mechanism can be reused indefinitely to generate noise for other runs by simply resampling new noise and applying the same decoder. The best known methods for computing the optimal factorizations scale as at least $\mathcal{O}(n^{3})$ (Yuan et al. 2016; Denisov et al. 2022). This optimization cost can become intractable when $n$ grows too large. Thus, in this section we focus on reducing optimization computation at a small decrease in the achievable privacy-utility tradeoff.

A prime candidate for this goal is the Discrete Fourier Transform (DFT) because there are known algorithms both for nearly-optimal private convolutions (Fawaz et al. 2013) which are intimately related to the DFT, and for efficient calculation of the DFT using the Fast Fourier Transform (FFT) (Cooley & Tukey 1965; Nussbaumer 1981). We present an FFT-based mechanism that reduces noise generation costs, prove rigorous DP guarantees for it, and show that these lead to near-optimal privacy-utility tradeoffs in the single-epoch setting. We provide two improvements over prior work (Fawaz et al. 2013): 1) extending the result to the multi-epoch and multi-dimensional setting and 2) providing explicit non-asymptotic analysis of the algorithm’s utility.

J.1 of Gray 2006 (restated in Section J.1) shows there exists a diagonal $\Sigma$ such that $\mathbf{A}_{\sf circ}=\mathbf{F}^{*}\Sigma\mathbf{F}$ for diagonal $\Sigma$, where $\mathbf{F}$ is the DFT matrix. Then, $\mathbf{A}_{\sf circ}$ can then be factorized as $\mathbf{A}_{\sf circ}=\mathbf{B}_{\sf circ}\mathbf{C}_{\sf circ}$ where $\mathbf{B}_{\sf circ}=\mathbf{F}^{*}\Sigma^{1/2}$ and $\mathbf{C}_{\sf circ}=\Sigma^{1/2}\mathbf{F}$.

The (complex-valued) matrix mechanism specified by the factorization above (and presented as Algorithm 1 of Appendix C) is empirically nearly optimal in the class of matrix-factorization-based mechanisms, as we show in Appendix J. Though we prove a simple zCDP guarantee for any participation schema having at most $k=\max_{\pi\in\Pi}|\pi|$ participations in 4.1, we can instead use our 2.1 when the participation schema is known in advance (as in our experiments with $(k,b)$-participation).

Under $k$-participation, Algorithm 1 satisfies $(k^{2}\rho)$-zCDP.

Observe from 2.1 and 2.2 that the privacy guarantee of our multi-participation adaptive setting is independent of the choice of decoder, $\mathbf{B}$. Thus, instead of taking $\mathbf{B}_{\sf circ}$ above, we take the optimal decoder using the Moore-Penrose pseudoinverse of a distributionally equivalent real-valued encoder, as discussed in Appendix K. This optimization leads to significant improvements in the privacy-utility tradeoff (see Fig. 9 in Appendix E).

Though we define the optimal FFT decoder as a pseudoinverse, observe that we do not need to optimize (or even compute) the decoder; by Appendix K, the problem is reduced to that of solving a highly structured linear system. However, we find that even suboptimal implementations using the pseudoinverse can still factorize a mechanism for $n=10,000$ in 146 minutes on a V100 GPU, remaining well within practical requirements since we need only generate a mechanism once, before it can be reused indefinitely for training. In contrast, computing optimal matrices becomes practically difficult near $n\approx 10,000$, taking $24$ hours to compute an effective factorization for $n=8,192$ using batch-priority cloud CPU resources. We remark that this regime of $n$ is highly practical, e.g., standard federated benchmarks use $n\approx 2000$ (Reddi et al. 2020) and our central image classification uses $n=2000$. In terms of noise generation, the FFT mechanism shows preferable asymptotic properties, scaling as $\mathcal{O}(dn\log^{2}{n})$. However, even the optimal matrix mechanism with runtime scaling as $\mathcal{O}(dn^{2})$, noise generation on a GPU (even with significantly suboptimal implementation) takes negligible time. Further, noise from our mechanisms can be pre-generated if needed, at a tradeoff in the (typically cheaper) disk space. We discuss these tradeoffs in Appendix C.

Empirical Evaluation

We compare four main mechanism classes: tree-based mechanisms (Honaker 2015), including ‘tree-completion’ of Kairouz et al. 2021; our FFT mechanism; our optimal factorizations; and DP-SGD (incorrectly) assuming amplification via Poisson subsampling (Abadi et al. 2016b).

The manner in which baselines from the extant literature map to this setting can be found in Section D.1. Since the matrix mechanism reduces privacy cost of training to that of the release of a single Gaussian mechanism, accounting in our case becomes quite simple; see Section D.2.

Our work can be viewed as a drop-in replacement for noise sampling in DP optimizers through one additional step in training. Concretely, the three main steps to create DP-SGD are to 1) compute the per-example gradients, 2) clip each one to some chosen threshold $\zeta$, then compute the average as $\mathbf{x}$, and 3) add noise $\mathbf{z}\sim\mathcal{N}(0,\sigma)$ to $\mathbf{x}$ where $\sigma,\zeta$ are calibrated for DP. In our work, we define an additional step 4) which uses Eq. 1 to generate the noise as $\mathbf{B}\mathbf{z}$ using $\mathbf{x}$ from step 2) and $\mathbf{z}$ from step 3) with $\sigma=1$. Because $\mathbf{A}=\mathbf{B}\mathbf{C}$ is the chosen optimizer workload (e.g., residual prefix-sum or momentum-sum corresponding with SGD and SGD-M), we may ignore operations on $\mathbf{x}$.

To generate the $\mathbf{B},\mathbf{C}$, we must solve Eq. 6. Equivalently we must use Algorithm 1 of Appendix C for the FFT mechanism or Appendix K for the Optimal FFT Decoder. Observe that solving this linear task minimizes the noise required for a formal DP guarantee; instead, the ML optimizer generates gradients $\mathbf{x}$ for the non-linear learning task. We can instantiate the sensitivity constraint set $\mathcal{D}$ to encode the exact $(k,b)$-participation schema used in training. However, an encoder $\mathbf{C}$ factorized for one value of $k$ can be applied for another, by simply computing its new sensitivity (due to our Section 2), though this will alter its privacy-utility tradeoffs. For example, our MF1, 6e in Fig. 3 uses this to extend Denisov et al. 2022 to $k>1$.

When applying a mechanism that is determined independently of the number of participations (e.g., FFT or tree aggregation) or extending an optimal mechanism for $k$ participations to a larger number of participations, the sensitivity may scale poorly in $k$. In such cases, it may actually have lower sensitivity to reuse a single encoder multiple times over the course of $n$ steps, and hence more favorable privacy-utility tradeoffs. We term this approach encoder stamping. This also provides a straightforward method for extending any factorization to handle more iterations without, e.g., re-optimizing Eq. 6. Combined with our Section 2, stamping lets us apply mechanisms from Denisov et al. 2022, e.g., MF($\mathbf{k{=}1,n{=}1000}$)$\mathbf{\times 2}$ in Fig. 1; mechanisms with stamping have “$\times s$” appended in this way. Discussion of stamping and its relation to existing literature are in Section D.4.

2 Example-level DP for an Image Classification Task

We train image classification models on CIFAR10 (Krizhevsky 2009) which has become a de facto standard for comparing DP ML algorithms—sufficiently easy for existing DP algorithms to achieve nontrivial accuracy, but not so simple so as to be unable differentiating approaches. Details on our full setup are in Appendix E; generally, they match those of Kairouz et al. 2021. Notably, we make improvements on their Online Honaker-based approach by not just completing the tree with virtual steps, but also zeroing out noise from virtual steps as detailed in Section D.3. We find this led to significant improvements around a few percentage points. For all matrix mechanisms except the Denisov et al. 2022 baseline and our Optimal MF$(k=20,n=2000)\times 1$, we optimize over the stamps $s$ by the losses in Section D.5, which we find well match the ordering in ML accuracy.

In contrast to Section 5.3, in this section we only compare factorizations of the prefix-sum matrix; we do not incorporate momentum or cooldown directly into the mechanisms, though we use both momentum and cooldown as postprocessing for the matrix-factorization-based mechanisms and report results for the best settings we find (with both). For DP-SGD, we report results for the best setting (no momentum, with cooldown). Details are in Appendix E.

First, we see that the optimal factorization for the target $(k,b)=(20,100)$ setting outperforms all other mechanisms across (nearly) all privacy levels, only slightly underperforming DP-SGD with amplification at $\varepsilon=2$. To the best of our knowledge, this represents the first empirical demonstration of an ML algorithm which is competitive with DP-SGD into this high-privacy regime, without any amplification by sampling. The FFT (optimal decoder) mechanism outperforms all baselines, again well toward the high-privacy regime at $\varepsilon\approx 4$. Though at a worse privacy-utility tradeoff compared with our multi-epoch optimal matrices, this mechanism shows promise for outperforming prior work when $n$ grows too large for generating optimal factorizations.

A major reason we outperform DP-SGD is because we take a fundamentally different approach that enables us to account and optimize for the specifics of the entire training procedure in a single shot. That is, DP-SGD views itself as the repeated independent application of a single mechanism on each iteration, using techniques like strong composition to obtain a finite DP guarantee across independent and arbitrary queries (gradient steps). Our matrix-factorization based ML training procedures, by contrast, consider the entire training procedure to be the application of a single mechanism, parameterized by factorizations of the matrix $\mathbf{A}$; we optimize over this space of mechanisms. We believe this is a more powerful class of mechanisms because it essentially enables us to (anti-) correlate noise across iterations to achieve minimal total squared error (DP-SGD can only use independent noise).

3 User-level DP for a Next Word Prediction Task

User-level DP for language models is an important real-word task (McMahan & Thakurta 2022). There is much history in the DP language modelling literature which we briefly describe in Section F.1. We use the standard benchmark: StackOverflow next-word prediction (Reddi et al. 2020). We use the same empirical setup as Kairouz et al. 2021 and Denisov et al. 2022 except notable changes below. Details are in Section F.2.

We conducted initial simulations which verified that two observations from Denisov et al. 2022 for the single-epoch setting extended to our multi-epoch and large-batch (1000 clients/round instead of 167) setting. First, linear server learning rate cooldown from $1\times$ to $0.05\times$ over the final 512 rounds offered a small improvement over constant rates (more so in higher-privacy regimes). Second, optimizing and factorizing with momentum and cooldown encoded in the query matrix $\mathbf{A}$, rather than applying both as post-processing, consistently offers a small benefit. See Section F.2 for details. Thus, we fix these preferable design choices here and compare the following algorithms.

All algorithms train for $6$ epochs $2052$ rounds, and a large-batch 1000 clients/round (better for DP training) unless otherwise noted. Honaker, 6e is the DP-FTRL algorithm of Kairouz et al. 2021, trained for 2048 rounds (a power of 2). MF1, 1e (Denisov et al. 2022) is the state-of-the-art for single-epoch training, using $167$ clients/round and $k{=}1$. MF1, 6e uses our Eq. 3 to take the (non-negative) $(k{=}1)$-optimized matrix of the previous approach but compute the sensitivity under $(k{=}6,b{=}342)$-participation, allowing us to train for 6 epochs (2048 rounds) with large batches. MF6, 6e is our approach directly optimizing the matrix factorization for $(k{=}6,b{=}342)$-participation via Eq. 6. DP-SGDM, 6e is the DP-FedAvg algorithm of McMahan et al. 2018, to 2052 rounds, and accounted with Poisson sampling—incorrect, though standard, as noted in Section 1. DP-GDM, 2052e estimates the infeasible-to-actually-run benchmark of full-batch gradient descent for 2052 rounds and 2052 epochs, or $342\times$ the computation cost of our 6 epoch runs. We compute the exact privacy cost, and estimate the accuracy from experiments with 1000 clients per round, following the methodology of Kairouz et al. 2021. This is essentially an upper bound on the best privacy-accuracy tradeoffs with unlimited computational resources.

We find that our MF6, 6e, is the best feasible private result at $25.25\%$ accuracy and $(17.7,10^{-6})$-DP. This exceeds the non-private baseline of $25.2\%$ accuracy reported by Kairouz et al. 2021, is within the margins of small hyperparameter tuning differences of our improved non-private baselines (25.43%) and private full-batch gradient descent ($25.31\%$). At $(8.84,10^{-6})$-DP, MF6, 6e achieves $24.94\%$ accuracy, substantially improving over the previous state-of-the-art at this privacy level given by MF1, 1e at $24.11\%$ accuracy, and considerably improves on both the accuracy and privacy of Honaker, 6e ($24.86\%$ accuracy at $(21.0,10^{-6})$-DP). In fact, we achieve better accuracy at $\varepsilon=8.84$ than prior methods achieve at $\varepsilon=17.7$. DP-SGDM, 6e is outperformed by MF6, 6e and MF1, 6e across all $\varepsilon$ values evaluated. Fig. 12 in Appendix F shows results for each learning rate separately, with numeric results in Tables 5 and 4.

Discussion and Conclusions

Our work significantly improves the privacy-utility tradeoffs in DP ML. Indeed, our work outperforms the state-of-the-art (and, DPSGD with amplification) by $\approx 5$ percentage points across many privacy levels—and as low as $\varepsilon\approx 2$—with practically implementable assumptions.

We compare our mechanisms with DP-SGD on a level-ground using well-performing but not state-of-the-art models and training protocols—e.g., very large models, augmentations before clipping (De et al. 2022), public data usage, and large batch sizes can all aid training. Many, if not all, of these techniques are applicable to our setting and can thus be used with our mechanisms to realize additional absolute performance gains, again, likely beyond the performances achieved by DP-SGD. Appendix G contains limitations and ethical considerations.

The authors would like to thank Adam Smith for helpful conversations in the formulation of the FFT approach. Shuang Song helped ensure correctness and comparability with Kairouz et al. 2021 on CIFAR. Galen Andrew and Sewoong Oh contributed their expertise in matrix calculus, and Ryan McKenna provided helpful feedback on the manuscript as well as contributing the concise counter-example presented in Section H.2.

References

Appendix A Summary of notation and terminology

The following table summarizes the notation used throughout the paper:

We utilize terminology from federated learning as well as standard centralized training, which generally map as follows:

Appendix B Generalized sensitivity as an operator norm

Eq. 2 shows that our generalized notion of sensitivity can be viewed directly as a particular operator norm. To see this, view $\mathbf{C}:\mathbf{V}_{1}\rightarrow\mathbf{V}_{2}$ as a linear operator from vector space $\mathbf{V}_{1}$ to $\mathbf{V}_{2}$. Then with $\|\cdot\|_{(1)}$ the vector norm on $\mathbf{V}_{1}$ and similarly for $\mathbf{V}_{2}$, an operator norm is defined as

the vector norm induced by $\mathfrak{D}$ (the fact that $\mathfrak{D}$ is a closed, convex, symmetric set ensures this is a norm). Note $\mathbf{u}\in\mathfrak{D}\ \Leftrightarrow\ \|\mathbf{u}\|_{\mathfrak{D}}\leq 1$. Thus, we have

Appendix C The FFT Mechanisms and Reducing Computation

We propose two FFT mechanisms. First, we propose the FFT mechanism which is described in Algorithm 1. This mechanism has the same computation complexity—no optimization costs and $\mathcal{O}(nd\log{n})$ noise generation—as the Honaker method used in Kairouz et al. 2021 but at a better privacy-utility tradeoff, as shown in Fig. 9. The privacy and utility analysis can be found in Appendix J.

The FFT Optimal Decoder (FFT_Opt_Dec) mechanism presented in Section 4 represents taking (a real-valued translation of) the encoder $\mathbf{C}$ from Algorithm 1 and using the optimal decoder, defined in terms of the Moore-Penrose pseudoinverse of $\mathbf{C}$. Similarly to Algorithm 1, there is no need to construct a literal matrix to multiply by in the case of noise defined by the optimal decoder; noise generation time of the mechanism, however, increases by a logarithmic factor to $\mathcal{O}(nd\log^{2}{n})$ (as discussed in Appendix K). This complexity is still feasible for many steps (which we will discuss below) and comes with significant utility benefits (see Fig. 1).

All the mechanisms we study scale as either $\mathcal{O}(n^{2})$ or $\mathcal{O}(n\cdot\polylog(n))$. For our $n=2000$ step environments, and even far beyond to $n\approx 10,000$, our algorithms can be efficiently realized on GPUs with runtime on the order of seconds per step (including computing and applying gradients and noise). The main challenge in these cases are storing the $n^{2}+nd$ coordinates in GPU memory (the former for the decoder matrix, the latter for the noise samples). Given that each of the $d$ coordinates of noise can be sampled independently, this algorithm is straightforwardly parallelizable and so work may be partitioned across many processors when needed. Noises could instead be pre-generated in an entirely separate process, and stored on disk, to be loaded into memory row-by-row concurrently with training.

Appendix D Mechanisms under consideration: baselines, subtleties, and losses.

Kairouz et al. 2021 and Denisov et al. 2022 both present approaches for ML model training which can be understood as instances of the matrix mechanism–the former grounded in the binary-tree mechanism as refined by Honaker 2015, and the latter explicitly optimizing a factorization under single participation. These two works yield two natural baselines:

Kairouz et al. 2021 explores ‘tree restarts’ (generalized as our notion of ‘stamps’, $s$, in Section 5) and the so-called ‘tree-completion trick’ for the Honaker estimator-from-below variant of the binary tree method for computing differentially private prefix sums; the matrix-factorization perspective on this estimator allows us to implement slightly optimized versions of these methods; see Sections D.3 and D.4.

Denisov et al. 2022 computes optimal factorizations of various optimization-related matrices, though only for a single epoch. For these matrices, we leverage the results in Section 2 to directly compute the sensitivity of the encoder matrices for multiple participations.

These two papers can be combined in other ways as well; e.g., the results of Denisov et al. 2022 show that the ‘fully efficient estimator’ of Honaker 2015 may be used as a drop-in replacement for the estimator from below in Kairouz et al. 2021. We focus on the two mechanisms specified above as the natural baselines for the present work.

D.2 Privacy Accounting

Privacy costs are computed as a single application of the Gaussian mechanism to $\mathbf{G}\mathbf{C}$ using the PLDAccountant provided by the Google DP Library https://github.com/google/differential-privacy. We also use this accountant to analyze DP-(S)GD baselines (which require more complex accounting), yielding a small improvements in $\varepsilon$ over the Renyi-DP accounting used in prior works.

D.3 Improvements to ‘Tree completion’ By Removing Noise from Virtual Steps

The “tree completion” trick of Kairouz et al. 2021 is used on the last step of any restart (in ours, stamp) to reduce the noise added on this step. This is achieved by adding virtual steps (with 0 inputs) until the final step of that level in the tree, because this noise will be the lowest in that level. In this section, we show how to further improve on this trick and that our matrix mechanisms make analyzing such tricks easier. Our implementations of the binary-tree baslines Honaker 2015; Kairouz et al. 2021 utilize these improvements.

For the online Honaker estimate, Honaker 2015 obtain a DP estimate for the release node $i\in[n]$ (representing the prefix sum until $i$) but summing the corresponding subtrees prior to this node. These are exactly the subtrees corresponding to the binary representation of this node (Honaker 2015). Then, the variance required to release node $i$, with subtrees of height $0,1,...,l_{i}-1$, is

Notice that reaching a new height in the tree decreases the variance needed. In Kairouz et al. 2021 and just before terminating on some non-power-of-two step $n^{\prime}<n$, they run $n-n^{\prime}$ virtual steps on zero gradients. This enables their mechanism to use the minimal noise for the power-of-two-step for the final real step.

However, notice that in both these cases, the methods assume that these virtual steps must be privatized. Indeed, they do not need to be because we know apriori that these steps are virtual, i.e., not corresponding to real gradients. Thus, in our methods we account for this in our mechanism and reduce the noise of the final step accordingly. Importantly, this can be computed without altering the asymptotic runtime and storage complexity of the algorithm: on the last step, the contributions of the virtual steps to the power-of-two noise can be calculated using, e.g., a second binary tree, and removed. This leads to a significant benefit in the loss as we observed in Table 1 in Section D.5.

We believe this oversight of prior works showcases the power of our matrix mechanism approach. Indeed, let $\mathbf{C}_{\text{tree}}$ be a matrix representing the (complete) binary tree used in the mechanisms of Honaker 2015; Kairouz et al. 2021 with $2^{\lceil\log_{2}(n)\rceil}$ leaves; for the sake of concreteness, assume this is the matrix constructed in Appendix C of Denisov et al. 2022. Let $\mathbf{B}_{\text{Hon}}$ represent the Honaker estimator-from-below; in our language, the decoder used by (Kairouz et al. 2021).

The tree completion trick of (Kairouz et al. 2016) can be understood as follows. The matrix $\mathbf{B}_{\text{Hon}}\mathbf{C}_{\text{tree}}$ is of size $2^{\lceil\log_{2}(n)\rceil}\times 2^{\lceil\log_{2}(n)\rceil}$. In the case that $n$ is not a power of two, the penultimate rows and columns of this matrix will go unused. However, for this factorization, the variance added by $\mathbf{B}_{\text{Hon}}$on the final row will be quite small, due to the binary tree’s redundancy in encoding estimates of this sum. The matrix we wish to factorize is a prefix-sum matrix $\mathbf{S}$ of size $n\times n$; this matrix can be expressed as any one of a family of transformations of the (potentially larger) product $\mathbf{B}_{\text{Hon}}\mathbf{C}_{\text{tree}}$:

One more optimization becomes clear when tree completion is formulated in this manner. Similar to our optimal decoder of Section 4, any decoder can be used without changing the sensitivity of the encoder. Noting that nonzero entries in the decoder increase our loss of Eq. 4, we can simply zero out the columns of the decoder corresponding to these virtual steps—this decreases the loss, preserves the same error in the DP estimate of the prefix sum (the inputs are 0), and maintains the same DP guarantee. In other words, we need not account for the noise, or the error it introduces, of virtual steps. We now provide a more rigorous explanation.

The image of $\mathbf{C}_{\text{tree}}\mathbf{E}$ can be contained in an axis-aligned subspace; effectively, the subspace corresponding to elements that may be nonzero in the binary tree when run for $n$ steps. In other words, the columns of the decoder corresponding to rows of the encoder that are removed via the projection need not be included: because the input is not processed. Therefore, denoting the projection onto this subspace by $\Pi$, we may write:

further reducing the variance of the decoder (without increasing the sensitivity of the encoder) in this incomplete binary tree case.

In our implementations of the online Honaker mechanism, we freely use these tricks, in addition to exact calculations of the sensitivity of the mechanism enabled by noting that the encoder is all-nonnegative and the observations of Section 2, leading to some improvements in these mechanisms over those in the existing literature. No changes to accounting are required, as privacy is inherited from the matrix mechanism perspective.

D.4 Stamping: Repeated mechanisms in the matrix-factorization setting

In stamping, we define a new encoder matrix as the Kronecker product of some given encoder $\mathbf{C}$ with an $s\times s$ identity matrix $\mathbf{I}$, creating a new $sn$-dimensional linear DP query mechanism. Assuming $\mathbf{C}$ is of shape $n\times n$, the resulting encoder is an $sn\times sn$ block-diagonal encoder matrix, formed by ‘repeating’ the matrix $\mathbf{C}$ along the diagonal.

Kairouz et al. 2021 explored ‘restarting’ their binary-tree based prefix sum estimation mechanism, treating the number of restarts used as a hyperparameter, and treating the result of a ‘completed’ application of the binary tree as fixed. For linear operators $\mathbf{A}$ with constant columns below the diagonal, this approach may be used to construct a new factorization from an existing one. This constant-columns property, or a block-based variant thereof, is required to simply treat the output from a ‘completed’ application of the existing mechanism as fixed; for a general matrix $\mathbf{A}$, there is no clear prefix property which can be leveraged for this purpose.

Taking the matrix view, one may construct a similar encoder/decoder pair for the prefix-sum matrix by reusing the initial decoder $\mathbf{B}$ on the block-diagonal and fixing the columns to simply repeat the final row of this decoder below the block-diagonal; notice that the constant-column property of the prefix sum matrix guarantees that this construction appropriately factorizes $\mathbf{A}$. The noise that the matrix mechanism thus constructed adds can be implemented as a ‘restarted’ tree mechanism; however, since we compute sensitivities exactly for decoders of this structure as in Section 2, the privacy properties of these mechanisms we construct to replicate the ‘restarts’ of (Kairouz et al. 2021) may not be identical to those presented there, where accounting is performed by composition.

The matrix-mechanism perspective additionally allows one to apply the ‘stamping’ construction to any linear operator $\mathbf{A}$ (e.g. the momentum matrix), where a reuse of fixed previous outputs is not possible. A ‘stamped’ factorization of any $\mathbf{A}$ may be obtained, for example, by matrix pseudoinversion: letting $\mathbf{B}=\mathbf{A}\left(\mathbf{C}\otimes\mathbf{I}\right)^{\dagger}$. This defines a legitimate factorization of any $\mathbf{A}$ requiring only suitable non-degeneracy assumptions on $\mathbf{C}$, and indeed represents the optimal decoder for the stamped encoder.

The pseudoinverse-based construction can, even in the prefix-sum case, be quite different from a construction designed to replicate ‘restarted’ mechanisms. For example, if $\mathbf{C}$ is a matrix representation of the binary-tree encoder and $\mathbf{I}=1$, then the resulting decoder matrix represents the full, rather than online, Honaker decoder (Honaker 2015); the validity of this mechanism in the adaptive streaming setting was only shown quite recently by Denisov et al. 2022.

For comparability with existing literature (and to preserve potential for an efficient implementation), however, all of the tree-based mechanisms we explore in the main body have decoders which replicate the setting of composition We show how the optimal binary-tree decoder differs from the online, composition-based decoder in Fig. 6. All other stamped mechanisms used the optimal decoder.

Considering instantiation enables us to directly analyze and minimize (over $s$) the stamped mechanism’s multi-epoch loss Eq. 4 without running compute intensive ML experiments. Indeed, we observe in Table 1 of Section D.5 that for mechanisms which were not explicitly optimized for the $(k,b)$participation setting, there exist stamped mechanisms ($s>1$) with much lower loss that correspondingly led to much more performant ML models (e.g., Figures 7 and 8 of Appendix D).

Interestingly, with our capacity to measure mechanisms at a single shot in the multi-epoch setting, we see a similar trend as was observed for restarts in Kairouz et al. 2021: that ‘stamped’ mechanisms have lower total loss than their non-stamped full-tree counterparts in the 20-epoch, 100 steps / epoch setting (see Table 1); training performance of these ‘stamped’ mechanisms on CIFAR10 can be found in Fig. 7. However, there is a significant improvement in our approach in that we can now directly tune this hyperparameter without the need to actually run ML training. This lets us reduce computation by only analyzing the loss of the generated matrices and then running the mechanism with the lowest loss.

D.5 Factorization losses and per-iterate variance

As a first measure of the privacy-utility tradeoff, we compare the losses of each mechanism from Eq. 4 for factorizations of the matrices under consideration.

We compare measured losses of several factorizations of the prefix-sum matrix for the $(k,b)=(20,100)$ setting of the CIFAR experiments in Table 1. We plot the per-iterate variance of the mechanisms in Fig. 1, along with several variations, in Fig. 4. In Fig. 5, we plot the per-iterate variance distribution of factorizations of the momentum and cooldown matrix (described in Section 5) at a fixed variance level, and with various privacies, computed for the $(k,b)=(6,342)$ participation setting.

Figs. 4 and 5 both demonstrate the effect of $(k,b)$-participations on the optimization problem. Particularly interesting to consider are the optimally-factorized matrices; in both cases, the epoch structure is clearly visible in the manner in which the mechanisms distribute variance. We see also the effect of ‘stamps’ $s$ in the variance distribution, effectively a proxy for the epoch structure directly accounted for by the optimal mechanisms.

Appendix E Details and additional experiments for CIFAR10.

We train image-classification models using the CIFAR10 dataset as hosted in tensorflow-datasets, containing 50,000 training and 10,000 test examples. We evaluate and compute test accuracies on the entire test set, following the open-sourced code of Kairouz et al. 2021. We reuse the network architecture, dataset processing and initialization strategies presented in Kairouz et al. 2021; in particular, the architecture we use can be found in their Table 2 (b).

We train all mechanisms for 20 epochs with batch size of 500, yielding 100 steps per epoch and 2000 total. After performing some small initial grid searches, we settled on using linear learning rate cooldown to $0.05\times$ the initial learning rate over the last 500 steps of training. We found this consistently improved utility for all mechanisms and privacy levels.

As mentioned in Section 5, for this 20-epoch training setup, we only compare factorizations of the prefix-sum matrix, and do not include any factorizations of matrices which incorporate momentum of learning rate cooldown directly in the mechanism itself (Denisov et al. 2022). We sweep over learning rates of values $(1\times 10^{i},2\times 10^{i},5\times 10^{i})$ for $i$ in $\{-2,-1\}$; for all mechanisms and noise levels, optimal values were in the interior of this sweep. We sweep over momentum values of $0,0.85,0.9,0.95$ though find nonzero momentum works best for all matrix mechanisms, and no momentum works best for DP-SGD at our scale as found previously by Kairouz et al. 2021.

For Honaker and FFT-based factorizations, there is no known a-priori way to choose the optimal number of $s$ for a given $(k,b)$ setting. Therefore we treat the value $s$ as a hyperparameter, and sweep across it, for $s\in\{1,2,5,10,20\}$. As can be seen in Table 1, the optimal $s$ for both of these factorizations was in the interior of this sweep. As shown, e.g., in Fig. 7, the training-time performance of these mechanisms matched the expected order for computed loss. This value $s$ represents an extra hyperarameter which must be set for the Honaker and FFT mechanisms; to the best of our knowledge, computing the loss for various instantiations of these mechanisms via Eq. 4 represents the only known method for setting this parameter other than simply training models.

We also apply our sensitivity analysis of Section 2 to the matrices of Denisov et al. 2022 which are optimized for $k=1$. In doing so, we can also optimize the number of stamps which we do. We report the best results as identified by the losses in Table 1.

Appendix F Additional StackOverflow Details

Language models trained on user data are an important real-world application of DP training, as these models can memorize their training data if appropriate mitigations are not applied (Carlini et al. 2019; Song & Shmatikov 2019; Carlini et al. 2021). Since one user might contribute 1000s of tokens (training examples) to a dataset, it is particularly important to consider user-level guarantees (McMahan et al. 2018). Building on the approach of Kairouz et al. 2021, Google recently announced the first-ever launch of a language model trained on user data with a formal user-level DP guarantee ($\rho=0.81$ zCDP), further demonstrating the importance of this application (McMahan & Thakurta 2022).

The StackOverflow next-word prediction task, introduced in (Reddi et al. 2020), has become a benchmark problem for DP training, and our experimental setup here fixes the same model and adapts hyperparmaeters from previous work including Kairouz et al. 2021; Denisov et al. 2022.

F.2 Hyperparameter tuning and initial experiments

All runs use server momentum 0.95, a client learning rate of 1.0, and a server learning-rate cooldown schedule for the last 25% of rounds. The clipping norm was fixed at $\zeta=1$. Zeroing outlier updates and using 1000 clients/round (6 epoch runs) allows the use of the higher server learning rates. MF1, 1e replicates the result of single-epoch training from Denisov et al. 2022; note that with 167 clients/round and this mechanism, the higher learning rate does not appear to help. Fig. 10 gives preliminary experimental results which informed the main experiments used in the paper. Note that the y-axis range (Test set accuracy) is highly compressed, and so the primary point of comparison is on epsilons. For example, Denisov et al. 2022 shows that cross-run variation of 0.002 or more is typical.

The 6 horizontal lines give test-set accuracy for various non-private training mechanisms. The “Unnoised MF” runs correspond to the same code path used for privacy, but without any noise addition. In particular, these use momentum with learning rate cooldown; the other unnoised runs use a standard FL implementation with momentum but a fixed learning rate schedule; “cpr=167“ corresponds to one epoch of training (167 clients/round), and “cpr=50” is 50 clients/rounds (only about 1/3 of an epoch). This last non-private baseline uses the best hyperparameters for FedAvgM from Reddi et al. 2020.

The two “Unnoised MF” runs with accuracies between 0.246 and 0.248 are functionally identical, and the line near 0.248 accuracy is the same except it does not use learning-rate cooldown. Thus, for the given learning rates, we see the higher-epsilon private runs are adding sufficiently small noise that the accuracy is essentially equivalent to unnoised baselines with the same hyperparameters. However, using larger learning rates can achieve accuracy over $25\%$, even with privacy as in the case of the MF-6-6 run, hence motivating the inclusion of larger learning rates in the main experiments.

The MF (Matrix Factorization) runs with “prefix” in the name correspond to computing an optimal factorization of the prefix-sum matrix (lower triangular matrix of ones) and then applying momentum (and possibly learning rate cooldown) as post-processing. The other MF runs directly factor the momentum or momentum+cooldown matrix.

F.3 Impact of zeroing-out large-norm updates

F.4 Complete results

In this section we give additional details on our main grid of experiments. Fig. 12 uses the same data as Fig. 3, but shows results for each learning rate individually. Tables 4 and 5 give the mean, minimum, maximum, and standard deviation of test-set accuracy corresponding to Fig. 12, as well as the number of replicated experiments (‘count’).

Table 3 gives the noise multipliers to achieve our various privacy targets $\varepsilon$. Due to a change in the accountant used, we have slightly different $\varepsilon$ targets around 8.8 for the different methods. Note the noise multipliers here are incomparable between the MF and (S)GD mechanisms in terms of the total noise introduced. For matrix factorization, we sample $\mathbf{Z}\sim\mathcal{N}(0,\zeta^{2}z^{2})$ for noise multiplier $z$; this noise is applied after mapping the raw gradients/updates $\mathbf{x}$ through the linear map $\mathbf{C}$ which is normalized so the total sensitivity is $\zeta$. For the (S)GD mechanisms, we add noise $\mathcal{N}(0,\zeta^{2}z^{2})$ independently to each model update. In all cases, noise is added to the sum of per-user updates, so the effective noise in the average update scales down with the number of clients per round.

F.5 Optimal matrix mechanisms

Appendix G Limitations and Ethical Considerations

The major limitation of our approach is the computation required to generate the optimal matrices. Though our optimal FFT decoder bridges the gap between the mechanisms without optimizer costs and our optimal mechanism, it still leaves some room for improvement in the privacy-utility tradeoff. We believe this is an important area for future work.

Our work aims to enable privacy-preserving ML with rigorous DP guarantees in broader settings via improved (state-of-the-art) privacy-utility tradeoffs. Though DP is the gold-standard in this area, it is currently impossible to train high-utility ML models with tight DP guarantees $<1$ in many, if not all, real-world settings without large amounts of public data. In this regime, it is not well-understood what privacy leakage may occur.

Appendix H Analysis for Section 2

All the entries of $\mathbf{C}^{\top}\mathbf{C}$ are non-negative.

The rows of $\mathbf{G}$ are all co-linear, $\mathbf{G}_{[i,:]}=u_{i}\cdot\mathbf{G}_{[1,:]}$ for $u_{i}\in\{-1,1\}$, $i>1$.

The rows of $\mathbf{G}$ are all orthogonal, $\langle\mathbf{G}_{[i,:]},\mathbf{G}_{[j,:]}\rangle=0$, $\forall i\neq j$, $i,j\in[k]$.

Furthermore, the following statements are also true without assuming conditions (1)-(3) above.

$\left\|\mathbf{C}\mathbf{G}\right\|_{F}\leq\frac{\pi}{2}$.

If we replace the condition on $\mathbf{G}$ to $\forall i\in[k],\|\mathbf{G}_{[i,:]}\|_{1}\leq 1$, then $\left\|\mathbf{C}\mathbf{G}\right\|_{F}\leq 1$.

In the following we prove each of the individual cases of Theorem H.1.

An equivalent argument is used in Denisov et al. 2022.

where Eq. 12 follows from the standard $\texttt{A.M.}\geq\texttt{G.M.}$ inequality.

We give a construction for a $\mathbf{u}$ that shows this. Observe

Observe we can choose $u_{1}=1$ and then $u_{i}=\text{sign}(b_{i})$ since $b_{i}$ depends only on $\mathbf{C}$ and the previously fixed $u_{j}$ for $j<i$, ensuring the double sum on the right is non-negative, completing the proof.

This argument is essentially due to (Denisov 2023); we inline here for completeness.

which is to say that $\|\mathbf{C}\|_{\infty,2}\leq 1$. Duality, therefore, implies that $\|\mathbf{C}^{\top}\|_{2,1}\leq 1$, and $\|\mathbf{X}\|_{\infty,1}\leq 1$.

Now, for $\mathbf{X}=\left(x_{ij}\right),$

where $K$ (known as the Grothendieck constant) is independent of $d$, but not known exactly in general.

However, in our case, our matrix $\mathbf{X}$ is symmetric by construction; in the case of symmetric $\mathbf{X}$ the constant $K$ may be replaced by $\pi/2$ (see Eq. 43 of (Khot & Naor 2011)), and this constant is known to be sharp (IE, the best constant which holds for all symmetric matrices $\mathbf{X}$ and in all dimensions $d$).

First notice that since $\left\|\mathbf{C}\mathbf{G}\right\|_{F}^{2}$ is a convex function, the maximum happens at the extreme points of the constraint set $\mathcal{G}=\{\mathbf{G}\mid\forall i\in[k],\|\mathbf{G}_{[i,:]}\|_{1}=1\}$. We use H.1 to identify the extreme points of $\mathcal{G}$.

The set of extreme points of the set of $k\times d$ row-stochastic matrices are precisely the set of row permutation matrices, i.e., set of the matrices with entries in $\{0,1\}^{k\times d}$ and each row has exactly one non-zero entry.

Notice that the constraint on any matrix $\mathbf{G}\in\mathcal{G}$ is oblivious to the sign, meaning, we can flip the sign of any set of entries in $\mathbf{G}$ and the new matrix will still be in $\mathcal{G}$. This along with H.1 immediately implies that the set $\mathcal{H}=\{H\in\{-1,0,1\}^{k\times d}:\forall i\in[k],\|\mathbf{H}_{[i,:]}\|_{0}=\left\|\mathbf{H}_{[i,:]}\right\|_{\infty}=1\}$ is the set of extreme points of the set $\mathcal{G}$. (If the set of extreme points of $\mathcal{G}$ is larger than $\mathcal{H}$, then the signs of any such extreme point can be flipped to create a new extreme point of row-stochastic matrices, which would violate H.1.)

It is not hard to observe that for any $H\in\mathcal{H}$, there exists an $\mathbf{u}_{\mathbf{H}}\{\pm 1\}^{k}$ s.t. $\left\|\mathbf{C}\mathbf{H}\right\|_{F}=\left\|\mathbf{C}\mathbf{u}_{\mathbf{H}}\right\|_{2}$. Since, $\left\|\mathbf{C}\mathbf{u}_{\mathbf{H}}\right\|_{2}\leq\max\limits_{\mathbf{u}\in\{\pm 1\}^{k}}\left\|\mathbf{C}\mathbf{u}\right\|_{2}$ for any choice of $\mathbf{u}_{\mathbf{H}}$, and the fact that $\max\limits_{\mathbf{G}}\left\|\mathbf{C}\mathbf{G}\right\|_{F}$ is reached at one of the matrices in $\mathcal{H}$, the claim in Theorem H.1 follows. ∎

H.2 A counterexample for general 𝐂\mathbf{C}

Direct calculation shows $\max_{\mathbf{u}\in\{\pm 1\}^{k}}\left\|\mathbf{C}\mathbf{u}\right\|_{2}=1$, but $\left\|\mathbf{C}\mathbf{G}\right\|_{F}=\sqrt{1.1}\approx 1.049$.

H.3 Proof of 2.1

When per-step contributions are bounded by $\zeta=1$, for any participation pattern $\Pi$ and dimensionality $d\geq 1$, when $\mathbf{C}^{\top}\mathbf{C}\geq 0$ elementwise, we have

For the other direction, for each $\pi\in\Pi$, we apply H.1 to the matrix $\mathbf{C}_{\pi}=\mathbf{C}_{[:,\pi]}$, and observe $\mathbf{C}^{\top}\mathbf{C}\geq 0$ is sufficient to imply $\mathbf{C}^{\top}_{\pi}\mathbf{C}_{\pi}\geq 0$. ∎

The condition $\mathbf{C}^{\top}\mathbf{C}\geq 0$ is sufficient but not in fact necessary for 2.1 to hold. In particular, for $(k,b)$-participation $\Pi$, the sub-matrices $\mathbf{C}^{\top}_{\pi}\mathbf{C}_{\pi}$ for $\pi\in\Pi$ “touch” only $k^{2}b$ entries of the $n^{2}=k^{2}b^{2}$ entries of $\mathbf{C}^{\top}\mathbf{C}$; the other entries of $\mathbf{C}^{\top}\mathbf{C}$ could in fact be negative. However, we did not need to use this observation for any of the matrices in our experiments.

H.4 Proof of 2.1

By assumption we have $\|\mathbf{u}\|\leq\sqrt{k}$, and so

Appendix I Analysis for Section 3

Then, for Lagrange multipliers $\mathbf{v}$ such that the $\mathbf{U}$ is full-rank, the Lagrange dual function $g$ can be expressed in closed form in terms of the Lagrange multipliers:

Recall that we have defined $\mathbf{H}_{\mathbf{v}}=[\mathbf{u}_{1},\dots,\mathbf{u}_{k}]\diag(\mathbf{v})^{\frac{1}{2}}$, and $\mathbf{U}=\mathbf{H}_{\mathbf{v}}\mathbf{H}_{\mathbf{v}}^{\top}$. Now, note:

Introducing Lagrange multipliers $\mathbf{v}_{\mathbf{u}}\geq 0$, for the problem Section I.1 we form the Lagrangian for positive-definite $\mathbf{X}$:

For fixed $\mathbf{v}$, any finite minimizer of $L$ for positive-definite $\mathbf{X}$ must correspond to a zero of this Lagrangian’s gradient. We then compute the gradient

$\mathbf{U}$ and $\mathbf{A}$ are full-rank by assumption; therefore I.2 is applicable, and Eq. 23 has a unique positive-definite zero (and indeed, the infimum in Section I.1 becomes a minimum):

Note that Eq. 23 also immediately implies that if $\mathbf{U}$ is not full-rank, then there is no finite positive-definite minimizer of $L$ in $\mathbf{X}$. Letting $g(\mathbf{v})=\min_{\mathbf{X}}L(\mathbf{X},\mathbf{v})$ be the Lagrange dual function and plugging back into Eq. 22, we have

In the same setup as 3.1, a maximizer of the dual ${\mathbf{v}^{\star}}$ must satisfy:

As noted in the remark after 3.1, any optimal setting of the dual variables must be in the interior of a neighborhood in which the representation Eq. 7 is valid. It is therefore permissible to differentiate this representation.

by defining $\mathbf{X}=\mathbf{U}^{-\frac{1}{2}}\big(\mathbf{U}^{\frac{1}{2}}\mathbf{A}^{\top}\mathbf{A}\mathbf{U}^{\frac{1}{2}}\big)^{\frac{1}{2}}\mathbf{U}^{-\frac{1}{2}}$ (recalling the usage of the symbol $\mathbf{X}$ in Eq. 24).

and we have the stated expression for the gradient of the dual function.

Now, at a maximizer of the dual function, this derivative must vanish. An equivalent condition is $\diagpart(\mathbf{H}^{\top}\mathbf{X}\mathbf{H})=\vec{1}$, and hence

so at the optimum ${\mathbf{v}^{\star}}$ in fact $g({\mathbf{v}^{\star}})=\sum_{\mathbf{u}}{\mathbf{v}^{\star}}_{\mathbf{u}}$, establishing the second claim of our result.

Again using the observation that $\diagpart(\mathbf{H}_{\mathbf{v}}^{\top}\mathbf{X}\mathbf{H}_{\mathbf{v}})=\vec{1}$ and so

Further, using the second claim of I.1, we can take

and multiplying this by $\mathbf{H}_{\mathbf{v}}^{\top}$ and $\mathbf{H}_{\mathbf{v}}$ on the left and right respectively yields

and so we conclude for the optimal Lagrange multiplier ${\mathbf{v}^{\star}}$,

I.2 Lemmas and Corollaries

where $\mathbf{e}_{i}$ is the $i^{th}$ standard basis vector.

Then for $\mathbf{X}$ satisfying our assumptions,

where the final inequality follows by the assumptions on $\mathbf{X}$.

Let $\mathbf{U},\mathbf{V}\in S^{n}_{++}$. Let $\mathbf{U}=\mathbf{U}_{L}\mathbf{U}_{R}$ be a factorization of $\mathbf{U}$ such that $\mathbf{U}_{R}\mathbf{V}\mathbf{U}_{L}$ is PSD, and the following equation defines a positive-definite matrix $\mathbf{X}$:

Then, this $\mathbf{X}$ solves the equation

Moreover, this positive-definite solution $\mathbf{X}$ is unique.

We will begin by showing that $\mathbf{X}$ as defined by Eq. 31 solves the equation $\mathbf{X}\mathbf{U}\mathbf{X}=\mathbf{V}$; then we will show that any two positive-definite representations of the form Eq. 31 are in fact identical.

Notice that the representation $\mathbf{U}=\mathbf{U}_{L}\mathbf{U}_{R}$ implies that $\text{rank}(\mathbf{U}_{L})\geq n$ and $\text{rank}(\mathbf{U}_{R})\geq n$. Therefore $\mathbf{U}_{L}\mathbf{U}_{L}^{\dagger}=\mathbf{I}=\mathbf{U}_{R}^{\dagger}\mathbf{U}_{R}$, as implied by the Moore definition of the Moore-Penrose pseudoinverse. So:

We claim that $\big(\mathbf{U}_{R}\mathbf{V}\mathbf{U}_{L}\big)^{\frac{1}{2}}\mathbf{P}_{R(\mathbf{U}_{L}^{\dagger})}=\mathbf{P}_{R(\mathbf{U}_{R})}\big(\mathbf{U}_{R}\mathbf{V}\mathbf{U}_{L}\big)^{\frac{1}{2}}=\big(\mathbf{U}_{R}\mathbf{V}\mathbf{U}_{L}\big)^{\frac{1}{2}}$. In both cases, this can be seen by multiplying on the left or right as appropriate by $\big(\mathbf{U}_{R}\mathbf{V}\mathbf{U}_{L}\big)^{\frac{1}{2}}$, and noting $\mathbf{P}_{R(\mathbf{U}_{R})}\mathbf{U}_{R}=\mathbf{U}_{R}$, $\mathbf{U}_{L}\mathbf{P}_{R(\mathbf{U}_{L}^{\dagger})}=\mathbf{U}_{L}$. Since all the terms here are symmetric, the appropriate equality follows by uniqueness of the symmetric matrix square root. Therefore:

The uniqueness of a positive-definite $\mathbf{X}$ solving $\mathbf{X}\mathbf{U}\mathbf{X}=\mathbf{V}$ follows from the uniqueness of the usual matrix square root. Indeed, assume $\mathbf{Y}$ positive-definite satisfies $\mathbf{Y}\mathbf{U}\mathbf{Y}=\mathbf{V}$. Then:

Since the positive-definite square root is uniquely determined, $\mathbf{U}^{1/2}\mathbf{Y}\mathbf{U}^{1/2}$ is uniquely determined. Since $\mathbf{U}$ is invertible, $\mathbf{Y}$ is uniquely determined as well, and we have $\mathbf{Y}=\mathbf{X}$. ∎

Two particular instantiations of I.2 are of interest. $\mathbf{X}$ as the matrix geometric mean of $\mathbf{U}^{-1}$ and $\mathbf{V}$ (taking $\mathbf{U}_{L}=\mathbf{U}_{R}=\sqrt{\mathbf{U}})$:

and assuming the representation $\mathbf{U}=\mathbf{H}_{\mathbf{v}}\mathbf{H}_{\mathbf{v}}^{\top}$:

By positive-definiteness of $\mathbf{U}$ and $\mathbf{V}$, Eq. 32 is clearly positive definite; Eq. 33 may be seen to be positive definite via the SVD of the pseudoinverses involved. Symmetry is again clear. Therefore both representations satisfy the assumptions of I.2. ∎

I.3 Impact of non-negativity constraints

We consider three variations of the constraints in Eq. 6. In the first, we have the default constraints

This requires enumerating the $b2^{k}$ corners of $\mathcal{D}^{1}_{\Pi}$, which is tractable for the small problems we consider here. Next, we consider the approach used for our experiments:

Note here we only have $b$ non-negative vectors $\mathbf{u}$, and so only $b+n^{2}$ constraints.

Finally, we consider and “in between” set of constraints. This is the minimal set of constraints that allows us to apply claim 2 of H.1. In order to define these constraints, let $\widehat{\Pi}=\big\{(i,j)|(i,j)\in\pi,i\neq j,\pi\in\Pi\big\}$, the set of pairs of distinct iterations in which the same user could possibly participate. For Claim 2 of H.1, because we apply the theorem to sub-matrices $\mathbf{C}_{[:,\pi]}$, a sufficient condition is $\mathbf{X}_{[i,j]}\geq 0$ for all $(i,j)\in\hat{\Pi}$. Then, the constraints become:

Since $|\widehat{\Pi}|=bk(k-1)$, we have $b+bk(k-1)<b+n^{2}$ constraints.

In Table 6 we show results for a tiny problem: $n=6$, $k=3$, and $b=2$. However, we have run additional experiments that align with the conclusions reported here for moderately larger problems (noting the full approach quickly becomes prohibitively expensive). For the prefix-sum workload, the choice of additional constraints does not impact the total error, as in all cases $\mathbf{X}\geq 0$. However, for the momentum workload (with momentum 0.95), we see small gaps, indicating that imposing the additional constraints does come at at a small (from 16.115 to 16.134) impact on the total error for this workload.

Appendix J Analysis for Section 4

We consider the special case where $\mathbf{A}$ is the prefix sum linear query matrix (lower-triangle matrix of ones). Then, we define the corresponding circulant matrix

It is straightforward to verify ${\mathbf{A}_{\sf circ}}_{[:n,:n]}=\mathbf{A}$.

In the following we provide the privacy guarantee and the main utility guarantee for the FFT mechanism defined in Algorithm 1.

Algorithm 1 is $\rho$-zCDP in the adaptive continuous release model.

Next, we analyze the utility of Algorithm 1 and show that it is nearly optimal in terms of the mean squared error (MSE) in the single-pass setting. First, we express the MSE in Theorem J.3 below.

The MSE achieved by Algorithm 1 using the real and imaginary components of $\widetilde{}\mathbf{z}$ is

In the following, we will have an explict expression for $\left\|\mathbf{v}^{\sf DFT}\right\|_{1}$ in terms of the problem parameters. Finally, we will argue that Theorem J.3 is nearly optimal.

The expected mean squared error (MSE) is given by the following:

Here, we show that Theorem J.3 is near-optimal in utility for the single-participation setting. To do this, we compare with a lower bound on the expected MSE of any factorization-based mechanism from Henzinger et al. 2022: $\frac{1}{2\rho\pi^{2}}\left(2+\ln\left(\frac{2n+1}{3}\right)+\frac{\ln(2n+1)}{2n}\right)^{2}$. We find that the though our analytical upper bound in Corollary J.1 is $\approx 6$x worse than the lower bound, the empirical noise added in Algorithm 1 closely tracks the lower bound to within a factor of $1.2$x—because it only adds the real part of the noise. Results are in Figure 14 of Appendix J.1.

Showing near-optimal utility via MSE experiments

J.2 Proof of J.2

Algorithm 1 is $\rho$-zCDP in the adaptive continuous release model.

First, consider the non-adaptive setting and the following mechanism, with parameters as defined in Algorithm 1,

We claim that this satisfies $\frac{\kappa^{2}\left\|\mathbf{v}^{\sf DFT}\right\|_{1}}{4n\sigma^{2}}$-zCDP. To see this, we proceed by bounding $\rho_{i}$ for each coordinate $i\in[2n]$ defined in Equation 39. For brevity, let $\mathbf{b}=\mathbf{F}\mathbf{x}_{\sf ext}$. Consider two neighboring data sets $\mathbf{g}$ and $\mathbf{g}^{\prime}$, correspondingly, $(\mathbf{b},\mathbf{x}_{\sf ext})$ and $(\mathbf{b}^{\prime},\mathbf{x}_{\sf ext}^{\prime})$. Then,

We will now prove zCDP guarantee independently for each of the $2n$ coordinates and then use standard zCDP composition (Bun & Steinke 2016). For any coordinate $a\in\{0,\ldots,2n-1\}$, adding noise $\left(\frac{\sigma}{\sqrt{|\mathbf{v}^{\sf DFT}[i]|}}\right)\cdot\mathcal{N}_{\sf complex}\left(0,1\right)$ to $\mathbf{b}[i]$ satisfies $\rho_{i}$-zCDP with $\rho_{i}=\frac{\kappa^{2}|\mathbf{v}^{\sf DFT}[i]|}{4n\sigma^{2}}$. Then by composition, we have that

Therefore, setting $\sigma^{2}=\frac{\kappa^{2}\left\|\mathbf{v}^{\sf DFT}\right\|_{1}}{4n\rho}$-satisfies a non-adaptive $\rho$-zCDP. Using the same $\sigma$, we prove the adaptive part using the same $\sigma$. We have the following from Equation 39.

Since, $\mathbf{w}$ in Equation 42 is spherical Gaussian, and the original query matrix $\mathbf{A}$ is lower triangular, by Theorem 2.1 in Denisov et al. 2022, the adaptive privacy guarantee follows. ∎

J.3 Proof of J.3

The MSE achieved by Algorithm 1 using the real and imaginary components of $\widetilde{}\mathbf{z}$ is

In equation 43, $\Sigma[:n,:n]$ refers to the top-left $n\times n$ submatrix of $\Sigma$. ∎

J.4 Proof of J.1

Under the same setting as Theorem J.3, the MSE for Algorithm 1 is the following

Recall the definition of DFT from Equation 38 and of $\mathbf{v}$ in Equation 37. It is immediate that $\mathbf{v}^{\sf DFT}=n$. For any $k\neq 0$, we have,

From Equation 44, we have that when $k>0$ is even, $\mathbf{v}^{\sf DFT}[k]=0$. For $k$ odd, we have

Combining these, the term $\left\|\mathbf{v}^{\sf DFT}\right\|_{1}$ is

J.5 Proof of 4.1

Under $k$ participation, Algorithm 1 satisfies $(k^{2}\rho)$-zCDP.

The proof goes exactly as Theorem J.2, except equation 40 gets replaced by the following:

Appendix K Two related FFT mechanisms.

The FFT mechanism presented in Section 4 can be understood as an application of a complex-valued matrix mechanism factorizing the prefix-sum matrix as

All these operations are linear; and since everything begins and ends in the real domain, this mechanism can be expressed as a real-valued mechanism. Therefore identical codepaths can be used for implementing experiments with the FFT, though notably without some special implementation of the mechanism, realizing the potential computational savings will not be immediate. In this small section, we translate this complex-valued mechanism into two real-valued mechanism which can be integrated with the code backing the rest of the paper. These mechanisms differ in their decoding matrix $\mathbf{B}$, and thus achieve different levels of loss. Both have efficient implementations, though with asymptotics differing by a logarithmic factor. We implement and experiment with both of these mechanisms, though we only report results from the mechanism with lower loss in the main body.

These two mechanisms share an encoding matrix:

which is real-valued by K.1. Note that the sensitivity of $\mathbf{C}_{\mathbf{F}}$ is identical to that of $\mathbf{C}$ for any notion of sensitivity expressible as Definition 1 due to the unitary of the Fourier transform. Since this matrix is of shape $[2n,n]$, there is choice in computing the decoder $\mathbf{B}$ such that $\mathbf{B}\mathbf{C}_{\mathbf{F}}$ represents the prefix-sum matrix. The two decoders we present below correspond to two subtly distinct mechanisms.

One natural translation of the analysis in Section 4 (indeed, a real-valued version of the precise operation described in Algorithm 1) may be computed by inserting a Fourier transform to match the inverse transform in $\mathbf{C}_{\mathbf{F}}$:

Clearly $\mathbf{B}_{\mathbf{F}}\mathbf{C}_{\mathbf{F}}=\mathbf{B}\mathbf{C}$, and $\mathbf{B}_{\mathbf{F}}$ real-valued by K.1.

For any $\mathfrak{D}$, the mechanism described in Section 4 is distributionally equivalent to an application of the real-valued matrix mechanism with the factorization $\left(\mathbf{B}_{\mathbf{F}},\mathbf{C}_{\mathbf{F}}\right)$, and satisfies the same privacy guarantees.

To show this result, by noting that $\mathbf{C}_{\mathbf{F}}$ and $\mathbf{C}$ have the same sensitivity, it suffices to show that:

$\Re[\mathbf{F}^{*}\sqrt{\Sigma}\mathbf{z}]$ (for $\mathbf{z}$ a sample from an isotropic complex Gaussian) is distributionally equivalent to $\mathbf{P}\mathbf{F}^{*}\sqrt{\Sigma}\mathbf{F}\mathbf{b}$ for $\mathbf{b}$ a sample from a real (isotropic) Gaussian with the same variance.

This is a consequence of the distributional invariance of the Gaussian under unitary transformations:

Note that the efficiency of the mechanism described in Section 4 carries over immediately to this factorization $\left(\mathbf{B}_{\mathbf{F}},\mathbf{C}_{\mathbf{F}}\right)$; indeed, the capacity to compute the noise $\mathbf{B}_{\mathbf{F}}\mathbf{b}$ with complexity $n\log(n)$ may be reasoned to directly, in a similar manner.

This mechanism is not, however, the optimal one for the encoder $\mathbf{C}_{\mathbf{F}}$, and this subtlety has difficult downstream effects in integrating with real-valued factorization codepaths (e.g., see the discussion in Section D.4). We proceed to show that the optimal decoder can be used directly, at only a moderate loss of efficiency with sufficiently careful implementation.

As noted in the literature (e.g. Section 3 of Denisov et al. 2022), for a fixed encoder, the optimal decoder may always be computed in terms of an appropriate pseudoinversion of the encoder. Therefore, we may compute the optimal decoder for the encoder $\mathbf{C}_{\mathbf{F}}$, defining:

where $\mathbf{S}$ is the prefix-sum matrix. Since $\mathbf{C}_{\mathbf{F}}$ is real, its pseudoinverse is as well, and $\mathbf{B}_{\mathbf{F}\text{opt}}$ is also real-valued. Since $\mathbf{B}_{\mathbf{F}\text{opt}}$ can have no more variance than $\mathbf{B}_{\mathbf{F}}$, all the utility analysis of the DFT mechanism in Section 4 carries through as an upper bound for this factorization. Privacy of this mechanism is ensured by the fact that this mechanism reuses the encoder $\mathbf{C}_{\mathbf{F}}$. The major way in which these mechanisms operationally differ comes down to the cost of computing the noise vector $\mathbf{B}_{\mathbf{F}\text{opt}}\mathbf{b}$, where $\mathbf{b}$ represents a sample from an isotropic Gaussian distribution. Though we do not know of a complexity result which matches the decoder $\mathbf{B}_{\mathbf{F}}$, we will show that the complexity cost which must be paid is only logarithmically higher.

First, notice that the matrix $\mathbf{C}_{\mathbf{F}}$ is one-to-one; indeed, this is immediately implied by the factorization $\mathbf{S}=\mathbf{B}_{\mathbf{F}}\mathbf{C}_{\mathbf{F}}$. By Theorem 1.2.1 (P6) of (Campbell & Meyer 1979), any one-to-one matrix $\mathbf{T}$ admits the following representation for its pseudoinverse:

Now, the matrix $\mathbf{P}\mathbf{F}^{*}|\Sigma|\mathbf{F}\mathbf{E}$ is Toeplitz, since $\mathbf{F}^{*}|\Sigma|\mathbf{F}$ is circulant, and $\mathbf{P}$, $\mathbf{E}$ combine to select out the top-left $n\times n$ square of $\mathbf{F}^{*}|\Sigma|\mathbf{F}$. Notice that $\mathbf{P}\mathbf{F}^{*}|\Sigma|\mathbf{F}\mathbf{E}$ is not circulant, and cannot therefore be diagonalized by the $n$-dimensional Fourier transform.

The development of Section 4 yield the representation:

which implies that matrix-vector products with the matrix $\mathbf{S}$ may be computed in $n\log n$ time by the use of the FFT. Similarly, matrix-vector products with $\left(\mathbf{F}^{*}\sqrt{\Sigma}\mathbf{F}\mathbf{E}\right)^{*}$ may be computed in $n\log n$ time.

Therefore the computational cost of computing the mapping $\mathbf{b}\mapsto\mathbf{S}\mathbf{C}_{\mathbf{F}}^{\dagger}\mathbf{b}$ can be upper bounded by the maximum of $n\log n$ and the cost of computing the mapping $\mathbf{v}\mapsto\left(\mathbf{P}\mathbf{F}^{*}|\Sigma|\mathbf{F}\mathbf{E}\right)^{-1}\mathbf{v}$.

The cost of computing this mapping is, in turn, bounded by the cost of inverting a general (full-rank) Toeplitz system, since $\left(\mathbf{P}\mathbf{F}^{*}|\Sigma|\mathbf{F}\mathbf{E}\right)^{-1}\mathbf{v}$ may be alternatively characterized as the solution $\mathbf{x}$ to the equation $\mathbf{P}\mathbf{F}^{*}|\Sigma|\mathbf{F}\mathbf{E}\mathbf{x}=\mathbf{v}$. The computational cost of solving such a system is known to be $n\log^{2}(n)$; see, e.g., (de Hoog 1987). ∎

For a real-valued vector $\mathbf{v}$, let $\hat{\mathbf{v}}$ represent its discrete Fourier transform. If $\hat{\mathbf{v}}$ has no purely real, negative entries, then letting $\sqrt{\cdot}$ denote the (pointwise) principal branch of the square root and $\mathbf{F}$ the matrix representation of the Fourier transform, the matrix $\mathbf{F}^{*}\sqrt{\hat{\mathbf{v}}}\mathbf{F}$ is real-valued.

Conjugate symmetry of the DFT states that for a $j$-dimensional real-valued vector $\mathbf{x}$, $\hat{\mathbf{x}}[m]=\overline{\hat{\mathbf{x}}}[j-m]$, and that the converse also holds–that if $\hat{\mathbf{x}}$ has this symmetry, $\mathbf{x}$ is real-valued. This can be seen by examining the action of conjugation of $\mathbf{x}$ on the Fourier transform $\hat{\mathbf{x}}$.

Now, by the assumptions on $\hat{\mathbf{v}}$ and the choice of the principal branch of the square root These assumptions can be avoided, though at the cost of taking care in choosing the square root of the negative elements of $\hat{\mathbf{v}}$ to preserve the appropriate symmetry., if $\hat{\mathbf{v}}$ has this conjugate symmetry, so does $\sqrt{\hat{\mathbf{v}}}$. Therefore there is some real-valued vector $\mathbf{y}$ such that $\hat{\mathbf{y}}=\sqrt{\hat{\mathbf{v}}}$. The matrix $\mathbf{F}^{*}\sqrt{\hat{\mathbf{v}}}\mathbf{F}$ represents convolution with $\mathbf{y}$ in the standard basis, and hence is real-valued. ∎