Efficient and Near-Optimal Noise Generation for Streaming Differential Privacy
Krishnamurthy Dvijotham H. Brendan McMahan Krishna Pillutla Thomas Steinke Abhradeep Thakurta



    
    Alphabetical author order.

Efficient and Near-Optimal Noise Generation
for Streaming Differential Privacy
    
    Krishnamurthy (Dj) DvijothamGoogle DeepMind 
    H. Brendan McMahanGoogle Research 
    Krishna PillutlaIIT Madras. Work done at Google. 
    Thomas SteinkeGoogle DeepMind 
    Abhradeep ThakurtaGoogle DeepMind 
    
    April 29, 2024
=================================================================================================================================================================================================================================




    In the task of differentially private (DP) continual counting, we receive a stream of increments and our goal is to output an approximate running total of these increments, without revealing too much about any specific increment.
    Despite its simplicity, differentially private continual counting has attracted significant attention both in theory and in practice.
    Existing algorithms for differentially private continual counting are either inefficient in terms of their space usage or add an excessive amount of noise, inducing suboptimal utility.
    
    The most practical DP continual counting algorithms add carefully correlated Gaussian noise to the values. The task of choosing the covariance for this noise can be expressed in terms of factoring the lower-triangular matrix of ones (which computes prefix sums). We present two approaches from this class (for different parameter regimes) that achieve near-optimal utility for DP continual counting and only require logarithmic or polylogarithmic space (and time).
    
    Our first approach is based on a space-efficient streaming matrix multiplication algorithm for a class of Toeplitz matrices. 
    We show that to instantiate this algorithm for DP continual counting, it is sufficient to find a low-degree rational function that approximates the square root on a circle in the complex plane. We then apply and extend tools from  approximation theory to achieve this.
    We also derive efficient closed-forms for the objective function for arbitrarily many steps, and show direct numerical optimization yields a highly practical solution to the problem.
    Our second approach combines our first approach with a recursive construction similar to the binary tree mechanism.








§ INTRODUCTION


The simplest task in differentially private data analysis is counting.
This task has an equally simple (and optimal) algorithm: we simply add Laplace or Gaussian noise to the final count <cit.>. 

Continual counting <cit.> asks that we not only release the final count, but that we provide a running total, i.e., we release all the partial sums.
This task and its variants have proved to be a surprisingly complex problem in the theory of differential privacy <cit.>.

Continual counting is also critically important in practice. 
State-of-the-art methods for differentially private machine learning rely on (high-dimensional extensions of) continual counting algorithms to maintain their state in a differentially private manner <cit.>. (See <ref> for discussion of this application.)
This has motivated significant work on making differentially private continual counting as accurate as possible <cit.>.

Practical algorithms for continual counting add Gaussian noise to the partial sums. However, this noise is not independent: it must be carefully correlated across steps, as we will see. Choosing an appropriate multivariate Gaussian can be formulated as a matrix factorization problem, as we now describe.

To convey the key ideas in this introduction, we consider summing  scalars; this discussion readily extends to summing vectors.
Let x ∈ℝ^ contain the  terms we wish to sum (over  iterations or steps).
Let A ∈{0,1}^× be the lower-triangular all-ones matrix so that Ax is the vector of partial sums – i.e., (Ax)_k = ∑_i ≤ k x_i for each k. We factor A=BC, where B,C^T ∈ℝ^×'. 
The corresponding differentially private continual counting algorithm  is given by

    (x)  B(Cx+z) = Ax + Bz = A(x+C^† z),

where z 𝒩(0,σ^2I) is a vector of independent Gaussian noise.[C^†=A^-1B is the inverse of C whenever C is a square matrix. However, we will also consider settings where C is not a square matrix and so C^† is a suitable pseudo-inverse.]
Equivalently, the output of the mechanism can be written as (x)=Ax+z where z𝒩(0,Σ) with Σ = σ^2 B B^T.



  
Mechanism privacy: We assume that one person can change only one entry of the input x by at most 1, which means Cx can change by at most the maximum column norm of C, denoted  

    C_1 → 2max_j √(∑_i Cij^2).

This is the L_2 sensitivity of Cx. Then, following the standard privacy analysis of Gaussian noise addition <cit.>, to ensure that Cx+z is (ε,δ)-differentially private, we scale the noise z 𝒩(0,σ^2I) to have standard deviation  σ = ζC_1→2, where ζ = O(1/ε√(log(1/δ))) is a noise multiplier depending only on the privacy parameters (ε, δ).[Importantly, particularly for machine learning applications, this privacy analysis can be extended to the case where each x_k is adaptively chosen depending on the prefix sums 0, …, k-1 already released by the mechanism <cit.>.] By postprocessing, this implies (x)=B(Cx+z) is (ε,δ)-differentially private.




  
Mechanism utility:
The error of the mechanism  is given by (x) - Ax = Bz. Thus, the root mean squared error of the ith partial sum can be calculated as 

    √(𝔼_((x)-Ax)_i^2)=√(𝔼_z[(Bz)_i^2]) = √(σ^2  ∑_j B_i,j^2),

which scales with the norm of the corresponding row of B.
We take our notion of utility to be the maximum[We can also consider other measures of error such as sum of variances (rather than the max). This corresponds to a different norm: Bz_2^2 = B^2 ·σ^2.] such error over all  partial sums released, which scales with the maximum row norm 

    B_2 →∞max_i √(∑_j B_i,j^2).















Thus, the maximum root mean squared error of the mechanism's answers is

    max_i √(((x)-Ax)_i^2)   = max_i √(z 𝒩(0,σ^2I)(Bz)_i^2) = max_i √(σ^2 ∑_j B_i,j^2)
       = σB_2→∞ = ζB_2→∞C_1→2 =: ζ(B,C),

where we define

    (B, C) 
      B_2 →∞C_1→2
      = √(max_i ∑_j B_i,j^2 )√(max_j ∑_i C_i,j^2 ).


Since the error in <ref>  can be written as a product of the noise multiplier ζ (which does not depend on the factorization B,C) and (B,C) (which does not depend on the differential privacy parameters), the same factorization will minimize error for all settings of the privacy parameters. Thus, we suppress ζ for the remainder of this paper.  




  
Matrix factorization:
To summarize, our goal is to solve the matrix factorization problem

    minimize_B,C (B, C)       subject to  B C = A.

The optimal value of this objective is denoted γ_2(A) and is known as the gamma-two factorization norm of A <cit.>.

Prior work has obtained near-optimal factorizations of A (and related matrices). The optimal value is 

    γ_2(A) inf{(B,C) : B, C ∈ℝ^×, BC=A } = log()/π± O(1).

However, the existing factorizations are either far from optimal or do not readily yield practical algorithms because generating the noise B z is not sufficiently computationally efficient.
That is the problem we address in this paper.
The main computational challenge is that our algorithm cannot store x or z (or even B or C) in memory; this would require memory that is linear in the number of steps . 

We provide a matrix factorization mechanism for differentially private continual counting attaining near-optimal error with only polylogarithmic memory overhead.

To give a sense of scale, in private machine learning applications, we often need to compute >10^6 partial sums (each coming from an iteration of stochastic gradient descent or a related algorithm), over >10^10 dimensions (the size of the model and its gradients).[In this introductory discussion, we have considered sums over =1 dimension. Extending to >1 is essentially a matter of running  parallel copies of the algorithm. In our technical sections, we take x to be a matrix whose rows we wish to sum, rather than a vector.]
Thus we need to be able to compute M(x) without storing all of the data or noise in memory, as this would take >10^16 units of memory – i.e., over 40 petabytes if we use 32-bit floating point numbers.[Throughout this paper we treat real computation as atomic; e.g., a real number takes up one unit of memory. In practice, floating point implementations can compromise privacy <cit.>, but techniques exist to ensure differential privacy in discrete settings <cit.>.]



  
Streaming setting: Our algorithm receives the coordinates of the input x and noise z one at a time and outputs the approximate partial sums (x) one at a time. In particular, we must output each partial sum (x)_k ≈∑_i ≤ k x_i before receiving the next input x_k+1.

Ideally, the memory usage should not grow with the stream length  at all. 
If we use the trivial factorization given by B=I and C=A, then this corresponds to adding independent noise to each output. This would only require constant memory and constant time per output, but the noise scale (B,C)=√() is far from optimal. 
Thus the challenge is to simultaneously obtain near-optimal error and computational efficiency.

To overcome this challenge we must impose some structure on the factors B and C that allows efficient noise generation without significantly increasing the matrix factorization objective (B,C).
The three structures that have been considered in the literature are (i) lower triangular matrices, (ii) sparse matrices, and (iii) Toeplitz matrices – i.e., constant diagonals, ∀ i,j    C_i,j = c_i-j for some vector c – and combinations of these structures.


We focus primarily on lower triangular Toeplitz matrices. 
The streaming setting naturally corresponds to lower triangular matrices, as it ensures that the k-th output is only a function of the first k inputs.
Adding the Toeplitz constraint has a minimal impact on the matrix factorization objective (<ref>). Furthermore, Toeplitz matrices (and certain specializations we consider) have convenient mathematical and algorithmic properties. See <ref> and <ref> for further discussion about lower triangular Toeplitz factorizations versus other factorizations. We also leverage the sparsity in the B and C matrices derived from generalizations of the binary tree aggregation mechanism to show additional results.




 §.§ Our Contributions


Our main theoretical result is a lower triangular Toeplitz matrix factorization that is nearly optimal accompanied by an efficient streaming algorithm for generating the corresponding noise.


    For all ≥ 1 and all μ∈ (0,1)
    , there exists a lower triangular Toeplitz matrix factorization B, C ∈ℝ^× with the following properties.
    
        
  * Validity:
        BC=A, where A is the × lower-triangular all-ones matrix.
        
  * Near-optimality:
        
        (B,C) ≤𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾() + μ,
        
        
        where 
        
    𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾() = 1+∑_k=1^-1( 2^-2k2k  k)^2 ≤ 1 + 0.57722+log(n)/π
 is the optimal value of (B,C) over all lower triangular Toeplitz factorizations A=BC.
        
  * Efficiency:
        There exists a streaming algorithm that at each step k takes as input z_k and outputs (Bz)_k (or, equivalently, (C^-1z)_k) and runs in time and space O(log^2(/μ)).
    


In order to ensure that the error term μ in our result is o(1), we have space and time complexity d=Θ(log^2 ).
It is natural to wonder whether this computational complexity can be improved.
We show that it can be improved to O(log) at the expense of a weaker multiplicative near-optimality guarantee for the matrix factorization objective.


    For all ≥ 1, there exists a matrix factorization B, C^T ∈ℝ^×' with '=O() and the following properties.
    
        
  * Validity:
        BC=A, where A is the × lower-triangular all-ones matrix.
        
  * Near-optimality:
        (B,C) ≤ (1+o(1)) ·𝖮𝗉𝗍() ,
        where 𝖮𝗉𝗍() = log(n)/π± O(1) is the optimal value of (B,C) over all factorizations.
        
  * Efficiency:
        There exists a streaming algorithm that at each step k takes as input some coordinates of z (but never reads the same coordinate more than once) and outputs (Bz)_k and runs in space (and amortized time per iteration) O(log).
    


More generally, we can smoothly trade off between the matrix factorization objective and computational efficiency. That is, we can interpolate between <Ref>; see <Ref> for a general statement.

<Ref> does not generate a lower triangular Toeplitz factorization; in fact it does not even produce a square factorization. It also gives a weaker near-optimality guarantee. We leave it as an interesting open problem whether it is possible to improve on log^2 space complexity with a Toeplitz factorization or with (B,C) ≤𝖮𝗉𝗍()+O(1).







  
Practical algorithms:
Since our work is directly motivated by practical considerations, we also study the problem numerically.
Specifically, we show that the method behind <Ref> yields a factorization that is indistinguishable from optimal Toeplitz mechanism for practical purposes, but which has a highly efficient procedure for generating the noise. 
<Ref> shows the ratio of (B,C) for our algorithm over the optimal 𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾().

To prove <ref> we provide closed-form parameters for the matrix factorization. This is already reasonably practical, but we further optimize the parameters numerically to obtain even better algorithms.
In order to optimize the matrix factorization we (i) choose an appropriate parameterization for the class of factorizations B,C that we consider, (ii) give an efficiently computable expression for the objective (B,C) in terms of this parameterization, and (iii) show that this expression is differentiable. This allows us to numerically optimize our factorization. While we do not prove that this optimization procedure converges, in practice it yields significantly better solutions than the closed-form parameters of <ref> 

To show the power of our approach, in <ref>, we show that even with d=1 buffer, we can asymptotically attain (B,C)=O(n^1/6).




  
Lower bounds:
In <ref> we prove a lower bound that exactly characterizes 𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(). Specifically, we show that the lower triangular Toeplitz factorization of <cit.> is precisely optimal for this class.

Further, in <ref>, we develop numerical lower bounds  on the objective of any matrix factorization for various classes of matrices for any fixed size n. In particular, we give lower bounds for arbitrary lower triangular matrices, arbitrary Toeplitz matrices, and Toeplitz matrices that are realized from constant degree . 









 §.§ Our Algorithms & Techniques




  
: Buffered Linear Toeplitz Matrices (<ref>):

The starting point for our main result is the lower triangular Toeplitz factorization of <cit.>, which we show is in fact optimal in <Ref>.
This factorization is given by B_i,j=C_i,j=f_i-j where the sequence f_0, f_1, ⋯ is the coefficients of a Taylor series:

    1/√(1-x) = f_0 + f_1 x + f_2 x^2 + f_3 x^3 + ⋯.


This sequence satisfies the recurrence f_k=(1-1/2k) · f_k-1 for k ≥ 1 with f_0=1. Our task is to compute
(Bz)_k = ∑_j=0^k B_k,j· z_j = ∑_j=0^k f_k-j· z_j 
in a streaming fashion for k ∈ [] {0,1,⋯,-1}.[We zero-index sequences, vectors, and matrices throughout. See <ref> for a summary of symbols and notation.] Unfortunately, the factorization of <ref> does not seem to admit an efficient sampling procedure. 



Suppose the sequence instead satisfied a linear recurrence f_k = q · f_k-1 for all k ≥ 1.
Then

    (Bz)_k = ∑_j=0^k f_k-j· z_j = f_0 · z_k + ∑_j=0^k-1 (q · f_k-j-1) · z_j = f_0 · z_k + q · (Bz)_k-1.

This equation gives us an efficient algorithm: Given the previous output (Bz)_k-1 and the current input z_k, we can compute the current output (Bz)_k. The memory requirement of this algorithm is simply to store the previous output in a single memory buffer. 



Next suppose the sequence instead satisfied a linear recurrence of the form

    f_k = q_1 · f_k-1 + q_2 · f_k-2 + ⋯ + q_d · f_k-d.

As we explain next, this recurrence gives an algorithm where the memory requirement is to only store  memory buffers. 

The recurrence f_k = q f_k-1 implies the closed form f_k = q^k f_0.
Similarly, the recurrence in <Ref> implies a closed form expression f_k = u^T W^k v where W ∈ℝ^d × d is a matrix and u,v ∈ℝ^d are vectors.
This closed form is what we use in <Ref> for our <Ref>.
Specifically, we can extend <Ref> to this matrix-power closed form

    (Bz)_k    = ∑_j=0^k f_k-j· z_j 
        = ∑_j=0^k u^T W^k-j v · z_j
        
        = u^T S_k+1



for a suitable state vector S_k ∈^, stored in memory. Namely, we initialize S_0=0 and, at each iteration, our algorithm updates 
    S_k+1 = v · z_k + W S_k
 and then outputs (Bz)_k = u^T S_k+1.



We refer to the entries of S_k as the  buffers of our algorithm.[Recall that while we treat these as scalars here, in practical ML applications for example, each of these has size equal to the number of parameters of the model being trained, e.g. possibly > 10^9, so keeping  to a small constant is critical.] 

The updates to the buffers on each step are an arbitrary linear function of the previous step's buffers (WS_k) and the current input (v z_k), and the output on each step is an arbitrary linear combination  of the buffers (u^T S_k+1). Hence, we term this class “Buffered Linear Toeplitz matrices” (). We overload the acronym and write as a shorthand encompassing matrices, factorizations, and mechanisms. 



Unfortunately, the optimal factorization <cit.> does not satisfy a recurrence like <ref> and cannot be expressed as a . Hence, our approach is to approximate the optimal factorization using .



  
Designing via rational function approximation (<ref>):
If we view the Toeplitz sequence f_0, f_1, … as being defined by an ordinary generating function[We use the terms generating function, ordinary generating function, and Taylor series interchangeably; the generating function view of f emphasizes the sequence being generated and requires (only) a formal power series, while the Taylor series view emphasizes f is a real or complex function. As long as f is analytic in a non-empty open neighborhood of zero, no ambiguity is introduced by these two views, see e.g. Thm. 2.8 (Transfer principle) of <cit.>.]

 f as in <Ref>, it turns out that satisfying a linear recurrence as in <Ref> is equivalent to the function being rational (i.e., f(x) = p(x)/q(x) for polynomials p and q) of degree ≤.
(This equivalence is analogous to the fact that a real number is rational if and only if its decimal representation is repeating.)


This equivalence also suggests our first approach to designing : we need a low-degree rational approximation to the function f(x) = 1/√(1-x).

We appeal to known results in approximation theory. 
Specifically, it is known that the function x↦√(1-x) can be uniformly approximated on the unit complex disc {x ∈ℂ : |x| ≤ 1} with error η>0 by a rational function of degree d=O(log^2(1/η)) <cit.>. That is, there exists a rational function r of degree ≤ d, such that |r(x)-√(1-x)|≤η for all x∈ℂ with |x|≤ 1.

It “only” remains to translate this approximation guarantee back to the matrix factorization objective.
Parseval's identity allows us to bound the difference between the sequences of Taylor coefficients in terms of an integral:
Suppose f(x)=∑_k=0^∞ f_k x^k and f(x)=∑_k=0^∞f_k x^k. Then

    ∑_k=0^∞ |f_k-f_k|^2 = 1/2π∫_-π^π|f(x(θ))-f(x(θ))|^2  dθ   where   x(θ) = exp(√(-1)θ).

In our case, f(x)=1/√(1-x) corresponds to the optimal lower triangular Toeplitz factorization while f(x) is our rational approximation – either f(x)=1/r(x) or f(x)=r(x)/(1-x), where r(x) ≈√(1-x).
We are interested in a finite sum ∑_k=0^-1 |f_k-f_k|^2, rather than the infinite sum in <Ref> (which does not converge in our setting).
Thus we consider a weighted version of Parseval's identity where the integral goes around a circle in the complex plane centered at 0 with radius e^-1/2n. 
Note that we require the approximation guarantee to hold on the complex plane, not just the real line.
Once we have this bound on ∑_k=0^-1 |f_k-f_k|^2, the near optimality guarantee of <Ref> follows from the triangle inequality.

The time and space requirement of <Ref> is O(log^2 n). This dependence arises from the degree of the rational approximation.
A degree at least Ω(log^2(1/η)) is necessary for approximating √(1-x) with error ≤η even for real values x ∈ [-1,1] <cit.>.
Thus, unless we can exploit some slack in our analysis, it seems we require different techniques to bring the space down to O(log n).

The proof of <Ref> gives an explicit rational function approximation, which we can directly convert into a matrix factorization and feed into our algorithm. We term the coming from this approach s, with each choice of the degree  leading to a different (and successively better) approximation to the optimal Toeplitz factorization. 



  
Designing via direct optimization (<ref>):
Recall that our goal is to minimize (B,C)=(AC^-1,C).
In order to ensure computational efficiency, we restrict C to the class of .
Using approximation theory we can directly construct a near-optimal C. 
However, we can also approach this as an optimization problem and optimize the parameters that define C numerically.

This optimization is far from straightforward. The class of matrices can be parameterized in multiple ways. Converting the parameters for C into (AC^-1,C) is nontrivial to compute – much less optimize – when the size  is large.
Nevertheless, the class of is algebraically closed under multiplication and addition, and this structure combined with the connection to rational generating functions, provides powerful tools for reasoning about them.

We give a parameterization for the class of that allows us to efficiently compute  in time practically independent of the size , specifically , and also to compute gradients.
Being able to compute gradients allows us to optimize (AC^-1,C) numerically.



The first challenge is that we need to be able to effectively parameterize both C and A C^-1.
If C is a , then, as discussed above, its Toeplitz coefficients are given by the Taylor series for a rational generating function c(x) = p(x)/q(x) for polynomials p and q. Further, these coefficients have a simple closed-form expression given by <ref> (see also <ref>). 

Of course 1 / c(x) = q(x)/p(x) is also a rational function, and it in fact generates the Toeplitz coefficients of C (<ref>), and so C is also a . The lower-triangular matrix of ones, A, is trivially a , and so <ref> also implies B = A C is a .

In <ref>, <ref> shows that given a  C, we can (explicitly and in closed-form) derive the parameters of C, and hence a closed form for its Toeplitz coefficients.[More precisely, <ref> uses a parameterization of p and q which leads to closed forms for both C and C.]

Using these closed-form expressions for the Toeplitz coefficients enables us to directly compute (B, C) for B, C ∈^ in time , see <ref>.

This is immediately useful, for example in <ref> allowing us to plot the performance of our mechanisms for  up to 10^8 using only a few seconds of computer time.[We could have in fact scaled the  calculations for our mechanisms to arbitrary ; the bottleneck is in the computation of the exact  for <cit.>, which requires time ().]

More importantly, however, <ref> shows computing (B, C) is a differentiable function of the parameters of the B and C, and hence we can use a gradient-based optimization method such as L-BFGS to directly minimize  targeting a specific number of steps , where B and C are  with  buffers defined by 2 parameters. We term these mechanisms, and they perform extremely well in practice. For example, for = 10^7, a 4 has  that is 1.032× that of the optimal Toeplitz factorization <cit.>, and 7 is  1.001× optimal; for smaller  the results are even better, for example for = 10^4, 4 achieves 1.001× optimal. <ref> gives more complete results.





















  
Generalizations of the binary tree mechanism (<ref>):
The starting point for <Ref>  is the binary tree mechanism of <cit.>.
The binary tree mechanism can be viewed as a recursive construction of a matrix factorization. 
A recursion of depth ℓ yields a matrix factorization of size =2^ℓ and an algorithm running in time and space O(ℓ).
The binary tree mechanism does not produce a Toeplitz or square matrix factorization; the structure that it relies on for computational efficiency is sparsity.
The matrix factorization objective (B,C) for this construction is O(log) – that is, it is within a constant factor of optimal. Specifically, the binary tree mechanism is asymptotically a factor of π/log 2≈ 4.5 from optimal. This factor is significant in practice, as shown in <ref>. 

We combine the binary tree mechanism's recursive approach with <Ref> to get the best of both worlds – near-optimal constants and O(log n) space. This proves <Ref>.

We illustrate one step of the recursive construction using the following example for size n=6. We can decompose the 6 × 6 all-ones lower triangular matrix A^(6) into a sum of expressions involving a 2 × 2 all-ones lower triangular matrix A^(2) and a 3 × 3 all-ones lower triangular matrix A^(3):

    A^(6) =
        ([ 1 0 0 0 0 0; 1 1 0 0 0 0; 1 1 1 0 0 0; 1 1 1 1 0 0; 1 1 1 1 1 0; 1 1 1 1 1 1 ])
           =
        ([ 1 0 0 0 0 0; 1 1 0 0 0 0; 0 0 1 0 0 0; 0 0 1 1 0 0; 0 0 0 0 1 0; 0 0 0 0 1 1 ])
        + 
        ([ 0 0 0 0 0 0; 0 0 0 0 0 0; 1 1 0 0 0 0; 1 1 0 0 0 0; 1 1 1 1 0 0; 1 1 1 1 0 0;   ]) 
       = 
        ([ 1 0 0; 0 1 0; 0 0 1 ]) 
        ⊗([ 1 0; 1 1 ])
        +
        ([ 0 0 0; 1 0 0; 1 1 0 ]) 
        ⊗([ 1 1; 1 1 ]) 
       = I ⊗ A^(2) + ( S^(3) A^(3)) ⊗( 11^T ),

where ⊗ denotes the Kronecker product,[A key property of the Kronecker product is that (A · B) ⊗ (C · D) = (A ⊗ C) · (B ⊗ D).] I is the identity matrix, 1 is the all-ones vector, and S^(3) = ([ 0 0 0; 1 0 0; 0 1 0 ]) is a non-cyclic shift matrix.
<Ref> can be used to take factorizations of size 2 and 3 and combine them into a factorization of size 6. 
Namely, if A^(2)=B^(2)C^(2) and A^(3)=B^(3)C^(3), then 

    A^(6)   = I ⊗ A^(2) + ( S^(3) A^(3)) ⊗( 11^T ) <ref>
       = (I · I) ⊗ (B^(2)· C^(2)) + ( S^(3) B^(3)· C^(3)) ⊗( 1·1^T ) 
       = (I ⊗ B^(2)) · ( I ⊗ C^(2)) + ( S^(3) B^(3)⊗1) ·(C^(3)⊗1^T ) 
       = ( I ⊗ B^(2) |  S^(3) B^(3)⊗1)_=B^(6)·([ I ⊗ C^(2); C^(3)⊗1^T ])_=C^(6).

The factors B^(6) and C^(6) are non-square matrices represented as block matrices.
The reason we move to non-square matrices is that <Ref> decomposes A^(6) as a sum of two matrix products and we must re-express this as a single matrix product.

Roughly speaking, the binary tree mechanism corresponds to starting with a factorization of size 2 and applying a recursive step similar to the above ℓ-1 times to obtain a factorization of size 2^ℓ. 
Rather than starting with a factorization of size 2, we can start with a larger factorization of size _1 given by <Ref> and then repeat the recursive step above ℓ-1 times to obtain a factorization of size _1^ℓ.
Intuitively, by picking a larger factorization as the starting point we get closer to the optimal constant.
With careful analysis and the right choice of parameters, this yields <Ref>.
 
This recursive construction attains excellent asymptotics, but, for practical parameter regimes, we find that the approach is so good that we do not need the recursion.



 §.§ An Empirical Comparison of Mechanisms

In this section, we provide an empirical comparison of the primary DP mechanisms discussed, demonstrating their effectiveness in practical regimes.

<ref> compares mechanisms (and lower bounds) in terms of  from n=1 to n=10^6 iterations. The sub-optimality of the binary tree mechanism is immediately clear. 
This plot also shows that little is lost by the restriction from general matrix mechanisms to Toeplitz mechanisms, where <cit.> provide the optimal (but inefficient) construction. Our mechanisms essentially match this performance, while requiring time and memory (1) instead of ().

<ref> provides a detailed comparison of our mechanisms and the optimal Toeplitz mechanism. Several important conclusions can immediately be drawn: (1) For both and , increasing the degree (number of allowed memory buffers) increases performance. (2)  A larger number of steps  requires a higher number of buffers  for both of our approaches; this is expected and necessary, as shown by our theory, see <ref>. (3) A key point to emphasize is that each blue line for corresponds to a single mechanism (a fixed rational approximation); for , we compute an optimized matrix factorization of the given degree for each different . This specialization of the mechanism to the specific anticipated number of steps  is critical to the advantage enjoyed by this approach. For example, we see with only 5 buffers outperforms with 9 buffers across the full range of .

It is of course possible to run mechanisms for a different number of steps than the optimization targeted. <ref> (Left column) explores this. We construct three fixed factorizations, optimized for ∈{100, 1000, 10000}, and compare their performance (in terms of  relative to the optimal Toeplitz mechanism) across a range of steps, from 10^1 to 10^5. As expected, the mechanisms work best for the s for which they were optimized; however, the excess error is highly asymetric; a mechanism optimized for ^* will generally perform well for < ^* steps, but can quickly perform very badly when > ^*. This is expected when one considers that should intuitively be ensuring a good approximation of the optimal Toeplitz coefficients r_0, …, r_^*-1, but the approximation of the optimal Toeplitz coefficients for larger  can become arbitrarily bad.  We see this in <ref> (Middle column), where we compare the Toeplitz coefficients defining C and B to the optimal coefficients corresponding to the generating function 1/ √(1 - x). To emphasize this issue, we optimize for ^* = 100, and consider degree =2. The factorization provides a better approximation to the optimal coefficients for  up to 100 compared to 2 (which does not depend on ), and a generally worse approximation beyond that.  <ref> (Right column) shows that while both and correspond to “reasonable” approximations to √(1 -x ) (top), they distribute their errors in the approximation very differently (bottom).








§ BACKGROUND


Fix a dimension ∈ℕ and let 

    A = A^()( [ 1 0 0 ⋯ 0; 1 1 0 ⋯ 0; 1 1 1 ⋯ 0; ⋮ ⋮ ⋮ ⋱ ⋮; 1 1 1 ⋯ 1 ]) ∈{0,1}^×
 
be the all-ones lower triangular matrix given by A^()_i,j = 1  i ≥ j and A_i,j = 0  i < j for all i,j ∈ [].
Note that A is invertible and its inverse is a lower triangular bi-diagonal matrix:

    A^-1 = ( [  1  0  0  ⋯  0; -1  1  0  ⋯  0;  0 -1  1  ⋯  0;  ⋮  ⋱  ⋱  ⋱  ⋮;  0  ⋯  0 -1  1 ]) ∈{-1,0,1}^×.


Our goal is to compute a matrix factorization A=BC where B,C^T ∈ℝ^× that minimizes 

    (B,C) B_2 →∞·C_1 → 2√(max_i ∈ []∑_j ∈ [] B_i,j^2 )·√(max_j ∈ []∑_i ∈ [] C_i,j^2 ).



Simultaneously we want computational efficiency in the following sense. 
We need to generate samples from BZ ∈ℝ^× where Z ∈ℝ^× is a matrix with independent standard Gaussian entries. 
We want to do this in a streaming setting where we output one row at a time and the memory is limited. Ideally the memory should be linear in  and constant or logarithmic in .

We do not restrict B and C to be square matrices. 
Non-square factorizations may be advantageous from a computational perspective.
However, computational efficiency aside, we can assume without loss of generality that they are square (i.e., =) by taking a singular value decomposition (SVD) of the factors and discarding the rows/columns that do not correspond to a non-zero singular value.

It is also natural to restrict B and C to be lower triangular matrices, like A.
Lower triangular structure implies each intermediate output (Cx)_k is only a function of the inputs seen so far x_0,x_1,⋯,x_k and not future inputs x_k+1,⋯,x_n-1.
This is particularly valuable for the privacy analysis in the case where the input stream x_1, ⋯, x_k is chosen adaptively, such as in machine learning applications. In the adaptive setting, each input x_k may depend on the previous intermediate outputs (Cx)_0,(Cx)_1, ⋯, (Cx)_k-1. Lower triangular structure prevents a cyclic dependency.
Fortunately, we can assume the factors are lower triangular without loss of generality, see for example Prop 2.2 of <cit.>.



 §.§ Prior Work on Matrix Factorizations for Continual Counting


There has been a lot of work on factorizing the all-ones lower triangular matrix. 
However, most of that work optimizes the matrix factorization objective (B,C) with little regard for the computational efficiency of sampling Bz for z 𝒩(0,σ^2I).

The binary tree mechanism <cit.> gives a factorization which optimizes (B,C) up to constant factors and is efficiently computable. Prior to our work, this is the only known approach offering both efficient computation and some form of near-optimality.
The binary tree factorization can be expressed recursively as B^(1)=C^(1)=(1) and

    B^(2) = ([ B^()    0    0;    0 B^()   1⃗ ]) ∈{0,1}^2× (4-1)   and   
        C^(2) = ([ C^()    0;    0 C^(); 1⃗^T    0 ]) ∈{0,1}^(4-1) × 2.

Note that this only gives a factorization of A^() when  is a power of 2. By discarding rows/columns we can extend this to any . The objective value is given by 
    (B,C) = B^(2^ℓ)_2 →∞^2 = C^(2^ℓ)_1 → 2^2 = ℓ+1 = ⌈log/log 2⌉ + 1.

The recursive formulation naturally leads to an efficient sampling algorithm using O(ℓ) space.
Subsequent work has attempted to improve the constants <cit.>.

Observe that the binary tree factorization produces sparse matrices. Specifically, each row of B^(n) and column of C^(n) only has O(log n) nonzero entries.
Intuitively, this sparsity is what enables efficient computation; that is, (Bz)_k only depends on logarithmically many elements of z.
However, our experience suggests that sparsity requires non-square matrices and cannot produce near-optimal factorizations like . 

<cit.> give the following elegant explicit lower triangular Toeplitz factorization. 
Let f_0 = 1 and, for k ∈ℤ_≥ 1, let f_k = f_k-1· (1-1/2k) and f_-k=0. Equivalently, f_k = 4^-k2k  k≤1/√(π k) for k ≥ 1 <cit.>.
Define B,C ∈ℝ^× by B_i,j = C_i,j = f_i-j for all i,j ∈ [], which we denote B = C =((f_i)_i=0^n-1) or simply (f, n).
Then A=BC and[Unfortunately there appears to be an off-by-one error in the proof of the bound given by Theorem 1 of <cit.>; this has since been corrected.]

    (B,C) = B_2 →∞^2 = C_1 → 2^2 = ∑_k=0^-1 f_k^2 ≤ 1+log()+γ/π,

where γ≤ 0.57722 is the Euler-Mascheroni constant <cit.>.
This factorization turns out to be optimal among the class of lower triangular Toeplitz factorizations; we prove this in <ref>.

However, unlike the binary tree mechanism, these Toeplitz matrices are dense and we do not know how to efficiently generate noise according to this factorization in the streaming setting.

The γ_2 factorization norm is in fact a norm, i.e., it satisfies the triangle inequality <cit.>.
Furthermore, there is a dual characterization of the γ_2 factorization norm <cit.>. 

    γ_2(A)    inf{(B,C) : B,C ∈ℝ^×, A=BC }
       = sup{ P^1/2· A · Q^1/2_trace : 
        
        P,Q  non-negative diagonal matrices
     with 𝗍𝗋𝖺𝖼𝖾(P)=𝗍𝗋𝖺𝖼𝖾(Q)=1
        },

where M _trace = 𝗍𝗋𝖺𝖼𝖾(√(M^TM)) is the sum of the absolute singular values.

Using the triangle inequality, Mathias <cit.> gives a non-constructive upper bound 

    γ_2(A) ≤1/2 + 1/2∑_j=1^1/sin(π2j-1/2)≤ 1 + log ()/π.


On the other hand, using the dual characterization (<ref> with P=Q=1n I), <cit.> give a lower bound[Proposition 4.1 of <cit.> simply states an Ω(log) lower bound, but these sharper expressions can easily be extracted from the proof. <cit.> also gives a lower bound which differs from the upper bound in <ref> by less than an additive 1/2, but this seems slightly weaker than Equation <ref>.]

    γ_2(A) ≥1/2∑_j=1^1/sin(π2j-1/4+2)≥log(2+1)/π.

The upper bound in <ref> or <ref> and the lower bound in <ref> match up to a small additive constant. Thus the optimal value of the matrix factorization objective is γ_2(A) = log/π±𝖼𝗈𝗇𝗌𝗍𝖺𝗇𝗍.
Numerically, the maximum gap between the Toeplitz upper bound in <ref> and the lower bound in <ref> is less than 0.365.

Figure <ref> shows how these upper and lower bounds compare.
The binary tree mechanism has the advantage of computational efficiency, but as we can see, it is far from optimal for the objective.
The leading term for the binary tree is log/log 2, while the optimal leading term is log/π. Thus the binary tree is asymptotically suboptimal by a multiplicative factor of π/log 2 ≈ 4.5.



 §.§ Lower Triangular Toeplitz Factorizations versus General Factorizations
 



Our work focuses on lower triangular Toeplitz factorizations. In particular, <Ref> gives such a factorization and proves near-optimality with respect to this class. (Although we consider non-square factorizations in <Ref>, where we prove <Ref>.) 
This structure is essential to our algorithms in <Ref>.
However, it is natural to wonder how much we lose in terms of the objective (B,C) by restricting B and C to be lower triangular Toeplitz matrices.

In this subsection we discuss the gap between lower triangular Toeplitz factorizations and general factorizations.
We do not have an explicit construction for the optimal general matrix factorization or even a formula for optimal value γ_2(A). 
For lower triangular Toeplitz matrix factorizations the best factorization we have is that of <cit.>.
We show that this is in fact optimal among the class of lower triangular Toeplitz factorizations.
This is – to the best of our knowledge – a novel result which may be of independent interest.

We begin with a lemma stating some basic properties of the factorization of <cit.> and we give a proof for completeness.


    Define a sequence f_0, f_1, ⋯ by f_0=1 and f_k = (1-1/2k)f_k-1 for all k ≥ 1.
    Then, for all k ≥ 1, we have 
    
    f_k = 4^-k2k  k∈[ 1/√(π(k+1)), 1/√(π k)].

    For all integers n ≥ 0, we have 
    
    ∑_k=0^n f_k f_n-k = 1

    and
    
    γ+log(n) - 1/π≤∑_k=1^∞ f_k^2 ≤γ + log(n)/π,

    where 0.57721 ≤γ≤ 0.57722 is the Euler-Mascheroni constant.
    Furthermore, for all x∈ℂ with |x|<1, we have
    
    ∑_k=0^∞ f_k x^k = 1/√(1-x).



The first part of <ref> can be shown by induction:
For k ≥ 1,

    4^-k2k  k = 1/4· 4^-(k-1)2k-2  k-12k (2k-1)/k^2 = 2k(2k-1)/4k^2 f_k-1 = ( 1 - 1/2k) f_k-1 = f_k.

Next consider the derivatives of f(x) 1/√(1-x): For k ≥ 1, we have 
    f^(k)(x) = (1-x)^-k-1/2∏_ℓ=0^k-1(ℓ+1/2).

Thus 
    f^(k)(0)/k! = ∏_ℓ=0^k-11/ℓ+1(ℓ+1/2) = ∏_ℓ=0^k-1(2ℓ+1)/2(ℓ+1)(2ℓ+2)/2(ℓ+1) = (2k)!/(k!)^2 2^2k = 4^-k2k  k = f_k.

By Taylor's theorem, for all x ∈ℂ with |x|<1, 
    1/√(1-x) = f(x) = ∑_k=0^∞f^(k)(0)/k! x^k = ∑_k=0^∞ f_k x^k,
 as required to prove <ref>.
Furthermore, for all x ∈ℂ with |x|<1, we have 
    ∑_n=0^∞ x^n = 1/1-x = f(x)^2 = ( ∑_k=0^∞ f_k x^k )^2 = ∑_n=0^∞ x^n ∑_k=0^n f_k f_n-k.

Matching coefficients proves <ref>

To prove the second part of <ref>, we use standard bounds on the central binomial coefficient <cit.>: For all integers k ≥ 1, 
    4^k/√(π(k+1))≤2k  k≤4^k/√(π k),
 whence 1/√(π (k+1))≤ f_k ≤1/√(π k) . 
Next we use standard bounds on the harmonic numbers to prove <ref>: For all integers n ≥ 1,

    log(n) + γ + 1/2(n+1)≤∑_k=1^n 1/k≤log(n) + γ + 1/2n

where 0.57721 ≤γ≤ 0.57722 is the Euler-Mascheroni constant <cit.>.
It follows that, for all n ≥ 2,

    ∑_k=1^n-1 f_k^2 ≤∑_k=1^n-11/π k≤γ + log(n-1) + 1/2(n-1)/π≤γ + log(n)/π.

Similarly, for all n ≥ 2

    ∑_k=1^n-1 f_k^2 ≥∑_k=1^n-11/π(k+1) =  ∑_k=2^n 1/π k≥γ + log(n) + 1/(2n+2)-1/π.



Most importantly, we show that the factorization of <cit.> is optimal among the class of lower triangular Toeplitz factorizations:


Let ,∈ℝ^n be any pair of two vectors. Consider the following optimization problem:

    [                             min_,                                3l;                              s.t. ∀ k  <n∑_0≤ i,j<n
     i+j=kb_ic_j=1.;                                   ]

The minimum is achieved at ==f, where f constitute the coefficients derived in <cit.>, i.e., the coefficients given by f_0 = 1, f_i=1-1/2if_i-1 for all i ≥ 0.



The proof proceeds in 3 steps: We first reformulate the optimization as a quadratically constrained quadratic program, then derive a natural Lagrangian relaxation of it, and finally construct a primal dual solution to the Lagrangian relaxation that is globally optimal with the primal solution ==f.

We begin by noting the following properties of f:

    
  * f_i > 0   ∀ i ∈{0, …, n-1}.
    
  * Let f denote the lower triangular Toeplitz matrix with first column equal to f. Then we have ff=A.
    
  * ∑_i=0^inff_i x_i = 1/√(1-x)  ∀ x ∈ [0, 1].

These properties were established in <cit.> and also are shown in <ref>.



  
Reformulation as a quadratically constrained quadratic program:
We begin by observing that the optimization problem can be rewritten as a quadratic optimization as follows:

    = min_ν > 0ν/2b^2 + 1/2νc^2 = min_ν > 01/2ν b^2 + 1/2c/ν^2  .

Furthermore, for any b, c that satisfy the constraints of <ref>, ν b, c/ν also satisfy the constraints (since the constraints only involve products of coordinates of b, c which remain preserved if b is multiplied by ν and c by 1/ν and are hence feasible. If S denotes the feasible set of <ref>, we have 

    min_,∈ S = min_,∈ Smin_ν > 01/2ν b^2 + 1/2c/ν^2 = min_, ∈ S1/2^2 + 1/2^2.


Thus, we can simply write the overall optimization problem as 

    [                        min_,                 3l1/2^2 + ^2;                         s.t. ∀ k  <n∑_0≤ i ≤ kb_ic_k-i=1.;                              ]





  
Lagrangian relaxation:
We write the Lagrangian relaxation of this optimization problem introducing dual variables λ_0, ⋯, λ_n-1 corresponding to the constraints:

    Lb,c,λ = 1/2^2 + ^2    - ∑_k=0^n-1λ_k ∑_0≤ i ≤ kb_ic_k-i-1
       =1/2^2 + ^2 - Γλ + λ
        = 1/2[ b; c ][    I - Γλ;  -Γλ    I ][ b; c ] + λ ,

where Γλ is the symmetric matrix

    Γλ[   λ_0   λ_1   λ_2     … λ_n-1;   λ_1   λ_2   λ_3     …     0;   λ_2   λ_3   λ_4     …     0;     ⋮     ⋮     ⋮     ⋮     ⋮; λ_n-1     0     0     …     0 ].

The constrained optimization problem in <ref> is equivalent to the unconstrained min-max problem min_b,cmax_λ L(b,c,λ).
We will set b=c=f and exhibit a setting λ=λ^* such that ∇_b,c,λL(b,c,λ) = 0 and L(b,c,λ) is jointly convex in b,c and concave in λ. It follows that b=c=f is the optimal solution to the original problem, as required.
Since b=c=f is feasible, ∇_λ L(b,c,λ) = 0.
The Lagrangian is linear in λ and hence trivially concave.
The Lagrangian is convex in , if and only if 

    [    I - Γλ;  -Γλ    I ]≽ 0

or equivalently (by Schur complements) if Γλ_⋆≤ 1 where Q_⋆sup{Qx_2 : x_2 ≤ 1 } denotes the operator norm. 



  
Construction of primal-dual optimal solution:
Now it only remains to compute λ=λ^* such that ∇_b,c L(b,c,λ) = 0 for b=c=f and Γλ^2_⋆≤ 1.

For all i ∈ [n], we have ∂/∂ b_i L(b,c,λ) = b_i - ∑_k=i^n-1λ_k c_k-i and ∂/∂ c_i L(b,c,λ) = c_i - ∑_k=i^n-1λ_k b_k-i. Equivalently, ∇_b L(b,c,λ) = b - Γ(λ) c and ∇_c L(b,c,λ) = c - Γ(λ) b. 
Thus, in order to ensure ∇_b,c L(b,c,λ) = 0  when b=c=f, it suffices to set λ = λ^⋆, where λ^⋆ is chosen so as to solve 

    f = Γλ^⋆f  Pf = P Γλ^⋆f [ f_n-1; f_n-2;     ⋮;   f_0 ] = [ λ^⋆_n-1       0       0       …       0; λ^⋆_n-2 λ^⋆_n-1       0       …       0; λ^⋆_n-3 λ^⋆_n-2 λ^⋆_n-1       …       0;       ⋮       ⋮       ⋮       ⋮       ⋮;   λ^⋆_0   λ^⋆_1   λ^⋆_2       … λ^⋆_n-1 ][   f_0;   f_1;     ⋮; f_n-1 ]
    [ f_n-1; f_n-2;     ⋮;   f_0 ] = [   f_0     0     0     …     0;   f_1   f_0     0     …     0;   f_2   f_1   f_0     …     0;     ⋮     ⋮     ⋮     ⋮     ⋮; f_n-1 f_n-2 f_n-3     …   f_0 ][ λ^⋆_n-1; λ^⋆_n-2;       ⋮;   λ^⋆_0 ] ,

where P is the permutation matrix such that Px is the reversal of the vector x. Thus, we can solve for λ^⋆ to obtain

    P λ^⋆ = [ λ^⋆_n-1; λ^⋆_n-2;       ⋮;   λ^⋆_0 ] = [         f_0           0           0           …           0;     f_1-f_0         f_0           0           …           0;     f_2-f_1   f_1 - f_0         f_0           …           0;           ⋮           ⋮           ⋮           ⋮           ⋮; f_n-1-f_n-2 f_n-2-f_n-3 f_n-3-f_n-4           …         f_0 ][ f_n-1; f_n-2;     ⋮;   f_0 ],

where we used the fact that f is a feasible solution:

    [   f_0     0     0     …     0;   f_1   f_0     0     …     0;   f_2   f_1   f_0     …     0;     ⋮     ⋮     ⋮     ⋮     ⋮; f_n-1 f_n-2 f_n-3     …   f_0 ][   f_0     0     0     …     0;   f_1   f_0     0     …     0;   f_2   f_1   f_0     …     0;     ⋮     ⋮     ⋮     ⋮     ⋮; f_n-1 f_n-2 f_n-3     …   f_0 ] = [ 1 0 0 … 0; 1 1 0 … 0; 1 1 1 … 0; ⋮ ⋮ ⋮ ⋮ ⋮; 1 1 1 … 1 ]
    [   f_0     0     0     …     0;   f_1   f_0     0     …     0;   f_2   f_1   f_0     …     0;     ⋮     ⋮     ⋮     ⋮     ⋮; f_n-1 f_n-2 f_n-3     …   f_0 ]^-1 = [ 1 0 0 … 0; 1 1 0 … 0; 1 1 1 … 0; ⋮ ⋮ ⋮ ⋮ ⋮; 1 1 1 … 1 ]^-1[   f_0     0     0     …     0;   f_1   f_0     0     …     0;   f_2   f_1   f_0     …     0;     ⋮     ⋮     ⋮     ⋮     ⋮; f_n-1 f_n-2 f_n-3     …   f_0 ]
    [   f_0     0     0     …     0;   f_1   f_0     0     …     0;   f_2   f_1   f_0     …     0;     ⋮     ⋮     ⋮     ⋮     ⋮; f_n-1 f_n-2 f_n-3     …   f_0 ]^-1 = [  1  0  0  …  0; -1  1  0  …  0;  0 -1  1  …  0;  ⋮  ⋮  ⋮  ⋮  ⋮;  0  0  0  …  1 ][   f_0     0     0     …     0;   f_1   f_0     0     …     0;   f_2   f_1   f_0     …     0;     ⋮     ⋮     ⋮     ⋮     ⋮; f_n-1 f_n-2 f_n-3     …   f_0 ].


From the definition f_0=1 and f_i = f_i-1(1-1/2i), for all i ∈ [n]={0,⋯,n-1}, we have

    λ^⋆_n-1-i   = f_0 f_n-1-i - ∑_j=1^i f_j-f_j-1f_n-1-i+j
       = f_n-1-i - ∑_j=1^i f_j-11/2jf_n-1-i+j
       =f_n-1-i1-∑_j=1^if_j-11/2jf_n-1-i+j/f_n-1-i.


Since f_n-1-i+j/f_n-1-i = ∏_ℓ=1^j 1-1/2(n-i+ℓ)≤ 1, we have

    λ^⋆_n-1-i≥ f_n-1-i1-∑_j=1^if_j-11/2j > f_n-1-i1-∑_j=1^∞  f_j-11/2j = f_n-1-i1-∑_j=0^∞f_j/2(j+1).

Since ∑_j=0^∞ f_j x^j = 1/√(1-x) for all x ∈ [0 , 1), integrating both sides between the limits 0 and t, we have

    ∑_j=0^∞f_j/j+1 t^j+1=21-√(1-t).

Taking limits as t → 1, we obtain

    ∑_j=0^∞f_j/j+1=2.

Thus, from <ref>, we can conclude  that λ^⋆_n-1-i > 0 for all i ∈ [n] and hence Γλ^⋆^2 is a matrix with all entries (Γ(λ^⋆)^2)_i,j = ∑_k=0^n-1-max{i,j}λ^⋆_i+kλ^⋆_j+k strictly positive. Hence, by the Perron Frobenius theorem, we have that Γλ^⋆^2 has a unique eigenvector with all coordinates positive, and this corresponds to the maximum eigenvalue. Since f is an eigenvector of Γλ^⋆^2 by construction, it must be the unique leading eigenvector corresponding to the eigenvalue 1. Hence Γλ^⋆^2_* ≤ 1.

Thus, L, , λ is convex in ,, since Γλ^2_⋆≤ 1 Γλ_⋆≤ 1.

Thus, we have ^⋆, f, λ^⋆ is a stationary point of the Lagrangian that is convex in , and concave in λ, hence it is the global optimum. 




<Ref> fully characterizes the optimal lower triangular Toeplitz factorization.
We now compare the value of the matrix factorization objective for lower triangular Toeplitz factorizations versus general factorizations.
For convenience, we define a notation for these values:


Let ∈{0,1}^n× n be the lower triangular matrix of all ones.
Recall that C ∈ℝ^n × n is lower triangular Toeplitz if C_i,j=0 for all i < j and C_i,j = c_i-j for all i ≥ j, where c ∈ℝ^n is a vector.
Define 

    𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(n)    inf{(B,C) : B,C ∈ℝ^n × n lower triangular Toeplitz, BC=A }, 
    𝖮𝗉𝗍(n)    inf{(B,C) : B,C ∈ℝ^n × n , BC=A },

where (B,C) B_2→∞C_1 → 2 = √(max_i ∑_j B_i,j^2)√(max_j ∑_i C_i,j^2).

By <Ref> and <Ref>, for all integers n ≥ 1, we have

    1 + γ-1 + log(n)/π≤𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(n) = 1 + ∑_k=1^n-1( 2^-2k2k  k)^2 ≤ 1 + γ+log(n)/π.

The lower bound of <cit.> in <ref> states the optimal over all factorizations is 
    𝖮𝗉𝗍(n) = γ_2(A^(n)) ≥ 1 + log(2n+1)/π.

Combining these bounds gives a bound on the gap between the optimal general factorization and the optimal lower triangular Toeplitz factorization.

    The gap between the lower triangular Toeplitz and general optimal factorizations is 

    𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(n)-𝖮𝗉𝗍(n) ≤(1 + γ+log(n)/π)-(1 + log(2n+1)/π) ≤γ+log(2)/π≤ 0.405.


Numerically, we have 𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(n)-𝖮𝗉𝗍(n) ≤ 0.365 for all n ≥ 1.
<Ref> shows numerically how the various known upper and lower bounds compare for differing values of n.

We conclude this subsection by remarking that, without the Toeplitz constraint, we can assume that the factors are lower triangular without loss of generality.
This is a consequence of the fact that two factorizations A=BC and A=B̂Ĉ are functionally equivalent as long as BB^T = B̂B̂^T. To see why this fact is true, note that for z 𝒩(0,σ^2I), the distribution of Bz is 𝒩(0,σ^2 BB^T). So if BB^T = B̂B̂^T, then Bz and B̂ z have the exact same distribution. Note that this fact is specific to the Gaussian distribution.
Given this fact, we can take an arbitrary factorization A=BC and perform a Cholesky decomposition of BB^T to obtain a lower triangular B̂ satisfying BB^T = B̂B̂^T.
However, if we restrict to Toeplitz matrices, then we cannot assume they are lower triangular without loss of generality, since the Cholesky decomposition may not preserve this structure.




 §.§ Applications of Matrix Factorizations in Machine Learning
 

In the recent years, Matrix Factorization (MF) have found extensive use both in obtaining optimal private regret/stochastic convex optimization (SCO) guarantees <cit.>, and in practical deployments <cit.>. The main observation that maps MF to the (convex) optimization setup is that (unconstrained) stochastic graident descent update with constant learning rate can be viewed as computing prefix sums over adaptively chosen gradient vectors[There are more general mappings for constrained optimization, and adaptive optimizers like momentum methods, but a detailed discussion on those is tangential to this work.]. Differentially Private Follow-the-regularized-leader (DP-FTRL) <cit.> is an optimization algorithm that arises naturally from this mapping. DP-FTRL has found extensive usage in production grade deployment of DP learning. For example, all next-word-prediction models on Android Gboard is trained with DP-FTRL <cit.>. 

For a given set of ℓ_2-Lipschitz loss functions {f_1(θ),…,f_(θ)} (where each loss function f_i can be thought of as the optimization loss on a disjoint minibatch of data samples), the simplest version of DP-FTRL can be described as follows:



    θ_t← -η∑_i=0^t-1∇ f_i(θ_i) + ẑ_i.

In <ref>, η is the learning rate, and {ẑ_t} are the correlated Gaussian noise added to ensure that the computation of the sequence of {θ_t}'s are differentially private. Mapping to the MF view in <ref>, the adaptively chosen data set is set to x={∇ f_t(θ_t): t∈[]}, and the noise sequence {ẑ_t} is set to B z, where z ←(0,σ^2I). DP-FTRL captures a large class of noise addition mechanisms (parameterized by the specific factorization scheme A=BC chosen). In particular, it captures DP-SGD <cit.> (another heavily used private learning algorithm) as a special case <cit.>, and has much better privacy/utility trade-offs in general <cit.>. <cit.> showed that for obtaining the best regret/SCO guarantees from DP-FTRL style algorithms, it suffices to optimize for (B, C) to instantiate the noise mechanism in <ref>.

One major challenge in large deployments of DP-FTRL is that noise generation for an arbitrary matrix mechanism can be prohibitively expensive, as discussed in <ref>. Our work provides a practical approach that only requires a constant factor blow-up (usually at most 4-6 times) in memory as compared to that of DP-SGD, and at privacy/utility trade-off comparable to DP-FTRL with an optimal matrix-factorization mechanism. (DP-SGD adds independent Gaussian noise to the gradients at each state update.)






 §.§ Other Related Work




  
Beyond additive noise:
Our algorithms are based on adding Gaussian noise (with a carefully chosen covariance). Since the noise is additive, the error is independent of the number of data points.[The error of our algorithms depends on the number of iterations . This is usually closely related to the number of data points |x|, but, in general, these two quantities can be unrelated.] The expected maximum error is 
    ∀ x    (x)-Ax_∞ = max_k |((x)-Ax)_k|≤ O(log(n))^3/2,

while the maximum expected error is 

    ∀ x    max_k |((x)-Ax)_k|≤ O(log(n)).

The additional √(log(n)) term in <Ref> compared to <Ref> comes from taking a union bound over n Gaussians.
Note that, for simplicity, we suppress the privacy parameters; to achieve (ε,δ)-differential privacy, the error above scales with 1/εlog(1/δ); to achieve ρ-zCDP <cit.>, the error above scales with 1/√(ρ).

The bound in <ref> on the maximum expected error is tight up to constants in the streaming setting – in fact the lower bound holds for average expected error: 

    ∀ ∃ x    1/n∑_k |((x)-Ax)_k|≥Ω(log(n)),

assuming  is differentially private <cit.>.
However, the bound in <ref> on expected maximum error is not tight if we move beyond additive noise mechanisms <cit.>.
<cit.> show that we can obtain an improved error guarantee:[<cit.> stated this result for pure differential privacy, with log(|x|)^2 instead of log(|x|)^3/2.]

    ∃ ∀ x    (x)-Ax_∞≤ O(log(n) + log(|x|)^3/2),

where |x| is the number of data points.[In most of the differential privacy literature n is the number of data points / people, but we use n for the number of iterations; |x| should be read as the size of the dataset x.]
The high level idea is to “freeze” the output of  until the partial sum has changed significantly and only then update the output.
This works because deciding whether the partial sum has changed significantly is easier than evaluating it. Specifically, we can use the so-called “sparse vector” technique to decide when to update.
This algorithm is a reduction – when it does decide to update, it uses additive noise – the “win” is that we need to use additive noise for ≤ |x| updates, rather than n updates. Thus our improvements for additive noise mechanisms can be applied in this setting.

<ref> is only an improvement over <ref> in the sparse setting (|x| ≪ n). Machine learning applications are typically dense (|x| ≥ n).



  
Offline setting:
Our work is in the streaming setting, where we must output an approximation to each partial sum as soon as it is complete.
It is also natural to consider the offline setting, where we are given the entire input at once.
In the offline setting, the memory constraints that we focus on are not as relevant.

The offline version of continual counting is known as threshold queries or quantiles and this task has been studied extensively <cit.>.
The offline setting permits more sophisticated algorithms, which achieve asymptotically better error guarantees.
<cit.> give an (ε,δ)-differentially private algorithm  with error

    ∀ x    (x)-Ax_∞≤O(1/εlog^2(1/δ) log^*(n))

Note that offline algorithms are not applicable in the machine learning applications discussed in <Ref>. 



  
Pure differential privacy:
The early work on continual counting <cit.> works under pure (ε,0)-differential privacy, which rules out Gaussian noise.
Instead Laplace noise can be used. In this case, we must bound the L_1 sensitivity. So the objective for the matrix factorization is to minimize B_2→∞·C_1 → 1, where C_1 → 1 = max_j ∑_i |C_i,j| is the maximum 1-norm of a column. 
This subtly changes the problem. The binary tree factorization still works well with this objective. However, Toeplitz factorizations do not work as well. Intuitively, the binary tree factorization relies on sparsity which controls both the 1- and 2-norms, whereas Toeplitz factorizations produce dense matrices for which the the 1- and 2-norms differ substantially.



  
Factoring other matrices:
We exclusively look at factorizing the all-ones lower triangular matrix. This is a natural and well-motivated problem.
However, there is also work on factorizing other matrices.

Most closely related to the all-ones lower triangular matrix A, <cit.> gives exact expressions for γ_2(A+A^T) and γ_2(A+A^T+I).

There is a different line of work <cit.> considering factorizations of generic workload matrices, which have a completely different structure than the matrices we consider.



 §.§ Notation Summary

Before proceeding, we briefly summarize the key symbols and notation used throughout:






§ EFFICIENTLY-SAMPLEABLE FACTORIZATION VIA RATIONAL GENERATING FUNCTIONS
 

In this section we describe our main algorithmic tools. 
First (in <ref>) we present the view of lower triangular Toeplitz matrices in terms of generating functions or sequences. This is a convenient mathematical formalism for analysis <cit.>.
Second (in <ref>) we discuss the special case of rational generating functions or, equivalently, constant-recurrent sequences. Our results focus on this special case, as the additional structure is useful for our algorithms; this structure is where we depart from the prior literature.
Third (in <ref>) we present our algorithm for sampling from lower triangular Toeplitz matrices with rational generating functions (that is, the multiplication algorithm for matrices).
In subsequent sections (<ref>, <ref>) we instantiate these rational generating functions.



 §.§ Sequences, Lower Triangular Toeplitz Matrices, & Generating Functions
 

We begin by describing the generating function view of lower triangular Toeplitz matrices used by <cit.>. 

Our analysis moves fluidly between three different views of the same mathematical object:

    
  * The sequence f_k = f_0, f_1, f_2, … with f_k ∈ for each index k.
    
  * The ordinary generating function of the sequence, f: →,
      
    f(x) = ∑_k=0^∞ f_k x^k = ∑_k=0^∞f^(k)(0)/k! x^k.

    (It will be clear from context if f refers to the sequence entries or the generating function).
    
  * The family of  lower-triangular Toeplitz matrices M(f,) ∈ℝ^× for ∈ that are generated by f as
    
    ∀ i,j ∈ []           
        M(f,)_i,j{[     f_i-j if  i ≥ j;         0   if  i<j ]}
 
    When the specific  is unimportant or clear from context, we write M(f) to refer to any matrix in this family (or equivalently, the infinite-dimensional linear operator).








As a canonical example, we have the all-ones sequence,


    g_k = 1, 1, 1, 1, …, 
      
     g(x) = 1/1-x = ∑_k=0^∞ x^k, 
       
    A^(4) M(g, 4) = [ 1 0 0 0; 1 1 0 0; 1 1 1 0; 1 1 1 1;   ]






It is straightforward to see that addition (under the usual definitions) is equivalent across all three representations, and (slightly less obviously) the same fact holds for suitable definitions of multiplication. We denote the Cauchy product or convolution * for sequences f and g by

    f_k * g_k = h_k  where   
        h_k ∑_i=0^k f_i · g_k-i,

and write h = f * g as a shorthand.

Then, the mathematical structure summarized in the following lemma is key to our analysis:

  Let f, g, and h be ordinary generating functions with corresponding sequences. Then the Cauchy product of sequences, multiplication of the generating functions, and matrix multiplication are all equivalent. That is,
         
    (h = f * g ) 
               ⟺  ( h(x) = f(x) g(x) )
               ⟺  ( ∀   M(h,) = M(f,) · M(g,) ).

      Similarly the usual definitions of addition are equivalent:
      
    ( ∀ i    h_i = f_i + g_i )
               ⟺  ( h(x) = f(x) + g(x) )
               ⟺  ( ∀   M(h,) = M(f,) + M(g,) ).


























    The results for sequences and their ordinary generating functions are well known, see e.g. <cit.>, and the results for addition are immediate. Hence, it is sufficient to show h(x) = f(x) g(x) iff 
    ∀∈, M(h,) = M(f,) · M(g,). 

    Fix i,j ∈ []. 
    If i<j, then M(h,)_i,j = 0 and (M(f,) · M(g,))_i,j = ∑_k ∈ [] M(f,)_i,k· M(g,)_k,j = ∑_k 0, since, for all k ∈ [], either i<k, whence M(f,)_i,k=0, or k<j, whence M(g,)_k,j=0, (or both). 
    Thus we can assume i ≥ j.

    By the product rule and induction, we have the derivative h^(ℓ)(x) = ∑_k = 0^ℓℓ k f^(ℓ-k)(x) · g^(k)(x) for all ℓ≥ 0 and all applicable x.
    Thus
    
    M(h,)_i,j   = h^(i-j)(0)/(i-j)!
       = 1/(i-j)!∑_k = 0^i-ji-j  k f^(i-j-k)(0) · g^(k)(0) 
       = ∑_k = 0^i-jf^(i-j-k)(0)/(i-j-k)!·g^(k)(0)/k!
       = ∑_k' = j^if^(i-k')(0)/(i-k')!·g^(k'-j)(0)/(k'-j)!k'=k+j
       = ∑_k' = j^i M(f,)_i,k'· M(g,)_k',j
       = ∑_k' ∈ [] M(f,)_i,k'· M(g,)_k',ji<k'  M(f,)_i,k'=0, k'<j  M(g,)_k',j=0
       = ( M(f,) · M(g,) )_i,j.





Given <ref> and the above result for multiplication, in order to factorize =, it is natural to factorize by taking the square root of the generating function: Let f(x) = 1/√(1-x). Then f(x) · f(x) = 1/1-x = g(x) and, hence, M(f,) · M(f,) = M(g,) = A^().
Indeed, this is the factorization given by <cit.>, which is optimal among the class of lower triangular Toeplitz matrices (per <Ref>).

Given any generating function r(x), we can obtain a valid factorization = where = M(b,) for b(x)=r(x)/(1-x) and = M(c,) for c(x) = 1/r(x).
Our approach is to choose a generating function such that r(x) ≈√(1-x) but which also permits us to design an efficient algorithm.



 §.§ Rational Generating Functions and Constant-Recurrent Sequences
 

We will focus our attention on rational generating functions of the form

    r(x) = p(x)/q(x) = ∑_i=0^d-1 p_i x^i/1+∑_j=1^d q_j x^j .

These have convenient algorithmic properties. In this subsection we discuss several equivalent representations of rational functions.
As a starting point, a rational function of degree d can be represented as a ratio of polynomials, where the numerator p(x) has degree <d and the denominator has degree ≤ d.
Note that we normalize the denominator so that q(0)=q_0=1.[We can rescale both the numerator and denominator by a constant to make this assumption true. Thus this assumption is without loss of generality, unless the rational function has a pole at 0 – but we do not consider such functions, as they are not valid generating functions.]

<Ref> is a convenient mathematical representation, but for algorithmic purposes we need the sequence r_k of Taylor series coefficients

    r(x) = ∑_k=0^∞ r_k x^k.

Thus our first task is to map from <Ref> to <Ref> with a convenient representation that is amenable to efficient computation. 

The generating function r(x) being a rational function is equivalent to the sequence r_k being a constant-recurrent sequence (a.k.a. linear-recursive or C-finite sequence); see for example <cit.>.
The sequence terms can be computed by taking powers of a matrix; this is the representation that we will use in our algorithm.
We summarize this result in <ref> and provide a proof for completeness.



    Let 
    
    r(x) = p(x)/q(x) = ∑_i=0^d-1 p_i x^i/1+∑_j=1^d q_j x^j

    be a rational function of degree ≤ d.    
    
    As in <ref>, let 
    r(x) = ∑_k=0^∞ r_k x^k.
 That is, r(x) generates the sequence r_k. 
    Then, for all k ≥ d, this sequence satisfies the recurrence
    
    r_k = - ∑_j=1^d q_j r_k-j.

    And, for 0 ≤ k < d, we have r_k = p_k - ∑_j=1^k q_j r_k-j. 
    Furthermore, for all k, we have 
    
    r_k = u^T W^k v,

    where
    
    u ([ 1; 0; 0; 0; ⋮; 0; 0 ]) ∈ℝ^d × 1,     
            v ([   p_0;   p_1;   p_2;   p_3;     ⋮; p_d-2; p_d-1 ]) ∈ℝ^d × 1,

    and
    
    W ([   -q_1      1      0      0      ⋯      0      0;   -q_2      0      1      0      ⋯      0      0;   -q_3      0      0      1      ⋯      0      0;   -q_4      0      0      0      ⋯      0      0;      ⋮      ⋮      ⋮      ⋮      ⋱      ⋮      ⋮; -q_d-1      0      0      0      ⋯      0      1;   -q_d      0      0      0      ⋯      0      0 ]) ∈ℝ^d × d.



 
    We have 
    
    r(x) = p(x)/q(x) = ∑_i=0^d-1 p_i x^i/1+∑_j=1^d q_j x^j = ∑_k=0^∞ r_k x^k ,

    which rearranges to 
    
    p(x) = ∑_i=0^d-1 p_i x^i = q(x) · r(x) = ( 1+∑_j=1^d q_j x^j ) ·( ∑_k=0^∞ r_k x^k ) = ∑_k=0^∞( r_k + ∑_j=1^min{d,k} q_j r_k-j) · x^k.

    Matching coefficients gives p_k = r_k + ∑_j=1^k q_j r_k-j for 1 ≤ k ≤ d-1 and 0 = r_k + ∑_j=1^d q_j r_k-j for k ≥ d, which rearranges to give the recurrence.

    Next we turn to the matrix representation.
    Define p^(0)(x) = p(x) – i.e., p^(0)_i = p_i for 0 ≤ i ≤ d-1.
    For all k ≥ 1, let p_j^(k) p_j+1^(k-1) - p_0^(k-1) q_j+1 for 0 ≤ j ≤ d-2 and let p_d-1^(k) - p_0^(k-1) q_d.
    We can check that these values satisfy
    
    ∑_i=0^d-1 p_i^(k-1) x^i/1+∑_j=1^d q_j x^j = p_0^(k-1) + x ·∑_i=0^d-1p_i^(k) x^i/1+∑_j=1^d q_j x^j.

    The left hand side of <Ref> with k=1 is simply r(x); from this we can conclude that r_0 = p_0^(0) = p_0.
    Expanding <Ref> for k=1 and k=2 gives r_1 = p_0^(1). Specifically, we have
    
    r(x) = ∑_i=0^d-1 p_i^(0) x^i/1+∑_j=1^d q_j x^j = p_0^(0) + x ·∑_i=0^d-1p_i^(1) x^i/1+∑_j=1^d q_j x^j = p_0^(0) + x ·( p_0^(1) + x ·∑_i=0^d-1p_i^(2) x^i/1+∑_j=1^d q_j x^j).

    We can iteratively repeat this expansion to extract the entire Taylor series.
    (This algorithm is known as polynomial long division.)
    That is, we can show by induction that, for all k, we have
    
    r(x) = ∑_ℓ=0^k-1 p_0^(ℓ) x^ℓ + x^k ·∑_i=0^d-1p_i^(k) x^i/1+∑_j=1^d q_j x^j.

    Thus r_k = p_0^(k) for all k.
    
    Now we can write the recursive definition of p^(k)_j in matrix form: For all k ≥ 1, we have
    
    ([   p_0^(k);   p_1^(k);         ⋮; p_d-1^(k) ])
             = W ·([   p_0^(k-1);   p_1^(k-1);           ⋮; p_d-1^(k-1) ]),

    where W ∈ℝ^d × d is given in Equation <ref>.
    Given this, induction shows that, for all k ≥ 0,
    
    ([   p_0^(k);   p_1^(k);         ⋮; p_d-1^(k) ])
             = W^k ·([   p_0^(0);   p_1^(0);         ⋮; p_d-1^(0) ])
             = W^k · v

    and, hence, r_k = p_0^(k) = u^T · W^k · v, as required.


We make some remarks about <Ref>:

    
  * <ref> assumes (p)<d and (q)≤ d, but it allows for the possibility that some of the coefficients p_j,q_j (including q_d) may be zero. Thus by padding, we can accommodate arbitrary degree of both the numerator and denominator.
    
    
  * Lemma <ref> is in fact an “if and only if” – that is, any Taylor series satisfying either of the conclusions of the Lemma must be a rational function (see e.g. <cit.>). However, we do not use the converse of the result.
    
    
  * The eigenvalues of W can be related to the poles of the corresponding rational function. Specifically,
    
    (λ I - W) = 0 λ^d + ∑_j=1^d q_j λ^d-j = 0 ( q(1λ) = 0  or λ=0=q_d ).

    
    
  * Furthermore, the eigenvectors of W can be written in terms of the poles: Suppose q(1λ) = 1 + ∑_j=1^d q_j (1λ)^j = 0. Let u^T = (λ^d-1, λ^d-2, ⋯, λ, 1) ∈ℝ^1 × d (or λ=0 and q_d=0). Then u^T W = λ· u^T or, equivalently, W^T u = λ· u.
    
    
  * If the eigenvalues of W are distinct (and real), then the above remark produces a complete (real) eigenbasis, and so we can diagonalize the matrix W (over ). A diagonal matrix W is particularly convenient for our streaming algorithm (<Ref>).    
    
    
  * If the eigenvalues of W are contained in the unit circle (i.e., |λ|≤ 1), then the computation of W^k is numerically stable. By Equation <ref>, this is equivalent to the poles of the rational function having magnitude ≥ 1.
    
    
  * Suppose the rational function r has distinct poles 1/θ_i, and can be presented as[Diagonalizing the matrix W given by <Ref> corresponds to finding a representation of this form.]
    
    r(x) = ∑_i ∈ []ω_i/1 - θ_i · x

    via a partial-fraction decomposition (instead of as a ratio of polynomials as in <Ref>).
    (Indeed, our rational function in <Ref> is presented in this form in <ref>.)
    Then, we have a simpler closed form for the sequence,
    
    r_k  = ∑_i ∈ []ω_i θ_i^k.

    This immediately yields a matrix representation: W is a diagonal matrix with W_i,i=θ_i for all i, and the vectors u and v simply need to satisfy u_i v_i = ω_i for all i.
    This matrix representation may be more convenient than the one given by <Ref>; we use this approach extensively in <ref>.

    
  * The form of the matrix W in <Ref> is known as a “companion matrix” <cit.>. Namely, the matrix W is a companion to the denominator polynomial q(x) since the roots of the polynomial correspond to the eigenvalues of the matrix (per <ref>).


Our choice in how to represent r(x) has some ramifications for the design of the noise generation (multiplication) algorithm in the next section:


    The representation in <Ref> assumes the numerator has lower degree than the denominator, i.e., (p)<d versus (q)≤ d. But, in this work, we typically consider rational functions with the numerator and denominator both having the same degree, i.e., (p)=(q)=d.
    This case can either be handled by padding (i.e., increment d ↦ d+1 and set q_d+1=0) or by including an additive constant:
    
    r(x) = p(x)/q(x) = t + p(x)/q(x) 
              = p(x) + tq(x)/q(x),
                where   
              t = p_d/q_d  and  p(x) = p(x) - t q(x),

    so that (p)≤ d-1.  We then have
    
    r_k = u^T W^k v + t 𝕀[k=0]

    with (u, W, v) as given in <ref> applied to p(x)/q(x). Here 𝕀 is the indicator function taking value 1 if the condition holds and 0 otherwise.        
    We prefer the representation with an additive constant, as it is more efficient to implement algorithmically.

We note the indicator function can also be written as 𝕀[k=0] = 0^k for k ≥ 0.
Thus we can incorporate the t 𝕀[k=0] into the matrix representation by appending a row/column of 0s to W, appending 1 to u and t to v:

    ( [ u; 1 ])^T ( [   W   0; 0^T   0 ])^k ( [ v; t ])
        = ( [ u; 1 ])^T ( [ W^k   0; 0^T 0^k ]) ( [ v; t ])
        
        = u^T W^k v + t ·𝕀[k=0].

This increases the dimension d by 1, which is equivalent to padding. However, for efficient streaming algorithm we introduce next, we can save one memory buffer by working directly with the r_k = u^T W^k v + t 𝕀[k=0] representation rather than the pure matrix representation in <ref>.




 §.§ Efficient Sampling via Multiplication
 

The generating function view outlined in <Ref> gives us a matrix factorization from any function r(x). Namely, A^() = M(b,) · M(c,) for b(x)=r(x)/(1-x) and c(x)=1/r(x). This turns out to be a good matrix factorization as long as r(x) ≈√(1-x).
The other desideratum is being able to efficiently sample noise according to the matrix factorization.

In this subsection we show that if r(x) is a rational function of low degree, then there exists an efficient sampling algorithm or, equivalently, an efficient algorithm for streaming multiplication by the lower-triangular Toeplitz matrix M(r(x),).
This algorithm relies on the representation in <ref>.
The final missing ingredient, which we provide in <Ref>, is to instantiate r(x).


Our algorithmic task is as follows. Fix a lower-triangular Toeplitz matrix M(r,) ∈ℝ^× with a generating function r. We are given as input a matrix Z ∈ℝ^× and must produce as output Z M(r,) · Z ∈ℝ^×. 
The obvious algorithm for this task is to use standard matrix multiplication; this would take O(^2) time and O() space. Another approach is to exploit the Toeplitz structure of M(r,) which makes the matrix multiplication a convolution; this means it can be accelerated using the fast Fourier transform to take O(log) time, but this is not a streaming algorithm and it still requires O() space. In order to reduce the space required we need an algorithm that is tailored to the streaming setting.

For streaming prefix sum applications (and private deep learning applications like DP-SGD and DP-FTRL in particular), the input Z consists of  independent samples from a Gaussian and Z = M(r,) · Z is the correlated noise we add to our private learning procedure – each row corresponds to one training iteration and each column to one parameter of the model.
In particular, the rows Z can be sampled as needed – since they are independent noise – and the output Z can be returned one row at a time. 

This is our streaming setting: At each iteration k, our algorithm receives as input the next row Z_k,·∈ℝ^ and must output the next row Z_k,·∈ℝ^.
Our goal is to develop an algorithm that runs in time and space O() per iteration (i.e., O() time in total).
Reducing the space usage relies on the fact that we do not need to store all of Z or Z in the streaming setting. Accomplishing this relies on the structure of the lower triangular Toeplitz matrix M(r,) ∈ℝ^d × d generated by a rational function 

    r(x) = t + ∑_i=0^d-1 p_i x^i/1+∑_j=1^d q_j x^j = r_0 + r_1 x + r_2 x^2 + r_3 x^3 + ⋯.

Following <ref> and <ref>, we use the representation of the sequence r_k in terms of matrix powers: r_k = u^T W^k v + t 𝕀[k=0].  This is the representation that we use for our algorithm. 






















For a given via a rational generating function r in the matrix representation of <ref>, <ref> in fact computes Ẑ = M(r) Z in row-by-row streaming fashion:

    <Ref> taskes as input a stream Z_0,·,Z_1,·, ⋯, Z_-1,·∈ℝ^1 × and outputs a stream Z_0,·,Z_1,·, ⋯, Z_-1,·∈ℝ^1 ×. 
    It takes parameters u,v ∈ℝ^d × 1 and W ∈ℝ^d × d and t ∈ℝ.
    At each iteration k, output satisfies 
    
    Z_k,· = ∑_j=0^k r_k-j· Z_j,·  = t Z_k,· + ∑_j=0^k u^T W^k-j v · Z_j,· ,
 
    where r_k = u W^k v + t 𝕀[k=0].
    The space usage is O(d  + d^2) and the runtime per iteration is dominated by a matrix multiplication W · S, where S ∈ℝ^d ×.


    By induction, for all k, we have
    
    S_k+1 = ∑_j=0^k W^k-j v Z_j,·

    and, hence,
    
    t Z_k,· + u^T S_k+1 = t Z_k,· + ∑_j=0^k u^T W^k-j v Z_j,· = ∑_j=0^k r_k-j Z_j,· = Z_k,·,

    where Z = M(r,) Z.
    
    The state of the algorithm after iteration k is given by S_k+1∈ℝ^d ×.
    The space usage of the algorithm is dominated by storing this state (d registers) plus storing the parameters u,v,W (2d+d^2 registers).
    More precisely, the only space that the algorithm requires is the state, plus whatever registers are required to perform the update S_k+1 = v Z_k,· + W · S_k. (The output computation Z_k,· = t Z_k,· + u^T S_k+1 requires  registers to store the output.)
    The exact number of registers required for this update depends on the structure of v and W. 
    If W is a diagonal matrix or has the structure given in <ref>, then the update can be performed “in place” and we only require  extra registers to store the input Z_k,·.
    In general, we can always perform this update using d additional registers.
    
    The running time of one iteration consists of the update S_k+1 = v Z_k,· + W · S_k and the matrix-vector product Z_k,· = u^T S_k+1.
    The matrix-matrix product will dominate over the matrix-vector products.
    The runtime of the matrix-matrix product depends on the structure of W. If W is sparse (e.g., if it is diagonal), then this is O(d) time.
    In general, we can perform this update in O(d^2) time.


We remark that, depending on the structure of W, we can reduce the space used by our algorithm. We use this optimization in practice to effectively reduce the parameter d by one.


    Suppose W has a zero-column (or, equivalently, a zero-row), i.e., ∃ j ∈ [n]  ∀ i ∈ [n]    W_i,j=0.
    Then <Ref> does not need to store the corresponding row of S.


Depending on the structure of W, we can often use one fewer buffer. If row i of W is all zeros (Remark about c_d = q_d = 0 as padding being the main example), Wi0 = 0⃗, then we need not store buffer Si: as during the “update” step we essentially zero buffer i. Instead, we add the contribution of the current Z_k that would have been carried in that buffer, u_i v_i Z_k, directly to the output Ẑk:.

Similarly, if column j of W is all zeros, W:j, then we never read buffer j in the buffer update step, and so we need not store buffer Sj:, again adding the current iteration contribution u_j v_j Z_k directly to Ẑk:.








§ FACTORIZATIONS VIA RATIONAL FUNCTION APPROXIMATION



In this section we prove <Ref>.
We follow the generating functions view presented in <Ref> and use the corresponding sampling algorithm (<ref>).
The proof is split into three steps. 
In <Ref> we recap how a generating function r yields a valid matrix factorization
 and prove that, if r(x) ≈√(1-x), then this yields an approximately optimal factorization.
In <Ref> we construct a low-degree rational function r that appropriately approximates the square root.
Finally we assemble the parts of the proof in <Ref>.






 §.§ Reduction to Approximating the Square Root


<Ref> shows that the optimal matrix factorization is given by ^()= where ==M(f,) and f(x)=1/√(1-x).
Suppose we have a generating function r(x) ≈√(1-x). 
Let b(x) r(x)/1-x≈1/√(1-x) = f(x) and c(x) 1/r(x)≈1/√(1-x). 
Since b(x) · c(x) = r(x)/1-x·1/r(x) = 1/1-x = g(x), we always obtain a valid factorization M(b,) · M(c,) = ^().
We can bound the objective value of this factorization in terms of the approximation r(x) ≈√(1-x).


    Let τ > 0 and ∈ℕ.
    Let r : C →ℂ, where C { z ∈ℂ : |z|<1 } is the open unit disc in the complex plane centered at zero.
    Let 
    γ_τsup{|r(x)-√(1-x)| : x ∈ℂ, |x|=exp(-τ)}.

    Let b(x) r(x)/1-x and c(x) 1/r(x) and f(x) = 1/√(1-x).
    Then M(b,) · M(c,) = M(c,) · M(b,) = A^() = M(f,)^2, where M(·,) is defined in Equation <ref> and A^() is defined in Equation <ref>.
    Moreover, if γ_τ≤(1-exp(-2τ)/4)^2, then
    
    M(b,)_2 →∞   ≤M(f,)_2 →∞ +  exp(τ) ·γ_τ/√(exp(2τ)-1),
    M(c,)_1 → 2   ≤M(f,)_1 → 2 + exp(τ) ·γ_τ/√((exp(2τ)-1)^2 - 2^7/2γ_τexp(4τ)).


We will end up setting τ=1/2. This ensures that exp(τ) ·γ_τ/√(exp(2τ)-1)≤√()·exp(1/2) ·γ_τ < 2 √(n)γ_τ and, if ≥ 6 and γ_τ≤1/32 ·^2, then exp(τ) ·γ_τ/√((exp(2τ)-1)^2 - 2^7/2γ_τexp(4τ))≤√(2)·√()·exp(1/2) ·γ_τ < 2 √(n)γ_τ.

Let b(x) = ∑_k=0^∞ b_k x^k, c(x) = ∑_k=0^∞ c_k x^k, and f(x) = ∑_k=0^∞ f_k x^k.
By linearity and the triangle inequality,

    M(b,)_2 →∞   = M(f,) + M(b-f,)_2 →∞≤M(f,)_2 →∞ + M(b-f,)_2 →∞ , 
    M(c,)_1 → 2   = M(f,) + M(c-f,)_1 → 2≤M(f,)_1 → 2 + M(c-f,)_1 → 2 .

Thus it suffices to bound M(b-f,)_2 →∞ = √(∑_k=0^-1 (b_k-f_k)^2) and M(c-f,)_1 → 2 = √(∑_k=0^-1 (c_k-f_k)^2).
We bound these using a weighted version of Parseval's identity.


    Let g,h : C →ℂ be analytic, where C { z ∈ℂ : |z|<1 } is the open unit disc in the complex plane centered at zero.
    Let g(x) = ∑_k=0^∞ g_k x^k and h(x) = ∑_k=0^∞ h_k x^k for x ∈ C, where g_k,h_k ∈ℂ.
    Let τ > 0. Then
    
    1/2π∫_-π^π |g(exp(√(-1)θ-τ)) - h(exp(√(-1)θ-τ))|^2 dθ   = ∑_k=0^∞ |g_k-h_k|^2 ·exp(-2 τ k) 
       ≥exp(-2 τ (-1)) ∑_k=0^-1 |g_k-h_k|^2 .



    Let w(x)=g(x)-h(x) = ∑_k=0^∞ w_k x^k, where w_k=g_k-h_k. Then 
    
    1/2π∫_-π^π |g(exp(√(-1)θ-τ)) - h(exp(√(-1)θ-τ))|^2 dθ
       = 1/2π∫_-π^π |w(exp(√(-1)θ-τ))|^2 dθ
       = 1/2π∫_-π^π w(exp(√(-1)θ-τ)) ·w(exp(√(-1)θ-τ))dθ
       = 1/2π∫_-π^π( ∑_k=0^∞ w_k · (exp(√(-1)θ-τ))^k ) ·( ∑_ℓ=0^∞ w_ℓ· (exp(√(-1)θ-τ))^ℓ)dθ
       = ∑_k,ℓ=0^∞1/2π∫_-π^π  w_k · (exp(√(-1)θ-τ))^k ·w_ℓ· (exp(-√(-1)θ-τ))^ℓdθ
       = ∑_k,ℓ=0^∞ w_k ·w_ℓ·exp(-(k+ℓ)τ) ·1/2π∫_-π^πexp((k-ℓ)√(-1)θ) dθ
       = ∑_k,ℓ=0^∞ w_k ·w_ℓ·exp(-(k+ℓ)τ) ·𝕀[k-ℓ=0] 
       = ∑_k=0^∞ |w_k|^2 ·exp(-2kτ) 
       = ∑_k=0^∞ |g_k-h_k|^2 ·exp(-2τ k).



Recall b(x) r(x)/1-x and f(x) = 1/√(1-x), whence b(x)-f(x) = r(x)-√(1-x)/1-x.
Now we have 

    M(b-f,)_2 →∞   = √(∑_k=0^-1 (b_k-f_k)^2)
       ≤exp(τ (-1)) √(1/2π∫_-π^π |b(exp(√(-1)θ-τ)) - f(exp(√(-1)θ-τ))|^2 dθ)
       = exp(τ (-1)) √(1/2π∫_-π^π|r(exp(√(-1)θ-τ)) - √(1-exp(√(-1)θ-τ)))|^2/|1-exp(√(-1)θ-τ))|^2dθ)
       ≤exp(τ (-1)) √(1/2π∫_-π^πγ_τ^2/|1-exp(√(-1)θ-τ))|^2dθ)
       = exp(τ (-1)) ·γ_τ·√(exp(τ)/4π∫_-π^π1/cosh(τ) + cos(θ)dθ)
       = exp(τ (-1)) ·γ_τ·√(exp(τ)/21/√(cosh(τ)^2-1))
       = exp(τ (-1)) ·γ_τ/√( 1-exp(-2τ) ) = exp(τ) ·γ_τ/√(exp(2τ) -1 )≤exp(τ) ·γ_τ/√( 2τ).

In the above we use the identity

    |1-exp(√(-1)θ-τ)|^2    = (1-e^-τcosθ)^2 + (e^-τsinθ)^2 = 1 - 2 e^-τcosθ + e^-2τ (cos^2 θ + sin^2 θ) 
        = e^-τ·( e^τ - 2 cosθ + e^-τ) = 2 e^-τ·( coshτ + cosθ)

and the integral 

    ∀σ>1      1/2π∫_-π^π1/σ+cos(θ)dθ = 1/2π∫_-π^π1/σ+sin(θ)dθ = 1/√(σ^2-1),

which follows from the derivative d/dθ[ 2/√(σ^2-1)tan^-1(1 + σtan (θ/2)/√(σ^2-1)) ] = 1/σ+sin(θ).

Similarly, recall c(x) 1/r(x) and f(x) = 1/√(1-x), whence c(x) - f(x) = √(1-x)-r(x)/r(x) ·√(1-x).
Assuming γ_τ < (1-exp(-2τ))^2/2^7/2, we have


    M(c-f,)_1 → 2
       = √(∑_k=0^-1 (c_k-f_k)^2)
       ≤exp(τ (-1)) √(1/2π∫_-π^π |c(exp(√(-1)θ-τ)) - f(exp(√(-1)θ-τ))|^2 dθ)
       = exp(τ (-1)) √(1/2π∫_-π^π|r(exp(√(-1)θ-τ)) - √(1-exp(√(-1)θ-τ))|^2/| r(exp(√(-1)θ-τ))  ·√(1-exp(√(-1)θ-τ))|^2dθ)
       = exp(τ (-1)) √(1/2π∫_-π^π|r(exp(√(-1)θ-τ)) - √(1-exp(√(-1)θ-τ))|^2/| 1-exp(√(-1)θ-τ) +(r(exp(√(-1)θ-τ)) -√(1-exp(√(-1)θ-τ))) ·√(1-exp(√(-1)θ-τ))|^2dθ)
       ≤exp(τ (-1)) √(1/2π∫_-π^πγ_τ^2/(| 1-exp(√(-1)θ-τ)|-γ_τ·√(|1-exp(√(-1)θ-τ)|))^2dθ)
       ≤exp(τ (-1)) √(1/2π∫_-π^πγ_τ^2/| 1-exp(√(-1)θ-τ)|^2 -2√(2)γ_τdθ)
       = exp(τ (-1)) ·γ_τ·√(exp(τ)/4π∫_-π^π1/cosh(τ) - √(2)exp(τ) γ_τ + cos(θ) dθ)
       = exp(τ (-1)) ·γ_τ·√(exp(τ)/21/√((cosh(τ) - √(2)exp(τ) γ_τ )^2-1))
       = exp(τ) ·γ_τ/√(2exp(τ)√((cosh(τ)-√(2)exp(τ)γ_τ)^2-1))
       = exp(τ) ·γ_τ/√((exp(2τ)(1-2√(2)γ_τ) + 1)^2-4exp(2τ))
       = exp(τ) ·γ_τ/√(1 -2(1+2^3/2γ_τ)exp(2τ) + (1-2^3/2γ_τ)^2exp(4τ))
       = exp(τ) ·γ_τ/√((exp(2τ)-1)^2 - 2^5/2γ_τexp(2τ) - 2^5/2γ_τexp(4τ) + 2^3 γ_τ^2 exp(4τ))
       ≤exp(τ) ·γ_τ/√((exp(2τ)-1)^2 - 2^5/2γ_τexp(2τ) ( 1 + exp(2τ) ) )
       ≤exp(τ) ·γ_τ/√((exp(2τ)-1)^2 - 2^7/2γ_τexp(4τ) ) .


Combining the bounds yields the result.





 §.§ Rational Approximation of √(x)
 

<Ref> reduce the problem to finding a low-degree rational function that uniformly approximates √(1-x) for x ∈ℂ with |x|≤ 1.
Fortunately, such functions are known to exist.
In this subsection we construct the required rational functions. The proof is included for completeness and because it provides an explicit construction with explicit bounds. 


For notational simplicity and consistency with the literature, we rescale so that we are approximating √(x) instead of √(1-x). To get an equivalent result, we must uniformly approximate √(x) for x ∈ℂ with |x-1/2|≤1/2.

<cit.> gave an explicit rational function that approximates √(x) for x ∈ [0,1]. Specifically, let 

    p(x) = ∏_k=0^d-1 (x + exp(-k/√(d)))     and     r(x) = √(x)· (p(√(x))-p(-√(x)))/p(√(x))+p(-√(x)).

Then r(x) is a rational function of degree ⌈ d/2 ⌉ – specifically, r(x) is the ratio where the numerator is a polynomial in x of degree ⌈ d/2 ⌉ and the denominator is a polynomial in x of degree ⌊ d/2 ⌋. Furthermore, 

    sup_x ∈ [0,1] |r(x)-√(x)| ≤ 3 ·exp(-√(d)).
 
The surprising aspect of this result is that the dependence on the degree is exponential rather than polynomial, as is the case if we restrict to polynomial approximations, rather than rational approximations.
Newman also showed that this result was optimal up to constant factors in the degree d and the approximation error.
Subsequent work <cit.> obtained optimal constants for this approximation.

In order to apply <Ref>, we need to extend Newman's result into the complex plane. We provide a proof adapted from that of <cit.>.


    Let d ∈ℕ with d ≥ 2. Then there exists a rational function r of degree d with real coefficients and real negative simple poles such that
    
    |r(x) - √(x)| ≤(4 + 4/π) ·exp(-π√(1/2⌊d-1/2⌋)) ≤ 6 ·exp(-π/2√(d-2))

    for all x ∈ℂ with (x) ≥ 0 and |x| ≤ 1.
    Specifically,
    
    r(x) = 2h/π∑_k=-d_-^d_+x ·exp(hk)/x+exp(2hk),    where    d_+ = ⌊d-1/2⌋,  d_- = ⌈d-1/2⌉,   and    h = π/√(2 d_+).



    We invoke <Ref>. Let d_- = ⌈(d-1)/2⌉ and d_+ = ⌊(d-1)/2⌋ and let h ∈ (0,4.9) be determined later. Let r(x) = r_d_+,d_-,h(x) be as in Equation <ref>. Clearly r(x) is a rational function of degree d_++d_-+1 = d with real coefficients and the poles are given by -exp(2hk) ∈ (-∞,0) for k∈{-d_-,-d_-+1,⋯,d_+}. It remains to simplify the approximation guarantee and set h.
    From Equation <ref>, for all x ∈ℂ with (x) ≥ 0 and |x| ≤ 1, we have
    
    |r(x)-√(x)|    ≤ 2 √(|x|)( 1/exp((1-c)π^2 / h) - 1 + 1/exp((1+c)π^2 / h)-1) 
                  + 2h/π (exp(h)-1)( |x| ·exp(-hd_+) + exp(-hd_-)) 
       ≤ 2 ( 1/exp((1-|c|)π^2 / h) - 1 + 1/exp((1+|c|)π^2 / h)-1) |x| ≤ 1& symmetry in c
                  + 4h exp(-hd_+)/π (exp(h)-1)d_- ≥ d_+
       ≤ 2 ( 1/exp(π^2 / 2h) - 1 + 1/exp(π^2 / h)-1) + 4/π·exp(-hd_+) ·h/exp(h)-1|c| = |1/π(x)| ≤1/2
       ≤ 4 ·exp(-π^2/2h) + 4/π·exp(-hd_+) 
       = (4 + 4/π) ·exp(-π√(d_+/2))  
       = (4 + 4/π) ·exp(-π√(1/2⌊d-1/2⌋)) 
       ≤ 6 ·exp(-π/2√(d-1)).

    The fourth inequality <ref> follows from two simple facts:  
    First, h ≥ 0 implies h/exp(h)-1≤ 1.
    Second, 1/a-1+1/a^2-1≤2/a when a = exp(π^2/2h) ≥ 1+√(3) or, equivalently, h ≤π^2/2 log(1+√(3))≈ 4.91.
    The penultimate equality <ref> follows by setting h = π/√(2d_+)≤π/√(2) < 2.5.



    Let d ≥ 2 be an integer. Then there exists a rational function r of degree d with real coefficients such that all poles are simple and > 1 and we have
    
    |r(x) - √(1-x)| ≤ 8  ·exp(-π/2√(d-2))

    for all x ∈ℂ with |x| ≤ 1.
    Specifically,
    
    r(x) = 2h√(2)/π∑_k=-d_-^d_+(1-x) ·exp(hk)/1-x+2·exp(2hk) = 2h√(2)/π∑_k=-d_-^d_+exp(hk) - 2 ·exp(3hk)/1+2·exp(2hk)-x

    where d_+ = ⌊(d-1)/2⌋, d_- = ⌈(d-1)/2⌉, and h = π/√(2d_+).


    Let r be the rational function of degree ≤ d promised by <Ref>.
    Define r(x) = √(2)· r(1-x2).
    Suppose x ∈ℂ with |x|≤ 1.
    Then (1-x2) = 1/2 -1/2(x) ≥ 0, |1-x2| ≤1+|x|2≤ 1, and 
    |r(x)-√(1-x)|  = √(2)·| r(1-x2)-√(1-x2)| ≤√(2)·(4 + 4/π) ·exp(-π√(1/2⌈d-1/2⌉)) ≤ 8 exp(-π/2√(d-1)).

    If x is a pole of r, then 1-x2 is a pole of r; since the poles of r are negative we have 1-x2 < 0 and, hence, x>1.


The general form of the approximation promised by <Ref> is given below.


    Let d_+,d_- ∈ℕ and h>0. For x ∈ℂ∖ [-exp(2d_+h),-exp(-2d_-h)], define 
    
    r_d_+,d_-,h(x) 2h/π∑_k=-d_-^d_+x ·exp(hk)/x+exp(2hk).

    Then, for all x ∈ℂ with (x) ≥ 0, we have
    
    |r_d_+,d_-,h(x)-√(x)|    ≤ 2 √(|x|)( 1/exp((1-c)π^2 / h) - 1 + 1/exp((1+c)π^2 / h)-1) 
             +  2h/π (exp(h)-1)( |x| ·exp(-hd_+) + exp(-hd_-)),

    where c = 1/π(x) ∈[-1/2,1/2].


    Fix x ∈ℂ with (x) > 0.[For the proof we assume (x)>0. The case (x)=0 follows by continuity.]
    Let (x) = tan^-1((x)/(x)) = cπ∈(-π/2,π/2).
    
    We begin with the Cauchy distribution integral 
    
    ∫_0^∞1/1+u^2du = ∫_0^∞( d/dutan^-1(u) ) du = π/2.

    We perform the variable substitutions u = v/√(x) and v=exp(s) and rearrange to obtain
    
    √(x) = 2/π∫_0^∞x/x+v^2dv = 2/π∫_-∞^∞x ·exp(s)/x+exp(2s)ds .

    Now we make the approximation
    
    √(x) = 2/π∫_-∞^∞x ·exp(s)/x+exp(2s)ds ≈2h/π∑_k ∈ℤx ·exp(hk)/x+exp(2hk)≈2h/π∑_k = -d_-^d_+x ·exp(hk)/x+exp(2hk) = r_d_+,d_-,h(x).

    
    Thus all that remains to complete the proof is to make the the two approximations precise. 
    Intuitively, the first approximation replaces the continuous integral with a discrete-but-infinite Riemann sum. If the function is smooth and h is small, this approximation should be good.
    The second approximation truncates the infinite sum, which is a good approximation as long as the tails of the function decay rapidly and hd_+ and hd_- are large.
    
    Define f(z) 2h/πx ·exp(hz)/x+exp(2hz).
    Let
    
    f(t)    ∫_-∞^∞ f(z) ·exp(-2π z t √(-1)) dz 
       = ∫_-∞^∞2h/πx ·exp(hz - 2π z t √(-1))/x + exp(2hz)dz 
       = 2h/π∫_-∞^∞exp(- 2π z t √(-1))/exp(-hz) + exp(hz-log x)dz 
       = 2h/π∫_-∞^∞exp(- 2π (y+log x/2h) t √(-1))/exp(-hy-1/2log x) + exp(hy-1/2log x)dy z = y+log x/2h
       = √(x)·2h/π·exp(-π t log x/h√(-1)) ·∫_-∞^∞exp(- 2π y t √(-1))/exp(-hy) + exp(hy)dy 
       =  √(x)·2h/π·exp(-π t log x/h√(-1)) ·1/2∫_-∞^∞exp(- 2π y t √(-1))/cosh(hy)dy 
       =  √(x)·2h/π·exp(-π t log x/h√(-1)) ·π/2h ·cosh(π^2 t / h)
       =  √(x)·exp(-π t log x/h√(-1)) ·1/cosh(π^2 t / h) .

    The penultimate equality (<ref>) is the Fourier transform of the hyperbolic secant <cit.>.
    The Poisson summation formula states that 
    ∑_k ∈ℤ f(k) = ∑_t ∈ℤf(t).

    Per Equation <ref>, we have 
    f(0) = ∫_-∞^∞ f(z) dz = ∫_-∞^∞2h/πx ·exp(hz)/x+exp(2hz)dz = ∫_-∞^∞2/πx ·exp(s)/x+exp(2s)ds = √(x) .

    Thus the Poisson summation formula allows us to analyze the first approximation in Equation <ref>.
    Namely,
    
    | √(x) - 2h/π∑_k ∈ℤx ·exp(hk)/x+exp(2hk)|
       = | f(0) - ∑_k ∈ℤ f(k) | 
       =| f(0) - ∑_t ∈ℤf(t) | 
       ≤∑_t=1^∞ |f(t)| +  |f(-t)| 
       = ∑_t=1^∞ |√(x)| ·( |exp(-π t log x/h√(-1))| + |exp(π t log x/h√(-1))| ) ·1/cosh(π^2 t / h)
       = ∑_t=1^∞ |√(x)| ·( exp((π t log x/h)) + exp((-π t log x/h)) ) ·1/cosh(π^2 t / h)|exp(v√(-1))|=exp((v√(-1)))=exp((-v))
       = √(|x|)∑_t=1^∞ 2cosh(πt/h(log x) ) ·1/cosh(π^2 t / h)
       = √(|x|)∑_t=1^∞ 2cosh(πt/h· cπ) ·1/cosh(π^2 t / h)(log x) = (x) = cπ
       = 2√(|x|)∑_t=1^∞exp(c π^2 t / h) + exp(-c π^2 t / h)/exp(π^2 t / h) + exp(-π^2 t / h)
       ≤  2√(|x|)∑_t=1^∞exp(cπ^2 t / h) + exp(-cπ^2 t / h)/exp(π^2 t / h) + 0
       = 2√(|x|)∑_t=1^∞exp(-(1-c)π^2 t / h) + exp(-(1+c)π^2 t / h) 
       = 2√(|x|)( 1/exp((1-c)π^2 / h) - 1 + 1/exp((1+c)π^2 / h)-1) .

    
    The final equality (<ref>) follows from the usual geometric series formula: If a>0, then ∑_t=1^∞exp(-at) = exp(-a)/1-exp(-a) = 1/exp(a)-1.
    
    Next we analyze the second approximation in Equation <ref>, which requires bounding the tails of f.
    On one side, for k ∈ℕ, we have 
    
    |f(k)|    = |2h/πx ·exp(hk)/x+exp(2hk)| 
       = 2h/π|x|/|x ·exp(-hk)+exp(hk)|
       = 2h/π|x|/√((x ·exp(-hk)+exp(hk))^2+(x ·exp(-hk)+exp(hk))^2)
       ≤2h/π|x|/|(x ·exp(-hk)+exp(hk)|
       = 2h/π|x|/|(x) ·exp(-hk)+exp(hk)|hk ∈ℝ
       = 2h/π|x|/(x) ·exp(-hk)+exp(hk)(x)>0
       ≤2h/π|x|/exp(hk).

    On the other side, for k ∈ℕ, we have
    
    |f(-k)|    = |2h/πx ·exp(-hk)/x+exp(-2hk)| 
       = 2h/πexp(-hk)/|1+exp(-2hk)/x|
       ≤2h/πexp(-hk)/|(1+exp(-2hk)/x)|
       = 2h/πexp(-hk)/|1+exp(-2hk)/(x)|(1/x)=1/(x)
       = 2h/πexp(-hk)/1+exp(-2hk)/(x)(x)>0
       ≤2h/πexp(-hk).

    Combining Equations <ref> and <ref>, we have
    
    |2h/π∑_k ∈ℤx ·exp(hk)/x+exp(2hk) - r_d_+,d_-,h(x)|
               =2h/π|∑_k ∈ℤx ·exp(hk)/x+exp(2hk) - ∑_k = -d_-^d_+x ·exp(hk)/x+exp(2hk)| 
       = |∑_k=d_++1^∞ f(k) + ∑_k=d_-+1^∞ f(-k) | 
       ≤∑_k=d_++1^∞ |f(k)| + ∑_k=d_-+1^∞ |f(-k)| 
       ≤2h/π∑_k=d_++1^∞|x|/exp(hk) + 2h/π∑_k=d_-+1^∞exp(-hk) 
       = 2h/π |x| exp(-h(d_++1))/1-exp(-h) + 2h/πexp(-h(d_-+1))/1-exp(-h)
       = 2h/π (exp(h)-1)( |x| ·exp(-hd_+) + exp(-hd_-)).

    Combining Equations <ref> and <ref> yields the result in Equation <ref>.




 §.§ Putting Everything Together
 

Now we can assemble the tools we have developed to prove our main result (<Ref>).
In <Ref> we connected near-optimal lower triangular Toeplitz matrix factorizations to approximating the square root function.
In <Ref> we connected low-degree rational functions to efficient streaming algorithms.
In <Ref> we gave a low-degree rational function that approximates the square root function.
The combination of these three tools gives a near-optimal matrix factorization with a corresponding efficient streaming algorithm.


    Let , d, ∈ℕ satisfy n ≥ 5 and d ≥ 2+(12+4log n/π)^2.
    There exists a lower triangular Toeplitz matrix factorization B, C, A=BC ∈ℝ^× with the following properties.
    
        
  * Near-Optimality: 
        Let B^*=C^* = M(f,) for f(x)=1/√(1-x) be the optimal lower triangular Toeplitz matrix factorization of size . Then
        
    B_2 →∞   ≤B^*_2 →∞ + 16 √(n)exp(-π/2√(d-2)),
    C_1 → 2   ≤C^*_1 → 2 + 16 √(n)exp(-π/2√(d-2)).

        
  * Efficiency:
        There is an algorithm that takes as input a stream Z_1,·,Z_2,·,⋯,Z_,·∈ℝ^1 × and outputs a stream Z_1,·,Z_2,·,⋯,Z_,·∈ℝ^1 × such that Z = C^-1 Z and the space and time (per iteration) required by the algorithm is O(d).
    

To obtain <Ref> in the introduction from <Ref>, we simply need to set d=O(log(n/μ))^2 such that d ≥ 2 + (12+4log n/π)^2 and 16 √(n)exp(-π/2√(d-2)) ≤μ/4+log(n).

    <Ref> gives us a rational function r of degree d with real coefficients such that 
    
    |r(x) - √(1-x)| ≤ 8  ·exp(-π/2√(d-2))

    for all x ∈ℂ with |x| ≤ 1.
    
    Let b(x) r(x)/1-x and c(x) 1/r(x) and f(x) 1/√(1-x) and g(x)  1/(1-x) = b(x)c(x).
    Let B  M(b,), C  M(c,) ∈ℝ^× be as in <Ref>.
    By <Ref>, this is a valid matrix factorization, as B · C = M(b,) · M(c,) = M(g,).
    
    Let τ > 0 and 
    γ_τsup{|r(x)-√(1-x)| : x ∈ℂ, |x|=exp(-τ)}≤ 8  ·exp(-π/2√(d-2)).

    By <Ref>, if γ_τ≤(1-exp(-2τ)/4)^2, then
    
    M(b,)_2 →∞   ≤M(f,)_2 →∞ +  exp(τ) ·γ_τ/√(exp(2τ)-1),
    M(c,)_1 → 2   ≤M(f,)_1 → 2 + exp(τ) ·γ_τ/√((exp(2τ)-1)^2 - 2^7/2γ_τexp(4τ)).

    Set τ=1/2. 
    Then exp(2τ)-1 ≥ 1/n and exp(τ n) = exp(1/2). 
    Thus 
    exp(τ) ·γ_τ/√(exp(2τ)-1)≤√()·exp(1/2) ·γ_τ < 2 √(n)γ_τ
 and, if ≥ 5 and γ_τ≤(1/6 ·)^2, then 
    exp(τ) ·γ_τ/√((exp(2τ)-1)^2 - 2^7/2γ_τexp(4τ))≤exp(1/2) ·γ_τ/√(1/n^2 - 2^7/2γ_τexp(2/n))≤exp(1/2) γ_τ/√(1/n^2 - 2^7/2exp(2/5) / 36n^2)  < 2 √(n)γ_τ.

    Assuming n ≥ 5 and d ≥ 2 + (2/πlog(8 ·(6n)^2 ))^2, we have γ_τ≤(1/6n)^2 ≤(1-exp(-2τ)/4)^2, which means all the above inequalities hold. Combining inequalities yields the first part of the result.
    
    Per <Ref>, we can write 
    r(x) = t + p(x)/q(x)
 where (p)<d and (q)≤ d and t ∈ℝ.
    <Ref> states that there exist u,v ∈ℝ^d and W ∈ℝ^d × d such that we can express 
    r(x) = ∑_k=0^∞ r_k x^k      for     r_k = u^T W^k v + t 𝕀[k=0].
 Furthermore W is sparse in the sense that each row has only two nonzero entries.
    Now we can feed u,v,W,t as parameters to <Ref>.
    <Ref> tells us that <Ref> receives the rows of Z ∈ℝ^× as a stream and returns the rows of Z∈ℝ^× as a stream, where Z = M(r,) Z = C^-1 Z. (The fact that M(r,) = C^-1 follows from <Ref>.)
    The space required by the algorithm is O(d) and the time per iteration is dominated by computing the matrix multiplication WS for S ∈ℝ^d ×. Since W is sparse, this takes O(d) time.





§ FACTORIZATIONS VIA DIRECT OPTIMIZATION
 

In this section, we begin (in <ref>) by showing the following result:

Let r be a rational generating function. Then, for the matrix factorization of  given by = (r(x)/(1-x),) and = (1/r(x),), we can compute in “closed-form” (more specifically, in time ) the sensitivity and error of this matrix factorization, namely ^2_2 →∞ and ^2_1 → 2.


The result follows from <ref>.


We remark that the log dependence in <ref> comes (only) from the need to compute quantities like θ^k for θ∈ and k ≤, which can be done via iterated squaring. For practical  and modern hardware, this term can essentially be ignored, making the running time practically independent of . If the roots of the numerator and denominator of r are already known (as assumed in <ref>) then  can be computed in time (^2 log(n)); and such roots can always be found in time (()) (omitting the dependence on the desired numerical accuracy).

In <ref> we then show how <ref> can be used to directly optimize the parameters of a factorization to minimize  for a specific .



 §.§ Fast and Exact Error Computation

In addition to enabling efficient streaming multiplication, the rational function structure of the matrix factorizations we consider allows us to compute the total error of these mechanisms efficiently for any . 


Let r(x) = p(x)/q(x) be a rational function where (p) = (q) = with q(0) = p(0) = 1. Further, suppose p and q have pairwise distinct real roots, in particular there exist θ_i, _i ∈ (0, 1] for i ∈ [] such that

    p(x)    = (1 - _0 x) (1 - _1 x) ⋯ (1 - _d-1 x), 
    
        q(x)    = (1 - θ_0 x) (1 - θ_1 x) ⋯ (1 - θ_d-1 x)

with θ_i θ_j and θ̂_i θ̂_j when i  j.

Then, there exist _i, _i ∈ for i ∈ [] such that 

    r(x) = 1 + x (∑_i=0^-1ω_i/1 - θ_i x)
        and    (x) 1/r(x) = 1 + x (∑_i=0^-1ω̂_i/1 - θ̂_i x).

These _i, _i can in fact be computed in closed form, <ref>. Further, the coefficients of the corresponding sequences can be computed directly as 

    _i = 
    
    1      i = 0 
    ∑_j ∈ []_j θ_j^i-1    i > 0.
      and  _i = 
    
    1      i = 0 
    ∑_j ∈ []_j _j^i-1    i > 0.





















Before giving a proof, we observe that <ref> yields a diagonalized companion matrix of the form of <ref>, and using the representation of <ref> we can take
t = 1 - ∑_i=0^-1_i/θ_i so

    _k = t𝕀[k=0] + u W^k v [ 1; 1; ⋮; 1 ][  θ_0    0    0    0;    0  θ_1    0    ⋮;    ⋮    0    ⋱    0;    0    0    … θ_-1;      ]^k
        [   _0/θ_0;   _1/θ_1;        ⋮; _-1/θ_-1 ],

and hence we can efficiently multiply by M(r) using <ref>; a similar construction with (, ) enables efficient multiplication by M(1/r) using only  buffers. 

This construction is particularly useful for optimization, but assumes we already have factored p and q; one can of course always find the roots of p and q if they are given in some other polynomial representation.[While we only need to consider real roots in our constructions and our optimization, the results of this section go through for complex roots as well.] Specifically, if one is instead given the matrix-power representation of r, <ref> of <ref> shows these can directly be converted to a parameterization of 1/r.





<ref> follows from <ref> and the Taylor series

    x/1 - θ x = ∑_i=1^∞ωθ^i-1  x^i.

It remains to show we can find _i to instantiate <ref> for r(x) (the case of _i for (x) is completely symmetric). 

Letting

    q^(-i)(x) q(x)/1 - θ_i x = ∏_j ∈ []
     j ≠ i (1 - θ_j x),

we need _i that satisfy

    1 + x (∑_i=0^-1ω_i/1 - θ_i x)
    = q(x) + ∑_i=0^-1_i x q^(-i)(x)/q(x)
    = p(x)/q(x).

Since p(0) = q(0) = 1, p(x) - q(x) is divisible by x, so f(x) 1/x(p(x) - q(x)) is a polynomial of degree - 1.  Hence, we wish to solve

    ∑_i=0^-1_i q^(-i)(x) = f(x)
 
for _0, …, _ - 1. Equating the polynomial coefficients, this is a linear system of  equations and  unknowns. It remains to show that this system has full-rank or, equivalently, that the set of polynomials {q^(-i)}_i ∈ [] forms a basis for the vector space of polynomials of degree ≤ - 1.


Let z ∏_i ∈ [] -θ_i and w_i (∏_j  i (θ_i - θ_j)) and define

    ℓ(x) q(x)/z
     = (x - θ_0) (x - θ_1) ⋯ (x - θ_-1)
       and  ℓ^(-i) (x) ℓ(x) w_i/x - θ_i.

Then {ℓ^(-i)}_i ∈ [] is exactly the Lagrange polynomial basis for interpolating points at θ_0, …, θ_-1;
further, we have

    q^(-i)(x) 
       =  q(x)/-θ_i (x - θ_i) 
      = -z/θ_i w_iℓ^(-i)(x),

and since the θ_i are distinct and non-zero, -θ_i z/w_i is non-zero. It follows from Lagrange interpolation that

    f(x) 
      = ∑_i ∈ [] f(θ_i) ℓ^(-i)(x)
      = ∑_i ∈ [] f(θ_i)-θ_i w_i/ z q^(-i)(x)

and so

    ω_i = f(θ_i)-θ_i w_i/z.



It will be useful to define the prefix sums of the geometric series by

    γ_n(θ)  1 + θ + θ^2 + … + θ^n-1 = ∑_i=0^n-1θ^i.

with γ_∞(θ) = 1 / (1 - θ) and γ_0(θ) = 0.
When θ < 1, we have of course

    γ_n(θ) = ∑_i=0^n-1θ^i = 1 -θ^n/1 - θ.

(This use of γ is unrelated to the γ_2(A) factorization norm discussed earlier.)


Let ^() = (s(x), ) for a rational function 

    (x) = 1 + x (∑_i=0^-1_i/1 - _i x),

for distinct _i as per the assumptions of <ref>. Then, there exists a closed form expression for the sensitivity of ^(), (^()) = ^()_1→2 allowing computation in time and space 
(^2 log).


The squared sensitivity can be computed as

    ^2() = ∑_i=0^ - 1_i^2 
         = 1 + ∑_i=0^n-2(∑_j ∈ []_j ^i_j)^2
               by <ref>
       = 1 + ∑_i=0^n-2(∑_j ∈ []∑_k ∈ []_j_k ^i_j^i_k)
       = 1 +  ∑_j ∈ []∑_k ∈ [](∑_i=0^n-2_j_k (_j_k)^i)
       = 1 +  ∑_j ∈ []∑_k ∈ []_j_k γ_n-1(_j_k),

and taking a square root completes the proof.











Let ^() = ^() (^()) = ^()(r(x), ) = (r(x)/(1-x), ) for a rational function

    r(x) = 1 + x (∑_i=0^-1ω_i/1 - θ_i x)

for distinct θ_i as per the assumptions of <ref>.
Then, there exists a closed form expression for the variance induced in the nth prefix sum,
^(n)^2_2 →∞ 
allowing computation in time and space 
(^2 log n).


In order to make our notation more compact, we will compute sensitivity for an arbitrary row index  of (r(x)/(1-x), +1).  Since multiplication by (1-x) corresponds to taking prefix sums, we know = (t) where t_0 = r_0 = 1 and for any n > 0,

    t_n ∑_i=0^n r_i 
        = 1 + ∑_j ∈ []∑_i=0^n-1_j ^i_j
        = 1 + ∑_j ∈ []_j γ_n(_j)

and we calculate

    t_n^2 
        =  ( 1 + ∑_j ∈ []_jγ_n(_j))^2 
        = 1 + 2∑_j ∈ []_j γ_n(_j) 
            + ∑_j ∈ []∑_k ∈ []_j γ_n(_j)_k γ_n(_k).

It is thus sufficient to compute

    n:_2^2 
           =∑_i=0^n-1 t_i^2 
       = 1 + ∑_i=1^n-1 t^2_i
       = 1 + ∑_i=1^n-1( 1 + 2∑_j ∈ []_j γ_i(_j) 
            + ∑_j ∈ []∑_k ∈ []_j γ_i(_j) _k γ_i(_k)  ).

where we used <ref> for the last step.

Bringing the sum over i inside the sums over the buffers j and k, it will be sufficient to consider the following terms. For any j ∈ [b],

    Γ_j 
            ∑_i=1^n-1_jγ_i(_j) 
         = ∑_i=1^n-1_j(1 -_j^i)/1-_j  assuming 0 ≤ < 1
       = _j/1 - _j((n - 1) - ∑_i=1^n-1_j^i) 
       = _j/1 - _j(n - γ_n(1, _j)))

Similarly, for any j ∈ [], k ∈ [],  let

    Γ_j,k∑_i=1^n-1   _j γ_i( _j) _k γ_i(_k) 
       = ∑_i=1^n-1_j(1 -_j^i)/1-_j_k(1 -_k^i)/1-_k      since 0 ≤ < 1
       = _j_k/(1 - _j)(1 - _k)∑_i=1^n-1(1 - _j^i - _k^i + (_j _k)^i) 
       = _j_k/(1 - _j)(1 - _k)∑_i=0^n-1(1 - _j^i - _k^i + (_j _k)^i) 
       = _j_k/(1 - _j)(1 - _k)(n - γ_n(1, _j) - γ_n(1, _k) + γ_n(1, _j _k)).

Putting everything together we have

    n:_2^2 = ∑_i=0^ - 1 t_i^2 = n + 2∑_j ∈ []Γ_j + ∑_j ∈ []∑_k ∈ []Γ_j,k.

Since <ref> computes the L_2-norm squared for the row indexed by n, and this quantity is non-decreasing in n, it follows that

    ^(n)_2 →∞ = √(n-1:_2^2),

which can be computed using <ref>.



Rearranging <ref> and in particular expanding the terms Γ_j (<ref>) and Γ_j, k (<ref>) reveals there exists constants α_0, α_1 ∈ and ω̃_i ∈, θ̃_i ∈ (0, 1] for i ∈ [d + d^2] such that

    ^(n)_2 →∞ = α_0 + α_1 · + ∑_i ∈ [d + d^2]ω̃θ̃^ -1

where in particular

    α_1
     = 1 + 2∑_j ∈ []_j/1 - _j 
          + ∑_j ∈ []∑_k ∈ []_j_k/(1 - _j)(1 - _k).

We must have α_1 ≥ 0 as for a valid factorization the norm must remain positive for all , and lim_n →∞ω̃_i θ̃^ -1 = 0 (or another constant ω̃_i if θ̃_i = 1). For all the we consider (and we believe all useful ), we find α_1 > 0. Since sensitivity, <ref>, is non-decreasing in , for any such , eventually  must grow linearly. Hence, one must either increase  with  or optimize the for a specific  in order to hope for good performance. In particular, by optimizing for a specific  (<ref>), we must keep α small enough that the α_1 term remains relatively small (say, of the same order as α_0). However, for larger  this term can quickly explode; we see exactly this behavior in <ref> (Left column). 





 §.§ Direct Optimization of 

Inspection of <ref> reveals that the function (θ, ) → is in fact differentiable. With this result, together with the automatic differentiation capabilities of JAX <cit.>, we can in fact directly optimize for rational functions of a given degree  which minimize  for a given number of steps .

Some care is needed in the implementation. In particular, we find:

    
  * Using  precision is necessary, particularly when optimizing for > 10^6.
    
  * As  becomes large, max_i ∈ []θ_i tends to 1.  <ref> can be numerically unstable in such situations; this can be remedied by switching to a Taylor-series approximation taken around θ=1 for θ sufficiently near (or equal to) 1.
    
  * The optimization is more stable (particularly for larger ) if we add a barrier function
    
    ℓ(θ, ω) = 10^-7∑_i ∈ [] -log(θ_i) -log(_i)

    to ensure these parameters remain strictly positive. 
    
  * For large  (e.g >10^7), larger numbers of buffers  can make the optimization less stable. For example, with our current implementation we find for = 10^8, optimizing for more than =6 buffers actually slightly decreases performance (while theoretically it should only help).
    
    
    
    
    
    
    
    
    
    
    
    

With the above setup, we find the L-BFGS optimizer works well. Convergence generally takes less than a second running on CPUs on a modern workstration[There is approximately 5-15 seconds of overhead for JAX to just-in-time compile the loss function; this is only incurred once even if mechanisms are optimized for many different .], even from a naive initialization (initialization from <ref> would likely require even fewer iterations). Numerical results for these mechanisms are given in <ref>.

Proving the convergence properties of this approach and deriving more robust optimization algorithms are interesting directions for future work.



 §.§ Derivation of the Inverse Parameterization

We now show that the matrix-power representation of a rational function r(x) (see <Ref>) can be converted to a matrix-power representation for its reciprocal 1/r(x). This is useful for our algorithm.
Note that the dimension of the representation increases by one.

    Let u,v∈ℝ^d and W ∈ℝ^d × d. Assume ⟨ u , v ⟩ = 1.[This assumption is only made for simplicity and holds without loss of generality, as we can always rescale the function by r(0) = u^Tv: That is, let r(x)  r(x)/r(0), apply the result, and rescale back 1/r(x) = r(0) · 1/r(x).]
    Define r(x) = ∑_k=0^∞ u^T W^k v · x^k.
    Then
    
    1/r(x)   = 1 - ∑_ℓ=1^∞ u^T W (W - v u^T W)^ℓ-1 v · x^ℓ
       = ∑_k=0^∞u^T W_0^k v· x^k

    and
    
    1/r(x) · (1-x) = ∑_k=0^∞u^T W_1^k v· x^k

    where, for both β∈{0,1},
    
    W_β   = ( [           β       0 ⋯ 0;           v W - v u^T W ]) ∈ℝ^(d+1) × (d+1), 
    u   = ( [       1; - W^T u ]) ∈ℝ^d+1, 
    v   = ( [ 1; 0; ⋮; 0 ]) ∈ℝ^d+1.



    To prove the first part of the result (<ref>) it suffices to show that
    
    r(x) ·( 1 - ∑_ℓ=1^∞ u^T W (W - v u^T W)^ℓ-1 v · x^ℓ) = 1.

    We have
    
    r(x) ·( 1 - ∑_ℓ=1^∞ u^T W (W - v u^T W)^ℓ-1 v · x^ℓ) 
       =   ∑_k=0^∞ u^T W^k v · x^k ·( 1 - ∑_ℓ=1^∞ u^T W (W - v u^T W)^ℓ-1 v · x^ℓ) 
       = ∑_i=0^∞( u^T W^i v - ∑_ℓ=1^i u^T W^i-ℓ v · u^T W (W - v u^T W)^ℓ-1 v) · x^i 
       = ∑_i=0^∞( u^T W^i v - u^T W (W - v u^T W)^i-1 v - ∑_ℓ=1^i-1 u^T W^i-ℓ v · u^T W (W - v u^T W)^ℓ-1 v) · x^i 
       = u^T v + ∑_i=1^∞ u^T ( W^i - W (W - v u^T W)^i-1 - ∑_ℓ=1^i-1 W^i-ℓ v · u^T W (W - v u^T W)^ℓ-1) v · x^i 
       = 1 + ∑_i=1^∞ u^T ( W^i-1 - (W-Wvu^T)^i-1 - ∑_ℓ=1^i-1 W^i-ℓ v · u^T (W - W v u^T)^ℓ-1) W v · x^i ,

    where the final equality uses the fact 
    W (W - v u^T W)^ℓ-1 = W ((I-vu^T)W)^ℓ-1 = (W(I-vu^T))^ℓ-1W = (W - W v u^T)^ℓ-1 W .

    Now it suffices to show that, for all i ≥ 1, we have
    
    W^i-1 = (W-Wvu^T)^i-1 + ∑_ℓ=1^i-1 W^i-ℓ v · u^T (W - W v u^T)^ℓ-1.

    We prove <Ref> by induction on i. For i=1, the equation is trivially true (both sides are the identity).
    Assuming <Ref> holds for a given i ≥ 1, we have
    
    W^i    = W · W^i-1 = ((W-Wvu^T)+Wvu^T) · W^i  
       = ((W-Wvu^T)+Wvu^T) · (W-Wvu^T)^i-1 + W ·∑_ℓ=1^i-1 W^i-ℓ v · u^T (W - W v u^T)^ℓ-1
       = (W-Wvu^T)^i + Wvu^T · (W-Wvu^T)^i-1 + ∑_ℓ=1^i-1 W^i+1-ℓ v · u^T (W - W v u^T)^ℓ-1
       = (W-Wvu^T)^i + ∑_ℓ=1^i W^i+1-ℓ v · u^T (W - W v u^T)^ℓ-1
             - W^i+1-i v · u^T (W-Wvu^T)^i-1 + Wvu^T · (W-Wvu^T)^i-1
       = (W-Wvu^T)^i + ∑_ℓ=1^i W^i+1-ℓ v · u^T (W - W v u^T)^ℓ-1,

    which establishes that <Ref> holds for the next value of i.
    
    To prove <Ref>, we must show that u^T W_0^0 v = 1 and u^T W_0^k v = - u^T W (W-vu^TW)^k-1 v for k ≥ 1.
    We immediately have u^T v = 1.
    By induction we can show that, for all k ≥ 1, we have 
    
    W_0^k = ( [               0           0 ⋯ 0; (W-vu^TW)^k-1 v     (W-vu^TW)^k ]) ∈ℝ^(d+1) × (d+1).

    From this it follows that u^T W_0^k v = - u^T W (W-vu^TW)^k-1 v for all k ≥ 1, as required.
    
    More generally, for all β∈ℝ and k ≥ 1,
    
    W_β^k = ( [                            β^k                          0 ⋯ 0; ∑_ℓ=0^k-1β^k-1-ℓ (W-vu^TW)^ℓ v                    (W-vu^TW)^k ]) ∈ℝ^(d+1) × (d+1).

    Setting β=1 still gives u^T W_1^0 v = 1 and now 
    u^T W_1^k v = 1 - ∑_ℓ=0^k-1 u^T W (W-vu^TW)^ℓ v = 1+ ∑_ℓ=0^k-1u^T W_0^ℓ+1v = ∑_j=0^ku^T W_0^jv
 for k ≥ 1.
    
    We have 1/r(x) = ∑_k=0^∞u^T W_0^k v· x^k and, hence, 
    1/r(x)·(1-x) = ( ∑_k=0^∞u^T W_0^k v· x^k) ·(∑_i=0^∞ x^i) = ∑_ℓ=0^∞ x^ℓ∑_j=0^ℓu^T W_0^j v = ∑_ℓ=0^∞u^T W_1^ℓv· x^ℓ,

    which completes the proof.







§ GENERALIZING THE BINARY TREE APPROACH
 




In this section we prove <Ref>. Compared to <Ref> this attains a weaker approximation to the optimal matrix factorization, but is asymptotically better in terms of space usage.

Our starting point is the binary tree mechanism of <cit.>.
The binary tree mechanism can be viewed as a recursive construction of a matrix factorization; see <Ref> for a precise statement.
A recursion of depth ℓ yields a matrix factorization of size =2^ℓ and an algorithm running in time and space O(ℓ).
The binary tree mechanism is within a constant factor of optimal; we will improve the constant factor to be arbitrarily close to optimal.
Our approach is to combine a recursive construction with the factorization we already developed to prove <Ref>.
We will start with a matrix factorization of size _1 and, after a recursion of depth ℓ, we obtain a matrix factorization of size _1^ℓ. 
The space required is O(ℓ·log^2 _1). Appropriately setting the parameters yields  <Ref>.






    First we define the non-cyclic shift matrix S^()∈{0,1}^× by S_i,j=1  i=j+1 for all i,j∈ [].
    If we multply the non-cyclic shift matrix with the lower-triangular all-ones matrix, it has the effect of zeroing out the diagonal and producing a strictly lower-triangular matrix.
    For example,
    
    S^(6) =
        ([ 0 0 0 0 0 0; 1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0 ])
            and    
        A^(6) S^(6) =  S^(6) A^(6) =
        ([ 0 0 0 0 0 0; 1 0 0 0 0 0; 1 1 0 0 0 0; 1 1 1 0 0 0; 1 1 1 1 0 0; 1 1 1 1 1 0 ]).

    The non-cyclic shift matrix corresponds to the generating function h(x)=x – i.e., S^()=M(h,).
    
    Given matrices M ∈ℝ^n × m and M' ∈ℝ^n' × m', define the Kronecker product 
    
    M ⊗ M'  = (
                [ M_1,1· M' M_1,2· M'         ⋯ M_1,m· M'; M_2,1· M' M_2,2· M'         ⋯ M_2,m· M';         ⋮         ⋮         ⋱         ⋮; M_n,1· M' M_n,2· M'         ⋯ M_n,m· M';           ]) ∈ℝ^(n · n') × (m · m')

    by (M ⊗ M')_i,j = M_i_1,j_1· M'_i_2,j_2 where i = (i_1-1) · n' + i_2 ∈ [n· n'] and j = (j_1-1) · m' + j_2 ∈ [m· m'].
    
    M ⊗ M'  = (
                [     M_0,0· M'     M_0,1· M'             ⋯   M_0,m-1· M';     M_1,0· M'     M_1,1· M'             ⋯   M_1,m-1· M';             ⋮             ⋮             ⋱             ⋮;   M_n-1,0· M'   M_n-1,1· M'             ⋯ M_n-1,m-1· M';               ]) ∈ℝ^(n · n') × (m · m')

    by (M ⊗ M')_i,j = M_i_1,j_1· M'_i_2,j_2 where i = i_1 · n' + i_2 ∈ [n· n'] and j = j_1 · m' + j_2 ∈ [m· m'].
    The Kronecker product has the “mixed-product property” with matrix multiplication: We have (A ⊗ B) (A' ⊗ B') = (AA') ⊗ (BB') whenever the matrix dimensions match so that AA' and BB' are well-defined.
    
 


 §.§ Recursive Matrix Factorization
   
Without further ado, we present the basis of our recursive construction. 
Suppose we are given two matrix factorizations A^(_1) = B_1 C_2 and A^(_2) = B_2 C_2. We will combine these into one matrix factorization for A^(_1 ·_2).

    For B_1 ∈ℝ^_1 ×_1' and B_2 ∈ℝ^_2 ×_2', define
    
    B_1B_2(   I^(_1)⊗ B_2  _∈ℝ^(_1_2) × (_1_2')    S^(_1) B_1 ⊗1^(_2) _∈ℝ^(_1_2) ×_1' ) ∈ℝ^(_1_2) × (_1_2'+_1'),

    where I^(_1) is the _1 ×_1 identity matrix and 1^(_2) is the all-ones column vector of length _2.
    For C_1 ∈ℝ^_1' ×_1 and C_2 ∈ℝ^_2' ×_2, define
    
    C_1C_2( [      I^(_1)⊗ C_2; C_1 ⊗ (1^(_2))^T ]) ∈ℝ^(_1_2'+_1')×(_1_2),

    where I^(_1) is the _1 ×_1 identitiy matrix and (1^(_2))^T is the all-ones row vector of length _2


    Let  B_1,C_1^T ∈ℝ^_1 ×_1' and B_2,C_2^T ∈ℝ^_2 ×_2'.
    Suppose B_1C_1=A^(_1) and B_2C_2=A^(_2), where A^() is the all-ones lower triangular matrix defined in <Ref>.
    Then 
    (B_1B_2) (C_1C_2) = A^(_1_2).

    Furthermore, 
    
    B_1B_2_2→∞^2    = S^(_1)B_1_2→∞^2 + B_2_2→∞^2 ≤B_1_2→∞^2 + B_2_2→∞^2, 
    C_1C_2_1→2^2    = C_1_1→2^2 + C_2_1→2^2 .



    We have
    
    (B_1B_2) (C_1C_2)
               = ( I^(_1)⊗ B_2    S^(_1) B_1 ⊗1^(_2)) ( [      I^(_1)⊗ C_2; C_1 ⊗ (1^(_2))^T ]) 
       = (I^(_1)⊗ B_2) (I^(_1)⊗ C_2) + (S^(_1) B_1 ⊗1^(_2)) (C_1 ⊗ (1^(_2))^T) 
       = (I^(_1) I^(_1)) ⊗ (B_2 C_2) + (S^(_1) B_1 C_1) ⊗ (1^(_2) (1^(_2))^T) 
       = I^(_1)⊗ A^(_2) + (S^(_1) A^(_2)) ⊗ (1^(_2) (1^(_2))^T) 
       = A^(_1_2).

    The last equality is best understood by looking at an example for _1=3 and _2=2:[If you don't like this proof by example, imagine some ellipses (⋯,⋮,⋱) inserted into the matrices and it will look like a general proof. Formally, we can do an index-by-index case analysis.]
    
    I^(3)⊗ A^(2) + ( S^(3) A^(3)) ⊗( 1^(2) (1^(2))^T ) 
       = 
            ([ 1 0 0; 0 1 0; 0 0 1 ]) 
            ⊗([ 1 0; 1 1 ])
            +
            ([ 0 0 0; 1 0 0; 1 1 0 ]) 
            ⊗([ 1 1; 1 1 ]) 
       =
            ([ 1 0 0 0 0 0; 1 1 0 0 0 0; 0 0 1 0 0 0; 0 0 1 1 0 0; 0 0 0 0 1 0; 0 0 0 0 1 1 ])
            + 
            ([ 0 0 0 0 0 0; 0 0 0 0 0 0; 1 1 0 0 0 0; 1 1 0 0 0 0; 1 1 1 1 0 0; 1 1 1 1 0 0;   ]) 
       =
            ([ 1 0 0 0 0 0; 1 1 0 0 0 0; 1 1 1 0 0 0; 1 1 1 1 0 0; 1 1 1 1 1 0; 1 1 1 1 1 1 ]) 
       = A^(6).

    
    Now we have
    
    B_1B_2_2→∞^2
               = max_i ∈ [_1_2]∑_j ∈ [_1_2'+_1'] (B_1B_2)_i,j^2 
       = max_i ∈ [_1_2]∑_j' ∈ [_1_2'] (I^(_1)⊗ B_2)_i,j'^2 + ∑_j”∈ [_1] (S^(_1) B_1 ⊗1^(_2))_i,j”^2 
       = max_i_1 ∈ [_1], i_2 ∈ [_2]∑_j_1' ∈ [_1], j_2' ∈ [_2'] (I^(_1)_i_1,j_1· (B_2)_i_2,j_2')^2 + ∑_j”∈ [_1] ((S^(_1) B_1)_i_1,j”·1^(_2))_j”)^2 
       = max_i_1 ∈ [_1], i_2 ∈ [_2]∑_ j_2' ∈ [_2'] (B_2)_i_2,j_2'^2 + ∑_j”∈ [_1] (S^(_1) B_1)_i_1,j” ^2 
       = max_i_1 ∈ [_1]∑_ j_2' ∈ [_2'] (B_2)_i_2,j_2'^2 + max_ i_2 ∈ [_2]∑_j”∈ [_1] (S^(_1) B_1)_i_1,j” ^2 
       = B_2_2→∞^2 + S^(_1)B_1_2→∞^2.

    Similarly, C_1C_2_1→2^2 = C_1_1→2^2 + C_2_1→2^2.
    Finally, 
    S^(_1)B_1_2→∞^2 = max_i ∈ [_1]∑_j ∈ [_1'] (S^(_1)B_1)_i,j^2  = max_i ∈ [_1-1]∑_j ∈ [_1'] (B_1)_i,j^2 ≤B_1_2→∞^2.



<Ref> can be applied recursively. We start with some base factorization B_1 C_1 = A^(_1) and inductively define a factorization satisfying B_ℓ C_ℓ = A^(_1^ℓ) and 
    B_ℓ_2→∞·C_ℓ_1→2≤ℓ·B_1_2→∞·C_1_1→2.



    Let B_1, C_1^T ∈ℝ^_1 ×'_1. 
    Following <Ref>, for ℓ≥ 2, define 
    
    B_ℓ = B_1B_ℓ-1    and    
            C_ℓ = C_1C_ℓ-1.

    Then, for all ℓ≥ 1, we have B_ℓ, C_ℓ^T ∈ℝ^_ℓ×'_ℓ, where _ℓ_1^ℓ and 
    '_ℓ_1' ·∑_k=0^ℓ-1_1^k = _1' ·_1^ℓ-1/_1-1

    If B_1 C_1 = A^(_1), then B_ℓ C_ℓ = A^(_ℓ) for all ℓ≥ 1.
    Furthermore,
    
    B_ℓ_2→∞   ≤√(ℓ)·B_1_2→∞,
    C_ℓ_1→2   = √(ℓ)·C_1_1→2.



    Perform induction on ℓ using <Ref> with B_2=B_ℓ-1 and C_2=C_ℓ-1.




 §.§ Recursive Algorithm


Now we show that the recursive matrix factorization in <Ref> corresponds to an efficient algorithm.
The recursive construction naturally corresponds to a recursive algorithm.
That is, for streaming multiplication by B_ℓ we must call a base algorithm for B_1 and recursively call the algorithm for B_ℓ-1.

Recall the streaming multiplication setting: We receive the rows of Z ∈ℝ^'_ℓ× one by one as a stream and we must output the rows of B_ℓ· Z ∈ℝ^_ℓ× as a stream.
Note that our matrix B_ℓ is not a lower triangular matrix; in fact, it is not even square. 
This makes the problem setting not entirely straightforward, as we do not have a one-to-one correspondence between streaming inputs and outputs.

In our application, the input Z is independent random noise. In particular, the order of the rows is unimportant. 
The time when we read each row is also not important as the noise can be generated as needed. The constraint is simply that we can only read each row once. (We can store the rows for later use, but that requires memory, so we want to avoid this.) 
In other words, the input is a stream of independent random noise, which we can read at any time, but we cannot “rewind” it.




Our recursive streaming algorithm is presented as <Ref>.
We will instantiate the base streaming algorithm 𝒜_1 with <Ref>.
The recursive algorithm runs many copies of 𝒜_1.
The streaming input is simply split among these copies. For simplicity, we do not belabor the precise order in which the input stream is read.



    For ℓ≥ 1, let B_ℓ,C_ℓ^T ∈ℝ^_ℓ×_ℓ' be as in <Ref>.
    Let 𝒜_1 be a streaming algorithm that takes as input a stream Z_0,·,Z_1,·,⋯,Z__ℓ-1'-1,·∈ℝ^1 × and outputs a stream Z_0,·,Z_1,·,⋯,Z__ℓ-1-1,·∈ℝ^1 × such that Z = B_1 Z.
    Let 𝒜_ℓ denote <Ref> instantiated with the base algorithm 𝒜_1 and recursion depth ℓ≥ 1.
    
    Then 𝒜_ℓ is a streaming algorithm that takes as input a stream Z_0,·,Z_1,·,⋯,Z_'_ℓ-1,·∈ℝ^1 × (in an arbitrary – but fixed – order) and outputs a stream Z_0,·,Z_1,·,⋯,Z__ℓ-1,·∈ℝ^1 × such that Z = B_1 Z.

    Moreover, the space required by 𝒜_ℓ is at most ℓ times the space required by 𝒜_1. And the time required by 𝒜_ℓ is _1^ℓ+1-1/_1-1≤ O(_ℓ) times the time required by 𝒜_1.


    We may inductively assume that the lemma holds for ℓ-1 in place of ℓ. That is, we assume 𝒜_ℓ-1 is a streaming algorithm that takes as input a stream Z_0,·,Z_1,·,⋯,Z_'_ℓ-1,·∈ℝ^1 × (in an arbitrary – but fixed – order) and outputs a stream Ž_0,·,Ž_1,·,⋯,Ž__ℓ-1-1,·∈ℝ^1 × such that Ž = B_ℓ-1 Z.
    Let 
    
    Z' = ( [    z'_0;    z'_0;       ⋮;    z'_0;    z'_1;    z'_1;       ⋮;    z'_1;       ⋮; z'__1-1; z'__1-1;       ⋮; z'__1-1 ]) = ( [    z'_0;    z'_1;       ⋮; z'__1-1 ]) ⊗1^(_ℓ-1)∈ℝ^_1_ℓ-1×

    and
    
    Z” = ([       z”_0,0;       z”_0,1;            ⋮;    z”_0,_ℓ-1;       z”_1,0;       z”_1,1;            ⋮;    z”_1,_ℓ-1;            ⋮;    z”__1-1,0;    z”__1-1,1;            ⋮; z”__1-1,_ℓ-1;  ])  ∈ℝ^_1_ℓ-1×.

    We have Z = Z' + Z”. We want to show that 
    
    Z = B_ℓ Z = (B_1B_ℓ) Z = ( I^(_1)⊗ B_2    S^(_1) B_1 ⊗1^(_2)) ( [   Z_0:_1'_ℓ-1,·; Z__1'_ℓ-1:'_ℓ,· ]),

    where Z__1'_ℓ-1:'_ℓ,·∈ℝ^'_1 × and Z_0:_1'_ℓ-1,·∈ℝ^_1 '_ℓ-1× are a partition of the rows of Z ∈ℝ^'_ℓ×. (Note that '_ℓ = _1'_ℓ-1+'_1. The notation Z_i:j,· denotes the submatrix of Z formed by rows i, i+1, ⋯, j-1 and all columns.)
    Thus it suffices to show that Z' = (S^(_1) B_1 ⊗1^(_2)) Z__1'_ℓ-1:'_ℓ,· and Z” = (I^(_1)⊗ B_2) Z_0:_1'_ℓ-1,·.
    By our inductive assumption on 𝒜_ℓ-1, for all i ∈ [_1], 
    ([    z”_i,0;    z”_i,1;         ⋮; z”_i,_ℓ-1 ]) = B_ℓ· Z_i_ℓ-1':(i+1)_ℓ-1',·.

    Thus 
    
    Z” = ([             B_ℓ· Z_0:_ℓ-1',·;        B_ℓ· Z__ℓ-1':2_ℓ-1',·;                            ⋮; B_ℓ· Z_(_1-1)_ℓ-1':_1_ℓ-1',· ])  = (I^(_1)⊗ B_ℓ-1) · Z_0:_1_ℓ-1,· ,

    as required.
    By our assumption about the base algorithm 𝒜_1, we have 
    
    ( [  z'_1;  z'_2;     ⋮; z'__1 ]) = B_1 · Z__1'_ℓ-1:_ℓ',·∈ℝ^_1 ×.

    (Note that z'__1 is never actually used in the algorithm.)
    With z'_0=(0^(m))^T, the non-cyclic shift S^(_1) gives
    
    ( [    z'_0;    z'_1;       ⋮; z'__1-1 ]) = S^(_1)· B_1 · Z__1'_ℓ-1:_ℓ',·∈ℝ^_1 ×.

    Thus 
    Z' = ( [    z'_0;    z'_1;       ⋮; z'__1-1 ]) ⊗1^(_ℓ-1) = ( S^(_1)· B_1 · Z__1'_ℓ-1:_ℓ',·) ⊗1^(_ℓ-1) = ( ( S^(_1)· B_1 ) ⊗1^(_ℓ-1)) · Z__1'_ℓ-1:_ℓ',·,

    which establishes the correctness of the algorithm.
    
    The space usage of 𝒜_ℓ is is comprised of that of 𝒜_1 plus that of 𝒜_ℓ-1. (Note that, although many copies of 𝒜_ℓ-1 are run, they are run sequentially, so the space usage does not add up.) Thus, by induction, the space usage of 𝒜_ℓ is ℓ times that of 𝒜_1.
    
    The runtime of 𝒜_ℓ is comprised of that of 𝒜_1 plus _1 times that of 𝒜_ℓ-1. By induction, this gives a runtime of 1 + _1 + _1^2 + ⋯ + _1^ℓ = _1^ℓ+1-1/_1-1 = '_ℓ/'_1 times that of 𝒜_1, as required.




 §.§ Combining Recursion with 


The algorithm in <Ref> requires a base factorization and a corresponding algorithm.
We can instantiate this with the rational function approximation algorithm from <Ref>.
This combination yields the following result.


    Let , ℓ≥ 1, _1 ≥ 5, and d ≥ 2+(12+4log_1/π)^2 be integers.
    Let =_1^ℓ and ' = _1^ℓ+1-_1/_1-1.
    Then there exist matrices B,C^T ∈ℝ^×' and a streaming algorithm 𝒜 satisfying the following.
    
        
  * Valid Matrix Factorization:
        BC=A^().
        
  * Near-Optimality:
        
    (B,C) ≤ℓ·(√(𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(_1)) + 16√(_1)exp(-π/2√(d-2)))^2,

        where 𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(_1) = 1+∑_k=1^_1-1( 2^-2k2k  k)^2 ≤ 1 + 0.57722+log(_1)/π.
        
  * Valid Algorithm:
        The algorithm 𝒜 takes as input a stream Z_0,·,Z_1,·,⋯,Z_'-1,·∈ℝ^1 × in some arbitrary-but-fixed order (but only reading each Z_i,· once) and outputs a stream Z_0,·,Z_1,·,⋯,Z_-1,·∈ℝ^1 × such that Z = B · Z.
        
  * Efficient Algorithm:
        The algorithm 𝒜 uses space O(ℓ d ) and has total runtime O( d ).
    


    <Ref> provides a base factorization B_1, C_1 ∈ℝ^_1 ×_1 and a base algorithm 𝒜_1 (<Ref>) for computing B_1Z in a streaming manner.
    We are guaranteed that
    
    B_1_2 →∞   ≤B_1^*_2 →∞ + 16 √(_1)exp(-π/2√(d-2)),
    C_1_1 → 2   ≤C_1^*_1 → 2 + 16 √(_1)exp(-π/2√(d-2)),

    where B_1^*=C_1^*=M(f,_1) for f(x)=1/√(1-x).
    By <Ref>, B_1^*,C_1^* is the optimal matrix factorization. In particular, 
    B_1^*_2 →∞ = C_1^*_1 → 2 = √(𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(_1)) = √(1+∑_k=1^_1-1( 2^-2k2k  k)^2)≤√(1 + γ+log(_1)/π),

    where γ≤ 0.57722 is the Euler-Mascheroni constant.
    The algorithm 𝒜_1 runs in space and time (per iteration) O(d).
     
    Now we apply the recursive algorithm with this base construction.
    By <Ref>, we obtain a matrix factorization B_ℓ,C_ℓ∈ℝ^×' such that B_ℓ C_ℓ = A^(n),
    
    B_ℓ_2→∞≤√(ℓ)B_1_2→∞≤√(ℓ)( √(𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(_1)) + 16 √(_1)exp(-π/2√(d-2)) ) ,

    and
    
    C_ℓ_1→2 = √(ℓ)C_1_1→2≤√(ℓ)(√(𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(_1)) + 16 √(_1)exp(-π/2√(d-2))).


    By <Ref>, <Ref> is a valid streaming algorithm for computing Z = B_ℓ Z and runs in space O(ℓ d ) and total time O(n d ), as required. 


Setting parameters in <Ref> yields our result.


    Let , _1, ≥ 1 be integers.
    Then there exist matrices B,C^T ∈ℝ^×' (for some '≤ O()) and a streaming algorithm 𝒜 satisfying the following.
    
        
  * Valid Matrix Factorization:
        BC=A^().
        
  * Near-Optimality:
        
    (B,C) ≤𝖮𝗉𝗍(n) ·(1 + O(1/log(_1))) + O(log(_1)),

        where 𝖮𝗉𝗍(n) = inf{(B_*,C_*) : B_*C_*=A^(n)} = log n/π± O(1) is the optimal value of the matrix factorization objective over all possible factorizations.
        
  * Valid Algorithm:
        The algorithm 𝒜 takes as input a stream Z_0,·,Z_1,·,⋯,Z_'-1,·∈ℝ^1 × in some arbitrary-but-fixed order (but only reading each Z_i,· once) and outputs a stream Z_0,·,Z_1,·,⋯,Z_-1,·∈ℝ^1 × such that Z = B · Z.
        
  * Efficient Algorithm:
        The algorithm 𝒜 uses space O(log(n)log(n_1)) and has total runtime O(log(n_1)^2 ).
    

To obtain <Ref> from <Ref> we simply need to set _1 to a value that is superconstant (so that 1/O(log(_1))=o(1)) but not too large (so that  O(log(n)log(n_1))=O(log(n))). For example, we can set _1=Θ(log()).

    Let _1 ≥ 5, d = ⌈ 2+(12+4log_1/π)^2 ⌉, ℓ = ⌈log()/log(_1)⌉, and ' = _1^ℓ+1-_1/_1-1 be integers.
    
    Then ≤_1^ℓ < ·_1. We will provide a matrix factorization of size _1^ℓ, which can be truncated to one of size .

    <Ref> provides the matrix factorization and algorithm. It gives the near-optimality guarantee
    
    B_ℓ_2→∞·C_ℓ_1→2≤ℓ·(√(𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(_1)) + 16√(_1)exp(-π/2√(d-2)))^2,

    where 𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(_1) = 1+∑_k=1^_1-1( 2^-2k2k  k)^2 ≤ 1 + γ+log(_1)/π.
    It only remains for us to simplify this guarantee.    
        
    We assume d ≥ 2 + (2log(16) + 3log(_1)/π)^2 so that 16 √(_1)exp(-π/2√(d-2)) ≤ 1/_1.

    Thus
    
    B_ℓ_2→∞C_ℓ_1→2   ≤ℓ·(√(𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(_1)) + 1/_1)^2 
       = ℓ·𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(_1) + ℓ/_1·(2√(𝖮𝗉𝗍𝖫𝖳𝖳𝗈𝖾(_1)) + 1/_1) 
       ≤ℓ·( 1 + γ + log(_1)/π) + ℓ/_1·(2√(_1) + 1/_1 ) 
       =  ℓlog(_1)/π  + ℓ· (1+γ/π) + ℓ·(2/√(_1) + 1/_1^2 ) 
       =  log(_1^ℓ)/π  + ℓ·(1+γ/π+ 2/√(_1) + 1/_1^2) 
       ≤log()+log(_1)/π  + ℓ·(1+γ/π+ 2/√(5) + 1/5^2) 
       ≤log()/π  + 3ℓ + log(_1)/π.

    Now we invoke the lower bound of <cit.> or <cit.> (see <Ref>), which tells us 
    𝖮𝗉𝗍() ≥log()/π -1,

    where 𝖮𝗉𝗍() is the optimal value of the objective function over all factorizations of size .
    Thus
    
    B_ℓ_2→∞C_ℓ_1→2   ≤𝖮𝗉𝗍() + 1 + 3ℓ + log(_1)/π
       ≤𝖮𝗉𝗍() + 4 + 3 log()/log(_1) + log(_1)/π
       ≤𝖮𝗉𝗍() + 4 + 3 1+π𝖮𝗉𝗍()/log(_1) + log(_1)/π
       = 𝖮𝗉𝗍() ·(1 + 3π/log(_1)) + 4 + 3/log(_1) + log(_1)/π .

    By <Ref>, the space usage of the algorithm is O(ℓ d ) = O(log()/log(_1) (log(_1))^2 ) = O(log()log(_1)) and the total runtime is O( d ) = O(log(_1)^2 ).




§ NUMERICAL LOWER BOUND ON OPTIMAL PERFORMANCE

In this section, we develop numerical lower bounds on the optimal matrix factorizations of the all-ones lower triangular matrix A^(n) for various classes of matrices. More precisely, given the sequence length n, we write down a semidefinite program that provides a lower bound on the lowest achievable error over a specific class of correlation matrices (for example, Toeplitz matrices or triangular matrices). 

In general, we consider mechanisms whose matrices can be written as 
C=c=∑_i _i c_i where c ∈^m denotes the parameters of the mechanism (to be optimized) and : ^m ↦^n × n is an linear operator and _i ∈^n × n are “basis” matrices. Any linear subspace of the space of n × n matrices (including lower triangular matrices and Toeplitz matrices) can be expressed in this form.

All results in this section derive from the following fundamental result which is a straightforward consequence of weak duality:


Let  denote the n × n lower triangular matrix with all entries equal to 1. 
For any c ∈^m, define the objective for the correlated mechanism with correlation matrix c by

    Fc=(c,c) = c_1 → 2^2c_2 →∞^2.

Further, for any set ℭ⊆^m, we have 

    min_c ∈ℭ Fc≥1/min_Γ_0, …, Γ_N ∈^n_+
    ∑_i _i^T Γ_i _i=1max_c_1↦ 2=1cc∑_i Γ_i

where ^n_+ denotes the space of symmetric positive semidefinite n × n matrices, _i denotes the i-th column of  and AB=traceAB denotes the standard inner product between symmetric matrices A, B ∈^N.

Furthermore, if the inner maximization can be relaxed to a concave maximization problem, i.e, there exists a compact convex set 𝒞⊂^N such that

    max_c ∈ℭ
    c_1↦ 2=1ccΓ = max_S ∈𝒞SΓ  ∀Γ∈^N_+

then the inequality in <ref> is an equality. 


Since the objective is invariant to scaling c by any positive constant, i.e., Fα c=Fc for all α∈, α > 0, we can set α=1/c_1 → 2 so that the optimization problem becomes

    min_c ∈ℭ
    c_1 → 2=1c_2→inf^2=min_c ∈ℭ
    c_1 → 2=1max_ic_i^2


From Lemma <ref>, we have that

    M = 1/min_Γ∈^N_+ 
    ΓΓM

for any M ∈^n_++. Applying this to M=cc and =_i (for i=0, …, n) and plugging this into (<ref>), we obtain

    max_ic_i^2 = max_i max_Γ_i ≽ 0 
    _i^T Γ_i _i=11/ccΓ_i=1/min_i min_Γ_i ≽ 0 
    _i^T Γ_i _i=1ccΓ_i

The denominator can be rewritten as

    min_ν∈Δ^n∑_i ν_imin_Γ_i ≽ 0 
    _i^T Γ_i _i=1ccΓ_i=min_ν∈Δ^n∑_i min_Γ_i ≽ 0 
    _i^T Γ_i _i=ν_iccΓ_i
       =min_Γ_0, …, Γ_N ∈^n_+
    ∑_i _i^T Γ_i _i=1cc∑_i Γ_i

where Δ^n is the simplex in n dimensions. Thus, the optimization problem (<ref>) can be rewritten as

    min_c_1 → 2≤ 11/min_Γ_0, …, Γ_N ∈^n_+
    ∑_i _i^T Γ_i _i=1cc∑_i Γ_i = 1/max_c_1 → 2≤ 1min_Γ_0, …, Γ_n ∈^n_+
    ∑_i _i^T Γ_i _i=1cc∑_i Γ_i

By weak duality, we have that the optimal value is bounded below by

    1/min_Γ_0, …, Γ_n ∈^n_+
    ∑_i _i^T Γ_i _i=1max_c_1 → 2≤ 1cc∑_i Γ_i

which gives <ref>.

Furthermore, if the inner maximization can be rewritten as a concave maximization problem, the overall problem is a convex-concave problem and strong duality holds by the Von-Neumann minimax theorem <cit.>, so that the order of the min and max can be interchanged without changing the optimal value. Hence, the lower bound is tight.




Let 𝒥 = {c: c ∈ℭ}. A tight lower bound in <ref> is achieved in the following 3 cases:

    
  * 𝒥 is the set of all lower triangular matrices.
    
  * 𝒥 is the set of all lower triangular Toeplitz matrices.
    
  * 𝒥 is the set of all lower triangular Toeplitz matrices formed from a degree 1 constant recurrent sequence.



1. If 𝒥 is the set of all lower triangular matrices, for any 𝕁∈𝒥, letting Δ=𝕁𝕁, we have that 

    max_𝕁_1 → 2≤ 1, 𝕁∈𝒥𝕁𝕁Γ≤max_Δ≽ 0 
    diagΔ≤ 1ΔΓ

since for any feasible 𝕁, Δ=𝕁𝕁 satisfies the constraints listed in the optimization problem on the right hand side (note that diagΔ refers to the vector of diagonal elements of the matrix Δ). 

Conversely, given any feasible solution Δ satisfying the constraints, one can compute its upper triangular Cholesky factorization (the regular Cholesky factorization multiplied by the matrix with anti-diagonals equal to 1), to obtain 𝕁 that satisfies the constraint 𝕁_1 → 2≤ 1 and Δ=𝕁𝕁. Hence the two optimization problems have equal optimal values, and hence the lower bound in theorem <ref> is tight.

2. If 𝒥 is the set of all n × n lower triangular Toeplitz matrices, let 𝕁c denote the lower triangular Toeplitz matrix whose first column is c. Then, we have 

    c_1→ 2 = c, 𝕁c=∑_i c_i H_i

where H_i is the matrix with 0s everywhere except the i-th subdiagonal equal to 1. Thus

    max_𝕁c_1 → 2≤ 1, 𝕁c∈𝒥𝕁𝕁Γ = max_c = 1𝕁c𝕁cΓ

and the right hand side is a standard trust region subproblem with a quadratic objective and a single quadratic constraint, that can be solved exactly via an SDP relaxation by the S-lemma <cit.>, i.e, there exists an affine map Q: ^n →^n such that the right hand side equals max_C ≽ 0, traceC=1QΓ C. Hence, the lower bound in <ref> is tight.


3. If 𝒥 is the set of all n × n lower triangular Toeplitz matrices formed from a degree 1 constant recurrent sequence, we have that 𝕁=𝕁c where c_i = αθ^i. where α, θ∈ are arbitrary real numbers. Then, the problem

    max_c = 1𝕁c𝕁cΓ=max_c𝕁c𝕁cΓ/cc

can be written as

    max_θ, αα^2 pθ;Γ/α^2 qθ = max_θpθ;Γ/qθ

where p, q are polynomials in θ and the coefficients of q are independent of Γ while the coefficients of p are affine functions of Γ. Using the fact that 
    pθ-β qθ≤ 0  pθ-β qθ is a sum of squares polynomial

we can write the maximization problem as being equal to

    min_pθ;Γ-β qθ is a sum of squares polynomialβ

Thus, the inner maximization can be written as a concave maximization problem of the form required for lower bound tightness in Theorem <ref> <cit.>.






 §.§ Lower bound over all lower triangular matrices





 §.§ Lower bound over Toeplitz matrices


Let C ∈^N × N denote the Toeplitz lower triangular correlation matrix and c denote its first column. Its clear that the entries of C are affine functions of C and we use C=c=∑_i _i c_i to denote that where : ^N ↦^N × N is an linear operator and 

    _i ∈^N × N
    = [ 0_N-i, N-i   0_N-i, i; 𝕀_N-i, N-i   0_N-i, i ]

𝕀_k is the k × k identity matrix and 0_m, n is the m × n matrix of 0s. The inf-norm error of the prefix sum estimate is given by

    Fc=c^2 c^2

where · denotes the Euclidean ℓ_2 norm in ^N and ∈^N is the vector with all entries equal to 1.




min_c ∈^N Fc is bounded below by the optimal value of the following semidefinite prgram:



    max_t ∈, C ∈^N      t 
    Subject to    [                   t                    ;                     ∑_i, j=1^N_i_j C_ij ]≽ 0 
        C ≽ 0 
       trace[C] = 1




Furthermore, the choice of c=c^⋆ developed in <cit.> which is defined by c^⋆_0 = 0 and c^⋆_i+1 = c^⋆_i1 - 1/2(i+1) is such that Fc^⋆ equals the optimal value of the optimization problem (<ref>)




 §.§ Lower bound over degree 1 constant recurrent sequence matrices


Further, if we impose additional constraints on c, for example, that c is an order d constant recurrence sequence, we can get tighter lower bounds on the optimal performance. In particular, for d=1, we have the following result:

Let ℂ_d denote the set of constant recurrent sequences of order d. Then, min_c ∈ℂ_1 Fc is bounded below by the optimal value of the semidefinite program


    min_t ∈, C ∈^N      t 
    Subject to    [                     t                      ;                       ∑_i, j=0^N-1_i_j C_ij ]≽ 0 
        C ≽ 0 
       trace[C] = 1 
        C_i,k-i = C_0, k  ∀ i ∈{0, …, k}, k ∈{1, …, N-1}






[Lower bound on matrix fractional function]
For any ∈^n and M ∈^n_++, we have

    M = 1/min_Γ∈^n_+ 
    Γ = 1ΓM



Let ∈^n be an arbitrary vector. 
Now, consider any positive definite matrix M ∈^n_+ and eigendecomposition M=∑_i λ_i u_iu_i for some orthogonal basis {u_i}. Further, expressed in this basis, suppose that =∑_i σ_i u_i. 

Then, we have

    M^-1 = ∑_i 1/λ_iu_i^2 = ∑_i σ_i^2/λ_i

Let MΓ denote the inner product in ^n, i.e., MΓ=traceMΓ.
Further, for any matrix Γ∈^n_+, we have

    M^-1MΓ = ∑_i σ_i^2/λ_i∑_i λ_i u_iΓ u_i =∑_i σ_i^2/λ_i∑_i λ_i Γ^1/2u_i^2
       ≥∑_i σ_iΓ^1/2u_i^2 ≥∑_i σ_i Γ^1/2 u_i^2 = ^T Γ

where the first inequality is an application of the Cauchy-Schwartz inequality and the second inequality is an application of the triangle inequality. 

Thus, for any M, Γ∈^n_+, we have

    M^-1≥^T Γ/MΓ

Further, choosing Γ=MM, the inequality is replaced with equality. Hence, we have 

    M^-1 = max_Γ≽ 0^T Γ/MΓ = max_Γ≽ 0 
    ^T Γ=11/MΓ





 §.§ Symmetrization

    
    The lower triangular matrix A is somewhat messy to work with. Instead we will work with a skew-symmetric matrix
    
    AA-A^T = ( [  0 -1 -1  ⋯ -1;  1  0 -1  ⋯ -1;  1  1  0  ⋯ -1;  ⋮  ⋮  ⋮  ⋱  ⋮;  1  1  1  ⋯  0 ]) ∈{-1,0,1}^m × m .

    We can obtain A from A by A = A + 1⃗1⃗^T + I/2, where 1⃗ is the all-ones vector (and, hence 1⃗1⃗^T is the all-ones matrix) and I is the identity matrix.
    Furthermore, given a factorization A=BC with B,C^T ∈ℝ^m × d, we can obtain a factorization A=BC, where
    
    B    = 1/√(2)( B | 1⃗ |  I ) ∈ℝ^m × (d+1+m),
    
            C    = 1/√(2)( [    C; 1⃗^T;    I ]) ∈ℝ^(d+1+m) × m.

    This reduction preserves the objective:
    
    B_2 →∞^2    = 1/2B_2 →∞^2 + 1, 
    C_1 → 2^2    = 1/2C_1 → 2^2 + 1.

    Note that, if B^T = C, then B^T=C.
    
    Mathias <cit.> gives an exact expression for the objective value of the optimal factorization of this skew-symmetric matrix:
    
    γ_2(A)    inf{B_2 →∞·C_1 → 2 : A=BC}
       = 1/m∑_j=1^m |1/tan(2j-1/2mπ)| = 2/m∑_j=1^⌊ m/2 ⌋1/tan(2j-1/2mπ)

    
    

 §.§ Combining Factorizations

    
    Given matrices M ∈ℝ^n × m and M' ∈ℝ^n' × m', define the Kronecker product 
    
    M ⊗ M'  = (
                [ M_1,1· M' M_1,2· M'         ⋯ M_1,m· M'; M_2,1· M' M_2,2· M'         ⋯ M_2,m· M';         ⋮         ⋮         ⋱         ⋮; M_n,1· M' M_n,2· M'         ⋯ M_n,m· M';           ]) ∈ℝ^(n · n') × (m · m')

    by (M ⊗ M')_i,j = M_i_1,j_1· M'_i_2,j_2 where i = (i_1-1) · n' + i_2 ∈ [n· n'] and j = (j_1-1) · m' + j_2 ∈ [m· m'].
    The Kronecker product has the property that (A ⊗ B) (A' ⊗ B') = (AA') ⊗ (BB').
    
    Suppose we are given two matrix factorizations A^(m_1) = B^(m_1)C^(m_1) and A^(m_2) = B^(m_2)C^(m_2). We will combine these into one matrix factorization for A^(m_1 · m_2).
    
    For  B ∈ℝ^m × d and B' ∈ℝ^m' × d', define
    
    BB'( B ⊗1⃗_m'| I_m ⊗ B' ) ∈ℝ^(m · m') × (d + m · d'),

    where 1⃗_m' is a m' × 1 all-ones column vector and I_m is the m × m identity matrix.
    For  C ∈ℝ^d × m and C' ∈ℝ^d' × m', define
    
    CC'( [ C ⊗1⃗_m'^T;   I_m ⊗ C' ]) ∈ℝ^(d+ m · d') × (m · m'),

    where 1⃗_m'^T is the 1 × m' all-ones row vector and I_m is the m × m identity matrix.
    
    We have 
    
    (BB') (CC')    = (B ⊗1⃗_m') (C ⊗1⃗m'^T) + (I_m ⊗ B') (I_m ⊗ C') 
       = (BC) ⊗ (1⃗_m'1⃗_m'^T) + I_m ⊗ (B'C'), 
    BB'_2 →∞^2    ≤B_2 →∞^2 + B'_2 →∞^2, 
    CC'_1 → 2^2    ≤C_1 → 2^2 + C'_1 → 2^2.

    
    Assume B^(m_1), (C^(m_1))^T ∈ℝ^m_1 × d_1 and B^(m_2), (C^(m_2))^T ∈ℝ^m_2 × d_2 satisfy A^(m_1) = B^(m_1)C^(m_1) and A^(m_2) = B^(m_2)C^(m_2).
    Let m = m_1 · m_2 and d=d_1+m_1· d_2.
    We will construct a factorization A^(m) = B^(m)C^(m) with B^(m), (C^(m))^T ∈ℝ^m × d.
    Namely, B^(m) = B^(m_1)B^(m_2) and C^(m) = C^(m_1)C^(m_2).
    
    
    

 §.§ Iterative Construction

    
    Suppose we have a base factorization A^1 = B^1C^1∈{-1,0,1}^m_1 × m_1, where B^1, (C^1)^T ∈ℝ^m_1 × d_1.
    Then we can iteratively construct larger and larger factorizations: 
    B^k = B^1B^k-1∈ℝ^m_k × d_k      and      C^k = C^1C^k-1∈ℝ^d_k × m_k,

    where m_k = m_1^k and d_k = d_1 ·m_1^k-1/m_1-1. 
    In terms of the objective <ref>, we have
    
    B^k_2 →∞ = √(k)B^1_2 →∞      and      C^k_1 → 2 = √(k)C^1_1 → 2.

    We can then de-symmetrize this to obtain a factorization of the lower triangular matrix A=BC ∈{0,1}^m_1^k × m_1^k with B,C^T ∈ℝ^m_1^k × (d_1 m_1^k-1/m_1-1+1+m_1^k) and
    
    B_2 →∞^2 = k/2B^1_2 →∞^2 + 1
               and C_1 → 2^2 = k/2C^1_1 → 2^2 + 1.

    
    In terms of computational efficiency, we can generate m_k = m_1^k samples in a streaming fashion by running k parallel copies of the sampling algorithm for the base case. Note that we will restart these parallel algorithms many times.  
    In terms of space complexity the base algorithm requires d_1 “registers” and thus the recursive algorithm requires d_1 · k = d_1 ·log(m_k)/log(m_1) registers, which is logarithmic in the desired number of outputs m_k. 
    
    

  §.§.§ Instantiation with 2 × 2 Base Factorization

    
    We instantiate the iterative construction with the base algorithm given by
    
    B^1 = C^1 = 1/√(2)([  1 -1;  1  1 ]).

    This is the optimal factorization <ref> for dimension 2.
    This gives B^k, (C^k)^T ∈ℝ^2^k × (2^k+1-2) with 
    B^k_2 →∞ = C^k_1 → 2 = √(k).

    We can de-symmetrize as in <ref> to obtain B^k, (C^k)^T ∈ℝ^2^k × (3 · 2^k - 1) with 
    B^k_2 →∞ = C^k_1 → 2 = √(1/2 k + 1).

    If we have a particular size m we are aiming for we can set k=⌈log_2 m ⌉ and discard superfluous rows/columns. Thus the value of our objective is
    
    B_2 →∞·C_1 → 2≤1/2⌈log m/log 2⌉ + 1 ≤log m/2log 2 + 3/2.

    This bound is similar to that of <cit.>. We obtain the leading term log m/2 log 2, while <cit.> obtain the leading term log m/π. The difference is a factor of π/2log 2≈ 2.266. In exchange, the advantage of our method is that it is computationally efficient to sample.
    Figure <ref> shows how this mechanism compares.
    
    
    

  §.§.§ Instantiation with 3 × 2 Base Factorization

    
    Using a larger base factorization (rather than 2 × 2) should achieve better constants.
    
    Here is something for 3 × 2 which is nearly optimal:
    
    B^1 = ( [    0   -1;    1 -1/2;    1  1/2 ])
                  and      C^1 = ( [    1  1/2 -1/2;    0    1    1 ]).

    We have B^1_2 →∞ = C^1_1 → 2 = √(5/4) and, thus, applying the iterative construction we get following objective value
    
    B_2 →∞·C_1 → 2≤1/2·5/4·⌈log m/log 3⌉ + 1 ≤5log m/8log 3 + 13/8.

    Comparing leading constants: For a 2 × 2 base we have 1/2 log 2≈ 0.721; for a 3 × 2 base we have 5/8 log 3≈ 0.569; for the analytic construction <cit.> we have 1/π≈ 0.318. The suboptimality ratio for 3 × 2 is 5/8 log 3 / 1/π≈ 1.787.
    Figure <ref> shows how this construction compares to the binary tree, as well as the near-optimal constructions.
    
    
    
    
    

  §.§.§ Instantiation with Larger Base Factorizations

    We can instantiate the iterative construction with larger and larger base factorizations.
    Equation <ref> gives the optimal objective value for the base factorization.[The bound given by <cit.> is not constructive, but numerically we can obtain an optimal factorization from the singular value decomposition of the skew-symmetric matrix. I.e., if A=U Σ V, where U and V are unitary and Σ is diagonal, then B=U√(Σ) and C=√(Σ)V gives an optimal factorization.]
    
    In general, if we start with an m_1 × m_1 base factorization, the objective value for factorizing A^(m)=BC is
    
    B_2 →∞·C_1 → 2≤ 1 + 1/2⌈log m/log m_1⌉·γ_2(A^(m_1)) = 1 +  ⌈log m/log m_1⌉·1/m_1∑_j=1^⌊ m_1/2 ⌋1/tan(2j-1/2m_1π)

    We focus on the leading constant γ_2(A^(m_1))/log m_1 = 1/m_1 log m_1∑_j=1^⌊ m_1/2 ⌋1/tan(2j-1/2m_1π). As m_1 increases this leading constant approaches the optimal 1/π.
    However, it approaches optimality very slowly. Table <ref> shows that even for a base factorization of size m_1=10^6, the leading constant is 6% larger than optimal.
    
    The downside of picking a larger base factorization is that we pay a larger additive constant term. In particular, we need to round m up to the next power of m_1, which could be excessive if m_1 is large.
    
    

 §.§ Conclusion

    
    Our iterative construction provides an algorithm that is similar in computational efficiency to the binary tree mechanism, but which has better constants.
    In principle, the leading constant can be made arbitrarily close to optimal.
    Table <ref> shows how the leading constant improves as we increase the size of the base factorization.
    However, increasing the size of the base factorization makes the terms other than the leading constant larger and the computational efficiency deteriorates too.
    
    In particular, as Figure <ref> shows, the objective value grows like a staircase, because we must round up the size to the next power of the base factorization size.
    These steps are a significant low order term.
    
    To avoid these limitations we can combine multiple base factorizations, rather than a single one. This will “smoothen” the staircase. E.g., rather than rounding up the size to 2^a or 3^b, we can round up to a size of the form 2^a · 3^b.
    
  


§ NUMERICAL LOWER BOUND ON OPTIMAL PERFORMANCE

In this section, we develop numerical lower bounds on the optimal matrix factorizations of the all-ones lower triangular matrix A^(n) for various classes of matrices. More precisely, given the sequence length n, we write down a semidefinite program that provides a lower bound on the lowest achievable error over a specific class of correlation matrices (for example, Toeplitz matrices or triangular matrices). 

In general, we consider mechanisms whose matrices can be written as 
C=c=∑_i _i c_i where c ∈^m denotes the parameters of the mechanism (to be optimized) and : ^m ↦^n × n is an linear operator and _i ∈^n × n are “basis” matrices. Any linear subspace of the space of n × n matrices (including lower triangular matrices and Toeplitz matrices) can be expressed in this form.

All results in this section derive from the following fundamental result which is a straightforward consequence of weak duality:


Let  denote the n × n lower triangular matrix with all entries equal to 1. 
For any c ∈^m, define the objective for the correlated mechanism with correlation matrix c by

    Fc=(c,c) = c_1 → 2^2c_2 →∞^2.

Further, for any set ℭ⊆^m, we have 

    min_c ∈ℭ Fc≥1/min_Γ_0, …, Γ_N ∈^n_+
    ∑_i _i^T Γ_i _i=1max_c_1↦ 2=1cc∑_i Γ_i

where ^n_+ denotes the space of symmetric positive semidefinite n × n matrices, _i denotes the i-th column of  and AB=traceAB denotes the standard inner product between symmetric matrices A, B ∈^N.

Furthermore, if the inner maximization can be relaxed to a concave maximization problem, i.e, there exists a compact convex set 𝒞⊂^N such that

    max_c ∈ℭ
    c_1↦ 2=1ccΓ = max_S ∈𝒞SΓ  ∀Γ∈^N_+

then the inequality in <ref> is an equality. 


Since the objective is invariant to scaling c by any positive constant, i.e., Fα c=Fc for all α∈, α > 0, we can set α=1/c_1 → 2 so that the optimization problem becomes

    min_c ∈ℭ
    c_1 → 2=1c_2→inf^2=min_c ∈ℭ
    c_1 → 2=1max_ic_i^2


From Lemma <ref>, we have that

    M = 1/min_Γ∈^N_+ 
    ΓΓM

for any M ∈^n_++. Applying this to M=cc and =_i (for i=0, …, n) and plugging this into (<ref>), we obtain

    max_ic_i^2 = max_i max_Γ_i ≽ 0 
    _i^T Γ_i _i=11/ccΓ_i=1/min_i min_Γ_i ≽ 0 
    _i^T Γ_i _i=1ccΓ_i

The denominator can be rewritten as

    min_ν∈Δ^n∑_i ν_imin_Γ_i ≽ 0 
    _i^T Γ_i _i=1ccΓ_i=min_ν∈Δ^n∑_i min_Γ_i ≽ 0 
    _i^T Γ_i _i=ν_iccΓ_i
       =min_Γ_0, …, Γ_N ∈^n_+
    ∑_i _i^T Γ_i _i=1cc∑_i Γ_i

where Δ^n is the simplex in n dimensions. Thus, the optimization problem (<ref>) can be rewritten as

    min_c_1 → 2≤ 11/min_Γ_0, …, Γ_N ∈^n_+
    ∑_i _i^T Γ_i _i=1cc∑_i Γ_i = 1/max_c_1 → 2≤ 1min_Γ_0, …, Γ_n ∈^n_+
    ∑_i _i^T Γ_i _i=1cc∑_i Γ_i

By weak duality, we have that the optimal value is bounded below by

    1/min_Γ_0, …, Γ_n ∈^n_+
    ∑_i _i^T Γ_i _i=1max_c_1 → 2≤ 1cc∑_i Γ_i

which gives <ref>.

Furthermore, if the inner maximization can be rewritten as a concave maximization problem, the overall problem is a convex-concave problem and strong duality holds by the Von-Neumann minimax theorem <cit.>, so that the order of the min and max can be interchanged without changing the optimal value. Hence, the lower bound is tight.




Let 𝒥 = {c: c ∈ℭ}. A tight lower bound in <ref> is achieved in the following 3 cases:

    
  * 𝒥 is the set of all lower triangular matrices.
    
  * 𝒥 is the set of all lower triangular Toeplitz matrices.
    
  * 𝒥 is the set of all lower triangular Toeplitz matrices formed from a degree 1 constant recurrent sequence.



1. If 𝒥 is the set of all lower triangular matrices, for any 𝕁∈𝒥, letting Δ=𝕁𝕁, we have that 

    max_𝕁_1 → 2≤ 1, 𝕁∈𝒥𝕁𝕁Γ≤max_Δ≽ 0 
    diagΔ≤ 1ΔΓ

since for any feasible 𝕁, Δ=𝕁𝕁 satisfies the constraints listed in the optimization problem on the right hand side (note that diagΔ refers to the vector of diagonal elements of the matrix Δ). 

Conversely, given any feasible solution Δ satisfying the constraints, one can compute its upper triangular Cholesky factorization (the regular Cholesky factorization multiplied by the matrix with anti-diagonals equal to 1), to obtain 𝕁 that satisfies the constraint 𝕁_1 → 2≤ 1 and Δ=𝕁𝕁. Hence the two optimization problems have equal optimal values, and hence the lower bound in theorem <ref> is tight.

2. If 𝒥 is the set of all n × n lower triangular Toeplitz matrices, let 𝕁c denote the lower triangular Toeplitz matrix whose first column is c. Then, we have 

    c_1→ 2 = c, 𝕁c=∑_i c_i H_i

where H_i is the matrix with 0s everywhere except the i-th subdiagonal equal to 1. Thus

    max_𝕁c_1 → 2≤ 1, 𝕁c∈𝒥𝕁𝕁Γ = max_c = 1𝕁c𝕁cΓ

and the right hand side is a standard trust region subproblem with a quadratic objective and a single quadratic constraint, that can be solved exactly via an SDP relaxation by the S-lemma <cit.>, i.e, there exists an affine map Q: ^n →^n such that the right hand side equals max_C ≽ 0, traceC=1QΓ C. Hence, the lower bound in <ref> is tight.


3. If 𝒥 is the set of all n × n lower triangular Toeplitz matrices formed from a degree 1 constant recurrent sequence, we have that 𝕁=𝕁c where c_i = αθ^i. where α, θ∈ are arbitrary real numbers. Then, the problem

    max_c = 1𝕁c𝕁cΓ=max_c𝕁c𝕁cΓ/cc

can be written as

    max_θ, αα^2 pθ;Γ/α^2 qθ = max_θpθ;Γ/qθ

where p, q are polynomials in θ and the coefficients of q are independent of Γ while the coefficients of p are affine functions of Γ. Using the fact that 
    pθ-β qθ≤ 0  pθ-β qθ is a sum of squares polynomial

we can write the maximization problem as being equal to

    min_pθ;Γ-β qθ is a sum of squares polynomialβ

Thus, the inner maximization can be written as a concave maximization problem of the form required for lower bound tightness in Theorem <ref> <cit.>.






 §.§ Lower bound over all lower triangular matrices





 §.§ Lower bound over Toeplitz matrices


Let C ∈^N × N denote the Toeplitz lower triangular correlation matrix and c denote its first column. Its clear that the entries of C are affine functions of C and we use C=c=∑_i _i c_i to denote that where : ^N ↦^N × N is an linear operator and 

    _i ∈^N × N
    = [ 0_N-i, N-i   0_N-i, i; 𝕀_N-i, N-i   0_N-i, i ]

𝕀_k is the k × k identity matrix and 0_m, n is the m × n matrix of 0s. The inf-norm error of the prefix sum estimate is given by

    Fc=c^2 c^2

where · denotes the Euclidean ℓ_2 norm in ^N and ∈^N is the vector with all entries equal to 1.




min_c ∈^N Fc is bounded below by the optimal value of the following semidefinite prgram:



    max_t ∈, C ∈^N      t 
    Subject to    [                   t                    ;                     ∑_i, j=1^N_i_j C_ij ]≽ 0 
        C ≽ 0 
       trace[C] = 1




Furthermore, the choice of c=c^⋆ developed in <cit.> which is defined by c^⋆_0 = 0 and c^⋆_i+1 = c^⋆_i1 - 1/2(i+1) is such that Fc^⋆ equals the optimal value of the optimization problem (<ref>)




 §.§ Lower bound over degree 1 constant recurrent sequence matrices


Further, if we impose additional constraints on c, for example, that c is an order d constant recurrence sequence, we can get tighter lower bounds on the optimal performance. In particular, for d=1, we have the following result:

Let ℂ_d denote the set of constant recurrent sequences of order d. Then, min_c ∈ℂ_1 Fc is bounded below by the optimal value of the semidefinite program


    min_t ∈, C ∈^N      t 
    Subject to    [                     t                      ;                       ∑_i, j=0^N-1_i_j C_ij ]≽ 0 
        C ≽ 0 
       trace[C] = 1 
        C_i,k-i = C_0, k  ∀ i ∈{0, …, k}, k ∈{1, …, N-1}






[Lower bound on matrix fractional function]
For any ∈^n and M ∈^n_++, we have

    M = 1/min_Γ∈^n_+ 
    Γ = 1ΓM



Let ∈^n be an arbitrary vector. 
Now, consider any positive definite matrix M ∈^n_+ and eigendecomposition M=∑_i λ_i u_iu_i for some orthogonal basis {u_i}. Further, expressed in this basis, suppose that =∑_i σ_i u_i. 

Then, we have

    M^-1 = ∑_i 1/λ_iu_i^2 = ∑_i σ_i^2/λ_i

Let MΓ denote the inner product in ^n, i.e., MΓ=traceMΓ.
Further, for any matrix Γ∈^n_+, we have

    M^-1MΓ = ∑_i σ_i^2/λ_i∑_i λ_i u_iΓ u_i =∑_i σ_i^2/λ_i∑_i λ_i Γ^1/2u_i^2
       ≥∑_i σ_iΓ^1/2u_i^2 ≥∑_i σ_i Γ^1/2 u_i^2 = ^T Γ

where the first inequality is an application of the Cauchy-Schwartz inequality and the second inequality is an application of the triangle inequality. 

Thus, for any M, Γ∈^n_+, we have

    M^-1≥^T Γ/MΓ

Further, choosing Γ=MM, the inequality is replaced with equality. Hence, we have 

    M^-1 = max_Γ≽ 0^T Γ/MΓ = max_Γ≽ 0 
    ^T Γ=11/MΓ










§ WITH D=1 BUFFERS





In this section, we show that it is possible to achieve (B,C) = O(n^1/6) error for a degree d=1 constant recurrent sequence, without appealing to the sophisticated machinery of rational function approximation that build the crux of the paper. Unfortunately, this approach does not naturally extend to degrees d>1. In comparison, adding independent noise (that is when d=0) results in an ℓ_∞ error of O(n).

We will use the usual parameterization from rest of the paper:

    ℳ(x) = ^(n) + Bz=( +^-1·)

In <Ref>,  is the lower triangular matrix of all ones,  is factorized into =·, where  and  are lower-triangular Toeplitz matrices. = =^-1 is a vector of (potentially) non-istropic Gaussian noise, and ∼(0,σ^2·𝕀_n× n). The data vector  has _∞≤ 1.
We consider the factorization A=BC, where B=AC^-1 and   and ^-1 are parameterized as follows. 
Let a,λ∈ [-1,1] with |λ-a^2|≤ 1.
Define

    = ([      1      0      0      0      ⋯      0;    a^2      1      0      0      ⋯      0;   a^2λ    a^2      1      0      ⋯      0; a^2λ^2   a^2λ    a^2      1      ⋯      0;      ⋮      ⋮      ⋮      ⋮      ⋮      ⋮ ]),   ^-1=([             1             0             0             0             ⋯             0;          -a^2             1             0             0             ⋯             0;   -a^2(λ-a^2)          -a^2             1             0             ⋯             0; -a^2(λ-a^2)^2   -a^2(λ-a^2)          -a^2             1             ⋯             0;             ⋮             ⋮             ⋮             ⋮             ⋮             ⋮ ]).

That is, C corresponds to the generating function c(x) = 1 + a^2x/1-λ x, which clearly has degree 1.
It can be checked that the expression for C^-1 is correct (via <ref>) as 1/c(x) = 1 - a^2x/1-(λ-a^2)x.

 
Setting λ=1-1/n^2/3, and a^2=1/n^1/3(1-1/n^1/3) in <ref> gives 

    C_1→2^2    ≤ O(1),
    B_2→∞=A C^-1_2→∞   ≤ O(n^1/6)





The maximum squared column norm of  is

    C_1→2^2 = 1+ a^4(1+λ^2+λ^4+⋯λ^2(n-2))≤ 1+a^4/1-λ^2.

Since 1-λ^2 = 2 n^-2/3 - n^-4/3 = Θ(n^-2/3) and a^4 = n^-2/3(1-n^-1/3)^2 = Θ(n^-2/3), we have C_1→2=Θ(1).


Let =[b_n-1,⋯,b_0] be the last row of =^-1. This is the longest row, so B_2→∞^2 = b_2^2. 
In the following, we bound each element of the vector . We have b_0=1 and, for any i∈{1,⋯ n-1},

    b_i= 1-a^2(1+(λ-a^2)+⋯+(λ-a^2)^i-2)=1-a^2(1-(λ-a^2)^i-1)/1+a^2-λ.


Thus

    B_2→∞^2    =1 + ∑_i=1^n-1(1-a^2(1-(λ-a^2)^i-1/1+a^2-λ)^2 
       = 1 + ∑_i=1^n-1((1-λ)+a^2(λ-a^2)^i-1/1+a^2-λ)^2 
       = 1 + 1/(1+a^2-λ)^2∑_i=1^n-1 (1-λ)^2 + 2(1-λ)a^2(λ-a^2)^i-1 + a^4(λ-a^2)^2(i-1)
       = 1 + (n-1)(1-λ)^2 + 2(1-λ)a^21-(λ-a^2)^n-1/1+a^2-λ+a^41-(λ-a^2)^2(n-1)/1-(λ-a^2)^2/(1+a^2-λ)^2
       ≤ 1 + (n-1)(1-λ)^2 + 2(1-λ)a^21/1+a^2-λ+a^41/1-(λ-a^2)^2/(1+a^2-λ)^2.

Now a^2 = n^-1/3-n^-2/3, 1-λ = -n^-2/3, 1-λ+a^2 = n^-1/3 and 1-(λ-a^2)^2 = 2n^-1/3-n^-2/3=Θ(n^-1/3).
Substituting these into <ref> gives

    B_2→∞^2    ≤ 1 + (n-1)n^-4/3 - 2 n^-2/3 (n^-1/3-n^-2/3)/n^-1/3 + (n^-1/3-n^-2/3)^2/(2n^-1/3-n^-2/3)/n^-2/3
       = 1 + n-1/n^2/3 - 2 + 2/n^1/3 + n^-2/3-2n^-1+n^-4/3/2n^-1 - n^-4/3
       = 1 + n-1/n^2/3 - 2 + 2/n^1/3 + n^1/3-2+n^-1/3/2 - n^-1/3
       = Θ(n^1/3)

as required.
For brevity, let 0<x=1-λ<1. We have, the following:

    2ρ·()≤(∑_i=2^n(x+a^2(1-(x+a^2))^i-1/x+a^2)^2_E+1)·(1+a^4/1-(1-x)^2)_F


Let y=(x+a^2) for brevity. We will assume that 0<1-y≤ 1, and we will enforce this choice when we make concrete instantiations of the parameters in terms of n.

    E   = ∑_i=2^n(x+a^2(1-y)^i-1/y)^2=∑_i=2^n(x/y+a^2/y·(1-y)^i-1)^2
       < 2(nx^2/y^2+a^4/y^2∑_n=1^n-1(1-y)^2i)≤ 2(nx^2/y^2+a^4/y^2·(1-y)^2/1-(1-y)^2)
       =2(nx^2/y^2+a^4/y^2·1/y(2-y))≤ 2(nx^2/y^2+a^4/y^3)≤ 2(nx^2/y^2+a^4/y^3).

We now bound F in <Ref>.

    F   =(1+a^4/1-(1-x)^2)=(1+a^4/x(2-x))≤(1+a^4/x).

Plugging the values of E and F respectively from <Ref> and <Ref> in <Ref>, we have 

    2ρ·()   ≤ 2(nx^2/y^2+a^4/y^3)(1+a^4/x)=2(nx^2/y^2+a^4/y^3+nx a^4/y^2+a^8/xy^3).

Let α=x/a^2+x. We have from <Ref> the following.

    2ρ·()   =O(nα^2+(1-α)^2/y+nα(1-α)^2y+(1-α)^4/α)
       ≤ O(nα^2+1/y+nα y+1/α).

Setting nα^2=1/y (to balance the first two terms) and α=1/n^1/3 (to optimize for α), we get the following:

    2ρ·()   = O(nα^2+1/α)=O(n^1/3).

The corresponding choice of the parameters in <Ref> are λ=1-1/n^2/3, and a^2=1/n^1/3(1-1/n^1/3). This completes the proof.







§ CONJECTURE


We give the following conjecture which improves the results in <ref>. If true, this would imply that can work with d = O(log n) space complexity. 


    For all n ∈ℕ there exists a rational function r : ℂ→ℂ of degree d=O(log n) such that the following hold.
    
    1/2π∫_-π^π| 1/r(x) - 1/√(1-x)|^2 dθ≤1/n    where     x = exp(√(-1)θ - 1/n).

    
    1/2π∫_-π^π| r(x)/1-x - 1/√(1-x)|^2 dθ≤1/n    where     x = exp(√(-1)θ - 1/n).


Compared to what we prove in <ref>, this conjecture is quantitatively stronger in that d=O(log n) instead of d = O(log^2 n) but qualitatively weaker in that (i) we prove an approximation guarantee for |x|≤exp(-1/n) rather than |x|≤1 and, (ii) rather than a uniform bound on |r(x)-√(1-x)|, we bound the integrals directly.



§ ACKNOWLEDGEMENTS

tocsectionAcknowledgements
We thank Andreas Terzis, Shuang Song, and Arun Ganesh for comments on the draft.
We thank Aleksandar Nikolov and Jalaj Upadhyay for helpful discussions about prior work.
We thank Mark Bun for making us aware of the results of <cit.>.


tocsectionReferences






