00footnotetext: Alphabetical author order.

# Efficient and Near-Optimal Noise Generation
for Streaming Differential Privacy

Krishnamurthy (Dj) Dvijotham Google DeepMind<EMAIL_ADDRESS>H. Brendan
McMahan Google Research<EMAIL_ADDRESS>Krishna Pillutla IIT Madras. Work
done at Google<EMAIL_ADDRESS>Thomas Steinke Google DeepMind
<EMAIL_ADDRESS>Abhradeep Thakurta Google DeepMind<EMAIL_ADDRESS>
###### Abstract

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.

###### Contents

  1. 1 Introduction
    1. 1.1 Our Contributions
    2. 1.2 Our Algorithms & Techniques
    3. 1.3 An Empirical Comparison of Mechanisms
  2. 2 Background
    1. 2.1 Prior Work on Matrix Factorizations for Continual Counting
    2. 2.2 Lower Triangular Toeplitz Factorizations versus General Factorizations
    3. 2.3 Applications of Matrix Factorizations in Machine Learning
    4. 2.4 Other Related Work
  3. 3 Efficiently-Sampleable Factorization via Rational Generating Functions
    1. 3.1 Sequences, Lower Triangular Toeplitz Matrices, & Generating Functions
    2. 3.2 Rational Generating Functions and Constant-Recurrent Sequences
    3. 3.3 Efficient Sampling via BLT Multiplication
  4. 4 Factorizations via Rational Function Approximation
    1. 4.1 Reduction to Approximating the Square Root
    2. 4.2 Rational Approximation of $\sqrt{x}$
    3. 4.3 Putting Everything Together
  5. 5 Factorizations via Direct Optimization
    1. 5.1 Fast and Exact Error Computation
    2. 5.2 Direct Optimization of BLTs
    3. 5.3 Derivation of the Inverse BLT Parameterization
  6. 6 Generalizing the Binary Tree Approach
    1. 6.1 Recursive Matrix Factorization
    2. 6.2 Recursive Algorithm
    3. 6.3 Combining Recursion with BLTs
  7. 7 Numerical Lower Bound on Optimal Performance
  8. 8 BLT with $d=1$ Buffers
  9. 9 Conjecture
  10. Acknowledgements

##  1 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 [DMNS06].

_Continual counting_ [DNPR10, CSS11] 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 [BNS13, BNSV15, BDRS18, ALMM19,
KLMNS20, GJK21, CLNSS23].

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 [KMSTTX21, ZTC22, KMSTTX21, DMRST22, CCMRT22,
CCGMMRGTX24, XZACCKMRZ23]. (See Section 2.3 for discussion of this
application.) This has motivated significant work on making differentially
private continual counting as accurate as possible [Mat93, MNT20, FHU22,
HUU23, HUU24, AP24].

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 to achieve good utility. 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 $n$ scalars;
this discussion readily extends to summing vectors. Let $x\in\mathbb{R}^{n}$
contain the $n$ terms we wish to sum (over $n$ iterations or steps). Let
$A\in\\{0,1\\}^{n\times n}$ be the lower-triangular all-ones matrix so that
$Ax$ is the vector of partial sums – i.e., $(Ax)_{k}=\sum_{i\leq k}x_{i}$ for
each $k$. We factor $A=BC$, where $B,C^{T}\in\mathbb{R}^{n\times n^{\prime}}$.
The corresponding differentially private continual counting algorithm
$\mathcal{M}$ is given by

$\mathcal{M}(x)\coloneqq B(Cx+z)=Ax+Bz=A(x+C^{\dagger}z),$ (1.1)

where $z\leftarrow\mathcal{N}(0,\sigma^{2}I)$ is a vector of independent
Gaussian noise.111$C^{\dagger}=A^{-1}B$ 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^{\dagger}$ is a suitable pseudo-inverse. Equivalently,
the output of the mechanism can be written as $\mathcal{M}(x)=Ax+\widehat{z}$
where $\widehat{z}\leftarrow\mathcal{N}(0,\Sigma)$ with
$\Sigma=\sigma^{2}BB^{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\to 2}\coloneqq\max_{j}\sqrt{\sum_{i}C_{i,j}^{2}}.$

This is the $L_{2}$ sensitivity of $Cx$. Then, following the standard privacy
analysis of Gaussian noise addition [DR14, Ste22], to ensure that $Cx+z$ is
$(\varepsilon,\delta)$-differentially private, we scale the noise
$z\leftarrow\mathcal{N}(0,\sigma^{2}I)$ to have standard deviation
$\sigma=\zeta\|C\|_{1\to 2}$, where
$\zeta=O(\frac{1}{\varepsilon}\sqrt{\log(1/\delta)})$ is a noise multiplier
depending only on the privacy parameters
$(\varepsilon,\delta)$.222Importantly, 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,\dots,k-1$
already released by the mechanism [DMRST22]. By postprocessing, this implies
$\mathcal{M}(x)=B(Cx+z)$ is $(\varepsilon,\delta)$-differentially private.

#### Mechanism utility:

The error of the mechanism $\mathcal{M}$ is given by $\mathcal{M}(x)-Ax=Bz$.
Thus, the root mean squared error of the $i$th partial sum can be calculated
as

$\displaystyle\sqrt{\mathbb{E}_{\mathcal{M}}{(\mathcal{M}(x)-Ax)_{i}^{2}}}=\sqrt{\mathbb{E}_{z}\left[(Bz)_{i}^{2}\right]}=\sqrt{\sigma^{2}\sum_{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 maximum333We can also consider other measures of error
such as sum of variances (rather than the max). This corresponds to a
different norm:
${\mathbb{E}\left[\|Bz\|_{2}^{2}\right]}=\sigma^{2}\,\left\|B\right\|_{F}^{2}$.
such error over all $n$ partial sums released, which scales with the maximum
row norm

$\|B\|_{2\to\infty}\coloneqq\max_{i}\sqrt{\sum_{j}B_{i,j}^{2}}.$

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

$\displaystyle\max_{i}\sqrt{{\underset{\mathcal{M}}{\mathbb{E}}\left[(\mathcal{M}(x)-Ax)_{i}^{2}\right]}}$
$\displaystyle=\max_{i}\sqrt{{\underset{z\leftarrow\mathcal{N}(0,\sigma^{2}I)}{\mathbb{E}}\left[(Bz)_{i}^{2}\right]}}=\max_{i}\sqrt{\sigma^{2}\sum_{j}B_{i,j}^{2}}$
$\displaystyle=\sigma\|B\|_{2\to\infty}=\zeta\|B\|_{2\to\infty}\|C\|_{1\to
2}=:\zeta\operatorname{MaxErr}(B,C),$ (1.2)

where we take $\sigma=\zeta\|C\|_{1\to 2}$ to ensure
$(\varepsilon,\delta)$-differential privacy, and we define

$\operatorname{MaxErr}(B,C)\coloneqq\|B\|_{2\to\infty}\|C\|_{1\to
2}=\sqrt{\max_{i}\sum_{j}B_{i,j}^{2}}\sqrt{\max_{j}\sum_{i}C_{i,j}^{2}}.$
(1.3)

Since the error in Section 1 can be written as a product of the noise
multiplier $\zeta$ (which does not depend on the factorization $B,C$) and
$\operatorname{MaxErr}(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 $\zeta$ for the
remainder of this paper.

#### Matrix factorization:

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

$\text{minimize}_{B,C}~{}\operatorname{MaxErr}(B,C)~{}~{}~{}~{}~{}\text{
subject to }BC=A.$ (1.4)

The optimal value of this objective is denoted $\gamma_{2}(A)$ and is known as
the _gamma-two factorization norm_ of $A$. Prior work has obtained near-
optimal factorizations of $A$ (and related matrices) [Mat93, MNT20]. The
optimal value is

$\gamma_{2}(A)\coloneqq\inf\big{\\{}\operatorname{MaxErr}(B,C):B,C\in\mathbb{R}^{n\times
n},BC=A\big{\\}}=\frac{\log(n)}{\pi}\pm O(1).$ (1.5)

However, the existing factorizations are either far from optimal in terms of
their utility or are impractical due to their high computational cost, which
we turn to next.

#### Computationally efficient noise generation:

The computational cost of implementing this matrix factorization approach for
differentially private continual counting is dominated by the matrix-vector
multiplication $Bz$, where $z\in\mathbb{R}^{n}$ consists of i.i.d. Gaussian
noise $z\leftarrow\mathcal{N}(0,\sigma^{2}I)$. The main computational
challenge is that our algorithm cannot store the increments $x$ or the seed
noise $z$ (or even the matrices $B$ or $C$) in memory; this would require
memory that is linear in the number of steps $n$, which can be prohibitively
large.

To give a sense of scale, in private machine learning applications, we often
need to compute $n>10^{6}$ partial sums (each coming from an iteration of
stochastic gradient descent or a related algorithm), over $m>10^{10}$
dimensions (the size of the model and its gradients).444In this introductory
discussion, we have considered sums over $m=1$ dimension. Extending to $m>1$
is essentially a matter of running $m$ 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 $\mathcal{M}(x)$
without storing all of the data or noise in memory, as this would take
$nm>10^{16}$ units of memory – i.e., over 40 petabytes if we use 32-bit
floating point numbers.555Throughout 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 [Mir12], but techniques exist to
ensure differential privacy in discrete settings [CKS20, KLS21].

In this work, we address the computationally efficiency of sampling the
correlated noise. Specifically, _we provide a matrix factorization mechanism
for differentially private continual counting attaining near-optimal error
with only polylogarithmic memory overhead._

#### Streaming setting:

Our algorithm receives the coordinates of the input $x$ and i.i.d. noise $z$
one at a time and outputs the approximate partial sums $\mathcal{M}(x)$ one at
a time. In particular, we must output each partial sum
$\left({\mathcal{M}(x)}\right)_{k}\approx\sum_{i\leq k}x_{i}$ before receiving
the next input $x_{k+1}$.

Ideally, the memory usage should not grow with the stream length $n$ 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
$\operatorname{MaxErr}(B,C)=\sqrt{n}$ 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 $\operatorname{MaxErr}(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, $\forall 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 (1.4). Furthermore, Toeplitz matrices
(and certain specializations we consider) have convenient mathematical and
algorithmic properties. See Figure 2 and Section 2.2 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.

###  1.1 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.

######  Theorem 1.1 (Main Result – Informal version of Theorem 4.6).

For each integer $n\geq 1$ and error parameter $\mu\in(0,1)$ , there exists a
lower triangular Toeplitz matrix factorization $B,C\in\mathbb{R}^{n\times n}$
with the following properties.

  * •

Validity: $BC=A$, where $A$ is the $n\times n$ lower-triangular all-ones
matrix.

  * •

Near-optimality: $\operatorname{MaxErr}(B,C)\leq\mathsf{OptLTToe}(n)+\mu,$
where

$\mathsf{OptLTToe}(n)=1+\sum_{k=1}^{n-1}\left(2^{-2k}{2k\choose
k}\right)^{2}\leq 1+\frac{0.57722+\log(n)}{\pi}$

is the optimal value of $\operatorname{MaxErr}(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^{-1}z)_{k}$) and
runs in time and space $O(\log^{2}(n/\mu))$.

In order to ensure that the error term $\mu$ in our result is $o(1)$, we have
space and time complexity $d=\Theta(\log^{2}n)$. It is natural to wonder
whether this computational complexity can be improved. We show that it can be
improved to $\widetilde{O}(\log n)$ at the expense of a weaker multiplicative
near-optimality guarantee for the matrix factorization objective.666 In this
problem, the multiplicative constants in the error have a larger practical
impact on the mechanism’s utility than an additive error. For instance, the
binary tree mechanism is suboptimal by a multiplicative factor of $\pi/\log
2\approx 4.5$, which yields a much worse $\operatorname{MaxErr}$ than the
factorization of [FHU22] (cf. Figure 2, left).

######  Theorem 1.2 (Secondary Result – Informal version of Theorem 6.6).

For each integer $n\geq 1$, there exists an integer $n^{\prime}=O(n)$ and a
matrix factorization $B,C^{T}\in\mathbb{R}^{n\times n^{\prime}}$ with the
following properties.

  * •

Validity: $BC=A$, where $A$ is the $n\times n$ lower-triangular all-ones
matrix.

  * •

Near-optimality: $\operatorname{MaxErr}(B,C)\leq(1+o(1))\cdot\mathsf{Opt}(n),$
where $\mathsf{Opt}(n)=\frac{\log(n)}{\pi}\pm O(1)$ is the optimal value of
$\operatorname{MaxErr}(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) $\widetilde{O}(\log n)$.

The optimal matrix factorization objective value over lower-triangular
Toeplitz factorizations $\mathsf{OptLTToe}(n)$ is a small additive constant
(specifically, $0.365$) away from the optimal over all factorizations
$\mathsf{Opt}(n)$ (cf. Corollary 2.4). Thus, Theorem 1.1’s approximation bound
on the class of lower triangular Toeplitz matrices can be directly compared to
Theorem 1.2’s approximation bound over all possible factorizations.

More generally, we can smoothly trade off between the matrix factorization
objective and computational efficiency. That is, we can interpolate between
Theorems 1.1 and 1.2; see Proposition 6.5 for a general statement.

Theorem 1.2 does not generate a lower triangular Toeplitz factorization; in
fact it does not even produce a square factorization. It also gives a weaker
multiplicative near-optimality guarantee. We leave it as an interesting open
problem whether it is possible to improve on $\log^{2}n$ space complexity with
a Toeplitz factorization or with
$\operatorname{MaxErr}(B,C)\leq\mathsf{Opt}(n)+O(1)$.

Figure 1:  Ratio of $\operatorname{MaxErr}(B,C)$ of our RA-BLT and Opt-BLT
mechanisms for different numbers of steps $n$ and degrees $d$ (which
corresponds directly to the space complexity) over that of the optimal
Toeplitz mechanism of [FHU22]. This illustrates that even with modest degree
$d$, we obtain very good $\operatorname{MaxErr}(B,C)$ even for large numbers
of steps $n$. For example, Opt-BLT with $d=5$ is within $1\%$ of optimal for
$n=10^{7}$ (we do not plot Opt-BLT for $d=9$).

#### Practical algorithms:

Since our work is directly motivated by practical considerations, we also
study the problem numerically. Specifically, we show that the method behind
Theorem 1.1 yields a factorization that is indistinguishable from optimal
Toeplitz mechanism for practical purposes, but which has a highly efficient
procedure for generating the noise. Figure 1 shows the ratio of
$\operatorname{MaxErr}(B,C)$ for our algorithm over the optimal
$\mathsf{OptLTToe}(n)$.

To prove Theorem 1.1 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
$\operatorname{MaxErr}(B,C)$ in terms of this parameterization, and (iii) show
that this expression is differentiable. This allows us to numerically optimize
our factorization using gradient-based methods. While we do not prove that
this optimization procedure converges, in practice it yields significantly
better solutions than the closed-form parameters used to prove Theorem 1.1.

To further illustrate the power of our approach, we show in Section 8 that
even with constant memory, we can asymptotically attain
$\operatorname{MaxErr}(B,C)=O(n^{1/6})$. This is already an improvement over
$\operatorname{MaxErr}(I,A)=\sqrt{n}$ obtained by the trivial factorization.

#### Lower bounds:

In Section 2.2 we prove a lower bound that exactly characterizes
$\mathsf{OptLTToe}(n)$. Specifically, we show that the lower triangular
Toeplitz factorization of [FHU22] is precisely optimal for this class.

Further, in Section 7, 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 correspond
to our algorithm with a specific constant memory constraint.

###  1.2 Our Algorithms & Techniques

The starting point for our main result is the lower triangular Toeplitz
factorization of [FHU22], which we show is in fact the optimal lower
triangular Toeplitz factorization in Section 2.2. This factorization is given
by $B_{i,j}=C_{i,j}=f_{i-j}$ where the sequence $f_{0},f_{1},\cdots$ is the
coefficients of a Taylor series777We overload notation and use $x$ both for
the input to the algorithm ($x\in\mathbb{R}^{n}$), and the indeterminate in
polynomials and generating functions; the meaning should be clear from
context.:

$\frac{1}{\sqrt{1-x}}=f_{0}+f_{1}x+f_{2}x^{2}+f_{3}x^{3}+\cdots.$ (1.6)

This sequence satisfies the recurrence $f_{k}=(1-1/2k)\cdot f_{k-1}$ for
$k\geq 1$ with $f_{0}=1$. Our task in this case is to compute the correlated
noise $(Bz)_{k}=\sum_{j=0}^{k}B_{k,j}\cdot z_{j}=\sum_{j=0}^{k}f_{k-j}\cdot
z_{j}$ in a streaming fashion for $k\in[n]\coloneqq\\{0,1,\cdots,n-1\\}$.888We
zero-index sequences, vectors, and matrices throughout. See Table 1 for a
summary of symbols and notation. Unfortunately, the factorization of Eq. 1.6
does not seem to admit an efficient sampling or noise generation procedure.

####  BLTs: Buffered Linear Toeplitz Matrices (§3):

Suppose the sequence instead satisfied a linear recurrence $f_{k}=q\cdot
f_{k-1}$ for all $k\geq 1$. Then

$(Bz)_{k}=\sum_{j=0}^{k}f_{k-j}\cdot z_{j}=f_{0}\cdot
z_{k}+\sum_{j=0}^{k-1}(q\cdot f_{k-j-1})\cdot z_{j}=f_{0}\cdot
z_{k}+q\cdot(Bz)_{k-1}.$ (1.7)

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}\cdot f_{k-1}+q_{2}\cdot f_{k-2}+\cdots+q_{d}\cdot f_{k-d}.$ (1.8)

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

The recurrence $f_{k}=qf_{k-1}$ implies the closed form $f_{k}=q^{k}f_{0}$.
Similarly, the recurrence in Equation 1.8 implies a closed form expression
$f_{k}=u^{T}W^{k}v$ where $W\in\mathbb{R}^{d\times d}$ is a matrix and
$u,v\in\mathbb{R}^{d}$ are vectors. This closed form is what we use in Section
3.3 for our Algorithm 1. Specifically, we can extend Equation 1.7 to this
matrix-power closed form

$\displaystyle(Bz)_{k}$ $\displaystyle=\sum_{j=0}^{k}f_{k-j}\cdot
z_{j}=\sum_{j=0}^{k}u^{T}W^{k-j}v\cdot z_{j}=u^{T}S_{k+1}$

for a suitable state vector $S_{k}\in\mathbb{R}^{d}$, stored in memory.
Namely, we initialize $S_{0}=0$ and, at each iteration, our algorithm updates

$S_{k+1}=v\cdot z_{k}+WS_{k}$

and then outputs $(Bz)_{k}=u^{T}S_{k+1}$. We refer to the entries of $S_{k}$
as the $d$ _buffers_ of our algorithm.999Recall 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
$m>10^{9}$, so keeping $d$ 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 ($vz_{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” (BLTs). We
overload the acronym and write BLTs as a shorthand encompassing matrices,
factorizations, and mechanisms.

Unfortunately, the optimal factorization [FHU22] does _not_ satisfy a
recurrence like Eq. 1.8 and cannot be expressed as a BLT. Hence, our approach
is to _approximate_ the optimal factorization using BLTs.

#### Designing BLTs via rational function approximation (§4):

If we view the Toeplitz sequence $f_{0},f_{1},\dots$ as being defined by an
ordinary generating function101010We 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 [KP11]. $f$ as in Equation 1.6, it turns
out that satisfying a linear recurrence as in Equation 1.8 is equivalent to
the function being rational with degree at most $d$, i.e., $f(x)=p(x)/q(x)$
for polynomials $p$ and $q$ of degree $\leq d$. (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 BLTs: we need a
low-degree rational _approximation_ to the function $f(x)={1}/{\sqrt{1-x}}$
from Eq. 1.6 that underlies the optimal factorization of [FHU22].

We appeal to known results in approximation theory. Specifically, it is known
that the function $x\mapsto\sqrt{1-x}$ can be uniformly approximated on the
unit complex disc $\\{x\in\mathbb{C}:|x|\leq 1\\}$ with error $\eta>0$ by a
rational function of degree $d=O(\log^{2}(1/\eta))$ [New64, GT19]. That is,
there exists a rational function $r$ of degree $\leq d$, such that
$|r(x)-\sqrt{1-x}|\leq\eta$ for all $x\in\mathbb{C}$ with $|x|\leq 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)=\sum_{k=0}^{\infty}f_{k}x^{k}$ and
$\widetilde{f}(x)=\sum_{k=0}^{\infty}\widetilde{f}_{k}x^{k}$. Then

$\sum_{k=0}^{\infty}|f_{k}-\widetilde{f}_{k}|^{2}=\frac{1}{2\pi}\int_{-\pi}^{\pi}\left|f(x(\theta))-\widetilde{f}(x(\theta))\right|^{2}\,\mathrm{d}\theta\,~{}\text{
where }~{}x(\theta)=\exp(\sqrt{-1}\theta).$ (1.9)

In our case, $f(x)={1}/{\sqrt{1-x}}$ corresponds to the optimal lower
triangular Toeplitz factorization while $\widetilde{f}(x)$ is our rational
approximation – either $\widetilde{f}(x)=1/r(x)$ or
$\widetilde{f}(x)=r(x)/(1-x)$, where $r(x)\approx\sqrt{1-x}$. We are
interested in a finite sum $\sum_{k=0}^{n-1}|f_{k}-\widetilde{f}_{k}|^{2}$,
rather than the infinite sum in Equation 1.9 (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
$\sum_{k=0}^{n-1}|f_{k}-\widetilde{f}_{k}|^{2}$, the near optimality guarantee
of Theorem 1.1 follows from the triangle inequality.

The time and space requirement of Theorem 1.1 is $O(\log^{2}n)$. This
dependence arises from the degree of the rational approximation. A degree at
least $\Omega(\log^{2}(1/\eta))$ is necessary for approximating $\sqrt{1-x}$
with error $\leq\eta$ even for real values $x\in[-1,1]$ [New64]. 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 Theorem 1.1 gives an explicit rational function approximation,
which we can directly convert into a matrix factorization and feed into our
algorithm. We term the BLTs coming from this approach RA-BLTs, with each
choice of the degree $d$ leading to a different (and successively better)
approximation to the optimal Toeplitz factorization.

#### Designing BLTs via direct optimization (§5):

Recall that our goal is to minimize
$\operatorname{MaxErr}(B,C)=\operatorname{MaxErr}(AC^{-1},C)$. In order to
ensure computational efficiency, we restrict the matrix $C$ to the class of
BLTs. While approximation theory lets us directly construct a near-optimal $C$
as outlined above, we can also approach this as an optimization problem and
optimize the BLT parameters that define $C$ numerically.

This optimization is far from straightforward. The class of BLT matrices can
be parameterized in multiple ways. Converting the parameters for $C$ into
$\operatorname{MaxErr}(AC^{-1},C)$ is nontrivial to compute – much less
optimize – when the size $n$ is large. Nevertheless, the class of BLTs 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 BLTs that allows us to efficiently
compute $\operatorname{MaxErr}$ in time practically independent of the size
$n$, specifically $\mathcal{O}(\operatorname{poly}(d)\log n)$, and also to
compute gradients. Being able to compute gradients allows us to optimize
$\operatorname{MaxErr}(AC^{-1},C)$ numerically.

The first challenge is that we need to be able to effectively parameterize
both $C$ and $AC^{-1}$. If $C$ is a BLT, 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 Lemma 3.2 (see also
Eq. 5.2). Of course $1/c(x)=q(x)/p(x)$ is also a rational function, and it in
fact generates the Toeplitz coefficients of $C^{-1}$ (Lemma 3.1), and so
$C^{-1}$ is also a BLT. The lower-triangular matrix of ones, $A$, is trivially
a BLT, and so Lemma 3.1 also implies $B=AC^{-1}$ is a BLT.

In Section 5, Lemma 5.2 shows that given a BLT $C$, we can (explicitly and in
closed-form) derive the BLT parameters of $C^{-1}$, and hence a closed form
for its Toeplitz coefficients.111111More precisely, Lemma 5.2 uses a
parameterization of $p$ and $q$ which leads to closed forms for both $C$ and
$C^{-1}$. Using these closed-form expressions for the Toeplitz coefficients
enables us to directly compute $\operatorname{MaxErr}(B,C)$ for
$B,C\in\mathbb{R}^{n\times n}$ in time $\mathcal{O}(\operatorname{poly}(d)\log
n)$, see Lemmas 5.3 and 5.4.

This is immediately useful, for example in Figure 2 allowing us to plot the
performance of our mechanisms for $n$ up to $10^{8}$ using only a few seconds
of computer time.121212We could have in fact scaled the
$\operatorname{MaxErr}$ calculations for our mechanisms to arbitrary $n$; the
bottleneck is in the computation of the exact $\operatorname{MaxErr}$ for
[FHU22], which requires time $\mathcal{O}(n)$.

More importantly, however, Section 5 shows computing
$\operatorname{MaxErr}(B,C)$ is a differentiable function of the parameters of
the BLTs $B$ and $C$, and hence we can use a gradient-based optimization
method such as L-BFGS to directly minimize $\operatorname{MaxErr}$ targeting a
specific number of steps $n$, where $B$ and $C$ are ${\texttt{BLT}}s$ with $d$
buffers defined by $2d$ parameters. We term these Opt-BLT mechanisms, and they
perform extremely well in practice. For example, for $n=10^{7}$, a
$\text{{Opt-BLT}}(d\\!=\\!4)$ has $\operatorname{MaxErr}$ that is
$1.032\\!\times$ that of the optimal Toeplitz factorization [FHU22], and
$\text{{Opt-BLT}}(d\\!=\\!7)$ is $1.001\\!\times$ optimal; for smaller $n$ the
results are even better, for example for $n=10^{4}$, $\text{{Opt-
BLT}}(d\\!=\\!4)$ achieves $1.001\\!\times$ optimal. Figure 2 gives more
complete results.

#### Generalizations of the binary tree mechanism (§6):

The starting point for Theorem 1.2 is the binary tree mechanism of [DNPR10,
CSS11]. The binary tree mechanism can be viewed as a recursive construction of
a matrix factorization. A recursion of depth $\ell$ yields a matrix
factorization of size $n=2^{\ell}$ and an algorithm running in time and space
$O(\ell)$. 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
$\operatorname{MaxErr}(B,C)$ for this construction is $O(\log n)$ – that is,
it is within a constant factor of optimal. Specifically, the binary tree
mechanism is asymptotically a factor of $\frac{\pi}{\log 2}\approx 4.5$ from
optimal. This factor is significant in practice, as shown in Figure 2.

We combine the binary tree mechanism’s recursive approach with Theorem 1.1 to
get the best of both worlds – near-optimal constants and $\widetilde{O}(\log
n)$ space. This proves Theorem 1.2.

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

$\displaystyle A^{(6)}=\left(\begin{array}[]{cccccc}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\end{array}\right)$
$\displaystyle=\left(\begin{array}[]{cccccc}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\end{array}\right)+\left(\begin{array}[]{cccccc}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\\\
\end{array}\right)$ (1.28) $\displaystyle=\left(\begin{array}[]{ccc}1&0&0\\\
0&1&0\\\ 0&0&1\end{array}\right)\otimes\left(\begin{array}[]{cc}1&0\\\
1&1\end{array}\right)+\left(\begin{array}[]{ccc}0&0&0\\\ 1&0&0\\\
1&1&0\end{array}\right)\otimes\left(\begin{array}[]{cc}1&1\\\
1&1\end{array}\right)$ (1.39) $\displaystyle=I\otimes
A^{(2)}+\left(S^{(3)}A^{(3)}\right)\otimes\left(\mathbf{1}\mathbf{1}^{T}\right),$

where $\otimes$ denotes the Kronecker product,131313A key property of the
Kronecker product is that $(A\cdot B)\otimes(C\cdot D)=(A\otimes
C)\cdot(B\otimes D)$. $I$ is the identity matrix, $\mathbf{1}$ is the all-ones
vector, and $S^{(3)}=\left(\begin{array}[]{ccc}0&0&0\\\ 1&0&0\\\
0&1&0\end{array}\right)$ is a non-cyclic shift matrix. Equation 1.28 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

$\displaystyle A^{(6)}$ $\displaystyle=I\otimes
A^{(2)}+\left(S^{(3)}A^{(3)}\right)\otimes\left(\mathbf{1}\mathbf{1}^{T}\right)$
(Eq. 1.28) $\displaystyle=(I\cdot I)\otimes(B^{(2)}\cdot
C^{(2)})+\left(S^{(3)}B^{(3)}\cdot
C^{(3)}\right)\otimes\left(\mathbf{1}\cdot\mathbf{1}^{T}\right)$
$\displaystyle=(I\otimes B^{(2)})\cdot(I\otimes
C^{(2)})+\left(S^{(3)}B^{(3)}\otimes\mathbf{1}\right)\cdot\left(C^{(3)}\otimes\mathbf{1}^{T}\right)$
$\displaystyle=\underbrace{\left(I\otimes
B^{(2)}~{}\mid~{}S^{(3)}B^{(3)}\otimes\mathbf{1}\right)}_{=B^{(6)}}\cdot\underbrace{\left(\begin{array}[]{c}I\otimes
C^{(2)}\\\ C^{(3)}\otimes\mathbf{1}^{T}\end{array}\right)}_{=C^{(6)}}.$ (1.42)

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 Equation
1.28 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
$\ell-1$ times to obtain a factorization of size $2^{\ell}$. Rather than
starting with a factorization of size $2$, we can start with a larger
factorization of size $n_{1}$ given by Theorem 1.1 and then repeat the
recursive step above $\ell-1$ times to obtain a factorization of size
$n_{1}^{\ell}$. 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 Theorem 1.2.

This recursive construction attains excellent asymptotics, but, for practical
parameter regimes, we find that the BLT approach is practically
indistinguishable from the optimal Toeplitz mechanism.

###  1.3 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.

Figure 2 compares mechanisms (and lower bounds) in terms of
$\operatorname{MaxErr}$ 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 [FHU22] provide the optimal (but inefficient)
construction. Our BLT mechanisms essentially match this performance, while
requiring time and memory $\widetilde{\mathcal{O}}(1)$ instead of
$\mathcal{O}(n)$.

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

It is of course possible to run Opt-BLT mechanisms for a different number of
steps than the optimization targeted. Figure 3 (Left column) explores this. We
construct three fixed Opt-BLT factorizations, optimized for
$n\in\\{100,1000,10000\\}$, and compare their performance (in terms of
$\operatorname{MaxErr}$ 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 $n$s for which they were optimized; however, the excess error is
highly asymetric; a mechanism optimized for $n^{*}$ will generally perform
well for $n<n^{*}$ steps, but can quickly perform very badly when $n>n^{*}$.
This is expected when one considers that Opt-BLT should intuitively be
ensuring a good approximation of the optimal Toeplitz coefficients
$r_{0},\dots,r_{n^{*}-1}$, but the approximation of the optimal Toeplitz
coefficients for larger $n$ can become arbitrarily bad. We see this in Figure
3 (Middle column), where we compare the Toeplitz coefficients defining $C$ and
$B$ to the optimal coefficients corresponding to the generating function
$1/\sqrt{1-x}$. To emphasize this issue, we optimize for $n^{*}=100$, and
consider degree $d=2$. The Opt-BLT factorization provides a better
approximation to the optimal coefficients for $n$ up to 100 compared to
$\text{{RA-BLT}}(d\\!=\\!2)$ (which does not depend on $n$), and a generally
worse approximation beyond that. Figure 3 (Right column) shows that while both
Opt-BLT and RA-BLT correspond to “reasonable” approximations to $\sqrt{1-x}$
(top), they distribute their errors in the approximation very differently
(bottom).

Figure 2:  Comparison of known upper and lower bounds for factorizations
$A=BC$ of the all-ones lower triangular matrix $A_{i,j}=\mathbb{I}[i\geq j]$.
Note that this includes non-Toeplitz factorizations. This illustrates that
there is a small gap between lower triangular Toeplitz factorizations and
general factorizations; furthermore this gap is asymptotically constant. Left:
Vertical axis is $\operatorname{MaxErr}(B,C)$. Right: Vertical axis is
$\operatorname{MaxErr}(B,C)-\mathsf{OptLTToe}(n)$.  Figure 3: (Left column)
Comparison of three fixed Opt-BLT mechanisms across a range of $n$ (extending
beyond the optimization targets). (Center) Comparing the first 200 Toeplitz
coefficients defining the $C$ and $B$ matrices for Opt-BLT and RA-BLT for
degree $d=2$, with the Opt-BLT mechanism optimized for $n^{*}=100$. (Right)
Differences in the generating functions for the Opt-BLT and RA-BLT
factorizations of the middle column.

##  2 Background

We start by formally setting up the problem where each stream entry is a
vector of dimension $m$ (as opposed to a scalar in the introduction). Fix a
stream length $n\in\mathbb{N}$ and let

$A=A^{(n)}\coloneqq\left(\begin{array}[]{ccccc}1&0&0&\cdots&0\\\
1&1&0&\cdots&0\\\ 1&1&1&\cdots&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
1&1&1&\cdots&1\end{array}\right)\in\\{0,1\\}^{n\times n}$ (2.1)

be the all-ones lower triangular matrix given by $A^{(n)}_{i,j}=1\iff i\geq j$
and $A_{i,j}=0\iff i<j$ for all $i,j\in[n]$. Note that $A$ is invertible and
its inverse is a lower triangular bi-diagonal matrix:

$A^{-1}=\left(\begin{array}[]{ccccc}1&0&0&\cdots&0\\\ -1&1&0&\cdots&0\\\
0&-1&1&\cdots&0\\\ \vdots&\ddots&\ddots&\ddots&\vdots\\\
0&\cdots&0&-1&1\end{array}\right)\in\\{-1,0,1\\}^{n\times n}.$ (2.2)

Our goal is to compute a matrix factorization $A=BC$ where
$B,C^{T}\in\mathbb{R}^{n\times n^{\prime}}$ that minimizes

$\operatorname{MaxErr}(B,C)\coloneqq\|B\|_{2\to\infty}\cdot\|C\|_{1\to
2}\coloneqq\sqrt{\max_{i\in[n]}\sum_{j\in[n^{\prime}]}B_{i,j}^{2}}\cdot\sqrt{\max_{j\in[n]}\sum_{i\in[n^{\prime}]}C_{i,j}^{2}}.$
(2.3)

Simultaneously we want computational efficiency in the following sense. We
need to generate samples from $BZ\in\mathbb{R}^{n\times m}$ where
$Z\in\mathbb{R}^{n^{\prime}\times m}$ 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 $m$ and constant or logarithmic in $n$.

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., $n^{\prime}=n$) 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},\cdots,x_{k}$ and not future inputs $x_{k+1},\cdots,x_{n-1}$.
This is particularly valuable for the privacy analysis in the case where the
input stream $x_{1},\cdots,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},\cdots,(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 [DMRST22].

###  2.1 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 $\operatorname{MaxErr}(B,C)$ with little regard for the
computational efficiency of sampling the correlated noise $Bz$ from i.i.d.
seed noise $z\leftarrow\mathcal{N}(0,\sigma^{2}I)$.

The binary tree mechanism [DNPR10, CSS11] implicitly constructs a
factorization which optimizes $\operatorname{MaxErr}(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)\in\\{0,1\\}^{1\times 1}$ and

$B^{(2n)}=\left(\begin{array}[]{ccc}B^{(n)}&0&0\\\
0&B^{(n)}&\vec{1}\end{array}\right)\in\\{0,1\\}^{2n\times(4n-1)}~{}~{}\text{
and }~{}~{}C^{(2n)}=\left(\begin{array}[]{cc}C^{(n)}&0\\\ 0&C^{(n)}\\\
\vec{1}^{T}&0\end{array}\right)\in\\{0,1\\}^{(4n-1)\times 2n}.$ (2.4)

Note that this only gives a factorization of $A^{(n)}$ when $n$ is a power of
$2$. By discarding rows/columns we can extend this to any $n$. The objective
value is given by

$\operatorname{MaxErr}(B,C)=\|B^{(2^{\ell})}\|_{2\to\infty}^{2}=\|C^{(2^{\ell})}\|_{1\to
2}^{2}=\ell+1=\left\lceil\frac{\log n}{\log 2}\right\rceil+1.$ (2.5)

The recursive formulation naturally leads to an efficient sampling algorithm
for the correlated noise $Bz$ using only $O(\ell)$ space. Subsequent work has
attempted to improve the constants [QYL13, Hon15, AP24].

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 BLTs.

[FHU22] give the following elegant explicit lower triangular Toeplitz
factorization. Let $f_{0}=1$ and, for $k\in\mathbb{Z}_{\geq 1}$, let
$f_{k}=f_{k-1}\cdot(1-1/2k)$ and $f_{-k}=0$. Equivalently,
$f_{k}=4^{-k}{2k\choose k}\leq\frac{1}{\sqrt{\pi k}}$ for $k\geq 1$ [elk].
Define $B,C\in\mathbb{R}^{n\times n}$ by $B_{i,j}=C_{i,j}=f_{i-j}$ for all
$i,j\in[n]$, which we denote $B=C=\operatorname{M}((f_{i})_{i=0}^{n-1})$ or
simply $\operatorname{M}(f,n)$. Then, we have that $A=BC$
and141414Unfortunately there appears to be an off-by-one error in the proof of
the bound given by Theorem 1 of [FHU22]; this has since been corrected.

$\operatorname{MaxErr}(B,C)=\|B\|_{2\to\infty}^{2}=\|C\|_{1\to
2}^{2}=\sum_{k=0}^{n-1}f_{k}^{2}\leq 1+\frac{\log(n)+\gamma}{\pi},$ (2.6)

where $\gamma\leq 0.57722$ is the Euler-Mascheroni constant [SW]. This
factorization turns out to be optimal among the class of lower triangular
Toeplitz factorizations; we prove this in Proposition 2.2.

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 $\gamma_{2}$ factorization norm is in fact a norm, i.e., it satisfies the
triangle inequality [TJ89]. Furthermore, there is a dual characterization of
the $\gamma_{2}$ factorization norm [LSŠ08, Theorem 9].

$\displaystyle\gamma_{2}(A)$
$\displaystyle\coloneqq\inf\\{\operatorname{MaxErr}(B,C):B,C\in\mathbb{R}^{n\times
n},A=BC\\}$ $\displaystyle=\sup\left\\{\|P^{1/2}\cdot A\cdot
Q^{1/2}\|_{\text{trace}}:\begin{aligned} P,Q\text{ non-negative diagonal
matrices}\\\ \text{ with
}\mathsf{trace}(P)=\mathsf{trace}(Q)=1\end{aligned}\right\\},$ (2.7)

where $\|M\|_{\text{trace}}=\mathsf{trace}(\sqrt{M^{T}M})$ is the sum of the
absolute singular values.

Using the triangle inequality, Mathias [Mat93, Corollary 3.5] gives a non-
constructive upper bound

$\gamma_{2}(A)\leq\frac{1}{2}+\frac{1}{2n}\sum_{j=1}^{n}\frac{1}{\sin\left(\pi\frac{2j-1}{2n}\right)}\leq
1+\frac{\log(n)}{\pi}.$ (2.8)

On the other hand, using the dual characterization (Eq. 2.7 with
$P=Q=\tfrac{1}{n}I$), [MNT20] give a lower bound151515Proposition 4.1 of
[MNT20] simply states an $\Omega(\log n)$ lower bound, but these sharper
expressions can easily be extracted from the proof. [Mat93] also gives a lower
bound which differs from the upper bound in Eq. 2.8 by less than an additive
$\frac{1}{2}$, but this seems slightly weaker than Equation 2.9.

$\gamma_{2}(A)\geq\frac{1}{2n}\sum_{j=1}^{n}\frac{1}{\sin\left(\pi\frac{2j-1}{4n+2}\right)}\geq\frac{\log(2n+1)}{\pi}.$
(2.9)

The upper bound in Eq. 2.6 or Eq. 2.8 and the lower bound in Eq. 2.9 match up
to a small _additive_ constant. Thus the optimal value of the matrix
factorization objective is $\gamma_{2}(A)=\frac{\log
n}{\pi}\pm\mathsf{constant}$. Numerically, the maximum gap between the
Toeplitz upper bound in Eq. 2.6 and the lower bound in Eq. 2.9 is less than
$0.365$.

Figure 2 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
$\frac{\log n}{\log 2}$, while the optimal leading term is $\frac{\log
n}{\pi}$. Thus the binary tree is asymptotically suboptimal by a
multiplicative factor of $\pi/\log 2\approx 4.5$.

###  2.2 Lower Triangular Toeplitz Factorizations versus General
Factorizations

Our work mainly focuses on lower triangular Toeplitz factorizations. In
particular, Theorem 1.1 gives such a factorization and proves near-optimality
with respect to this class. (Although we consider non-square factorizations in
Section 6, where we prove Theorem 1.2.) This structure is essential to our
algorithms in Sections 3, 4 and 5. However, it is natural to wonder how much
we lose in terms of the objective $\operatorname{MaxErr}(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 $\gamma_{2}(A)$. For lower triangular Toeplitz matrix
factorizations the best factorization we have is that of [FHU22]. 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
[FHU22] and we give a proof for completeness.

######  Lemma 2.1.

Define a sequence $f_{0},f_{1},\cdots$ by $f_{0}=1$ and
$f_{k}=(1-1/2k)f_{k-1}$ for all $k\geq 1$. Then, for all $k\geq 1$, we have

$f_{k}=4^{-k}{2k\choose
k}\in\left[\frac{1}{\sqrt{\pi(k+1)}},\frac{1}{\sqrt{\pi k}}\right].$ (2.10)

For all integers $n\geq 0$, we have

$\sum_{k=0}^{n}f_{k}f_{n-k}=1$ (2.11)

and

$\frac{\gamma+\log(n)-1}{\pi}\leq\sum_{k=1}^{\infty}f_{k}^{2}\leq\frac{\gamma+\log(n)}{\pi},$
(2.12)

where $0.57721\leq\gamma\leq 0.57722$ is the Euler-Mascheroni constant.
Furthermore, for all $x\in\mathbb{C}$ with $|x|<1$, we have

$\sum_{k=0}^{\infty}f_{k}x^{k}=\frac{1}{\sqrt{1-x}}.$ (2.13)

###### Proof.

The first part of Eq. 2.10 can be shown by induction: For $k\geq 1$,

$4^{-k}{2k\choose k}=\frac{1}{4}\cdot 4^{-(k-1)}{2k-2\choose
k-1}\frac{2k(2k-1)}{k^{2}}=\frac{2k(2k-1)}{4k^{2}}f_{k-1}=\left(1-\frac{1}{2k}\right)f_{k-1}=f_{k}.$

Next consider the derivatives of $f(x)\coloneqq\frac{1}{\sqrt{1-x}}$: For
$k\geq 1$, we have

$f^{(k)}(x)=(1-x)^{-k-1/2}\prod_{\ell=0}^{k-1}\left(\ell+\frac{1}{2}\right).$
(2.14)

Thus

$\frac{f^{(k)}(0)}{k!}=\prod_{\ell=0}^{k-1}\frac{1}{\ell+1}\left(\ell+\frac{1}{2}\right)=\prod_{\ell=0}^{k-1}\frac{(2\ell+1)}{2(\ell+1)}\frac{(2\ell+2)}{2(\ell+1)}=\frac{(2k)!}{(k!)^{2}2^{2k}}=4^{-k}{2k\choose
k}=f_{k}.$

By Taylor’s theorem, for all $x\in\mathbb{C}$ with $|x|<1$,

$\frac{1}{\sqrt{1-x}}=f(x)=\sum_{k=0}^{\infty}\frac{f^{(k)}(0)}{k!}x^{k}=\sum_{k=0}^{\infty}f_{k}x^{k},$

as required to prove Eq. 2.13. Furthermore, for all $x\in\mathbb{C}$ with
$|x|<1$, we have

$\sum_{n=0}^{\infty}x^{n}=\frac{1}{1-x}=f(x)^{2}=\left(\sum_{k=0}^{\infty}f_{k}x^{k}\right)^{2}=\sum_{n=0}^{\infty}x^{n}\sum_{k=0}^{n}f_{k}f_{n-k}.$

Matching coefficients proves Eq. 2.11

To prove the second part of Eq. 2.10, we use standard bounds on the central
binomial coefficient [elk, spe]: For all integers $k\geq 1$,

$\frac{4^{k}}{\sqrt{\pi(k+1)}}\leq{2k\choose k}\leq\frac{4^{k}}{\sqrt{\pi
k}},$ (2.15)

whence $\frac{1}{\sqrt{\pi(k+1)}}\leq f_{k}\leq\frac{1}{\sqrt{\pi k}}$ . Next
we use standard bounds on the harmonic numbers to prove Eq. 2.12: For all
integers $n\geq 1$,

$\log(n)+\gamma+\frac{1}{2(n+1)}\leq\sum_{k=1}^{n}\frac{1}{k}\leq\log(n)+\gamma+\frac{1}{2n}$
(2.16)

where $0.57721\leq\gamma\leq 0.57722$ is the Euler-Mascheroni constant [SW].
It follows that, for all $n\geq 2$,

$\sum_{k=1}^{n-1}f_{k}^{2}\leq\sum_{k=1}^{n-1}\frac{1}{\pi
k}\leq\frac{\gamma+\log(n-1)+1/2(n-1)}{\pi}\leq\frac{\gamma+\log(n)}{\pi}.$
(2.17)

Similarly, for all $n\geq 2$

$\sum_{k=1}^{n-1}f_{k}^{2}\geq\sum_{k=1}^{n-1}\frac{1}{\pi(k+1)}=\sum_{k=2}^{n}\frac{1}{\pi
k}\geq\frac{\gamma+\log(n)+1/(2n+2)-1}{\pi}.$ (2.18)

∎

Most importantly, we show that the factorization of [FHU22] is optimal among
the class of lower triangular Toeplitz factorizations:

######  Proposition 2.2 (Optimal lower triangular Toeplitz factorization).

For any integer $n\geq 1$, consider the optimization problem

$\begin{array}[]{rrclcl}\displaystyle\min_{{b},{c}\,\in\,\mathbb{R}^{n}}&\lx@intercol\left\|{b}\right\|_{2}\left\|{c}\right\|_{2}\hfil\lx@intercol\\\
\textrm{s.t.}&\forall k<n~{}~{}\sum_{i=0}^{k}b_{i}c_{k-i}=1.\\\ \end{array}$
(2.19)

The minimum is achieved at ${b}={c}=(f_{k})_{k=0}^{n-1}$, where
$f_{0},f_{1},\cdots$ is given by $f_{0}=1$ and
$f_{k}=\left({1-\frac{1}{2k}}\right)f_{k-1}$ for all $k\geq 1$, which are the
coefficients derived by [FHU22].

###### Proof.

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
${b}={c}=(f_{k})_{k=0}^{n-1}$.

We use the following properties of $f_{k}$’s established by [FHU22] and also
shown in Lemma 2.1:

  * •

$f_{k}>0\quad\forall k\in\\{0,\ldots,n-1\\}$.

  * •

Let $F\in\mathbb{R}^{n\times n}$ denote the lower triangular Toeplitz matrix
whose first column is $f_{0},f_{1},\ldots,f_{n-1}$. Then, we have $F^{2}=A$.

  * •

$\sum_{k=0}^{\infty}f_{k}x_{k}=\frac{1}{\sqrt{1-x}}\quad\forall x\in[0,1]$.

#### Reformulation as a quadratically constrained quadratic program:

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

$\displaystyle\left\|{b}\right\|_{2}\left\|{c}\right\|_{2}=\inf\left\\{\frac{\nu}{2}\left\|b\right\|_{2}^{2}+\frac{1}{2\nu}\left\|c\right\|_{2}^{2}:\nu>0\right\\}=\min\left\\{\frac{1}{2}\left\|\nu
b\right\|_{2}^{2}+\frac{1}{2}\left\|\frac{c}{\nu}\right\|_{2}^{2}:\nu\in\mathbb{R}\setminus\\{0\\}\right\\}.$

Furthermore, for any feasible $(b,c)$ that satisfy the constraints of Eq. 2.19
and any $\nu\neq 0$, the pair $(\nu b,\frac{c}{\nu})$ also satisfies the
constraints. If $S$ denotes the feasible set of Eq. 2.19, we have

$\min_{{b},{c}\in
S}\left\|{b}\right\|_{2}\left\|{c}\right\|_{2}=\min_{{b},{c}\in
S}\min_{\nu\in\mathbb{R}\setminus\\{0\\}}\frac{1}{2}\left\|\nu
b\right\|_{2}^{2}+\frac{1}{2}\left\|\frac{c}{\nu}\right\|_{2}^{2}=\min_{{b},{c}\in
S}\frac{1}{2}\left\|{b}\right\|_{2}^{2}+\frac{1}{2}\left\|{c}\right\|_{2}^{2}.$

Thus, we can simply write the overall optimization problem as

$\begin{array}[]{rrclcl}\displaystyle\min_{{b},{c}\in\mathbb{R}^{n}}&\lx@intercol\frac{1}{2}{\left({\left\|{b}\right\|_{2}^{2}+\left\|{c}\right\|_{2}^{2}}\right)}\hfil\lx@intercol\\\
\textrm{s.t.}&\forall k<n~{}~{}\sum\limits_{i=0}^{k}b_{i}c_{k-i}=1.\\\
\end{array}$ (2.20)

#### Lagrangian relaxation:

We write the Lagrangian relaxation of this optimization problem introducing
dual variables $\lambda_{0},\cdots,\lambda_{n-1}\in\mathbb{R}$ corresponding
to the constraints:

$\displaystyle
L\left({b,c,\lambda}\right)=\frac{1}{2}\left({\left\|{b}\right\|_{2}^{2}+\left\|{c}\right\|_{2}^{2}}\right)-\sum_{k=0}^{n-1}\lambda_{k}\left({\sum_{i=0}^{k}b_{i}c_{k-i}-1}\right)$
(2.21)
$\displaystyle=\frac{1}{2}\left({\left\|{b}\right\|_{2}^{2}+\left\|{c}\right\|_{2}^{2}}\right)-{{b}}^{\top}\Gamma\left({\lambda}\right){c}+{\mathbf{1}}^{\top}\lambda$
$\displaystyle=\frac{1}{2}{\begin{pmatrix}b\\\
c\end{pmatrix}}^{\top}\begin{pmatrix}I&-\Gamma\left({\lambda}\right)\\\
-\Gamma\left({\lambda}\right)&I\end{pmatrix}\begin{pmatrix}b\\\
c\end{pmatrix}+{\mathbf{1}}^{\top}\lambda,$

where $\Gamma\left({\lambda}\right)$ is the symmetric Hankel matrix

$\Gamma\left({\lambda}\right)\coloneqq\begin{pmatrix}\lambda_{0}&\lambda_{1}&\lambda_{2}&\ldots&\lambda_{n-1}\\\
\lambda_{1}&\lambda_{2}&\lambda_{3}&\ldots&0\\\
\lambda_{2}&\lambda_{3}&\lambda_{4}&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
\lambda_{n-1}&0&0&\ldots&0\end{pmatrix}.$

The constrained optimization problem in Eq. 2.20 is _equivalent_ to the
unconstrained min-max problem

$\min_{b,c\in\mathbb{R}^{n}}\max_{\lambda\in\mathbb{R}^{n}}L(b,c,\lambda).$
(2.22)

Letting $\vec{f}=(f_{0},\ldots,f_{n-1})\in\mathbb{R}^{n}$, we will set
$b=c=\vec{f}$ and exhibit a setting $\lambda=\lambda^{*}$ such that
$\nabla_{b,c,\lambda}L(b,c,\lambda)=0$ and $L(b,c,\lambda)$ is jointly convex
in $b,c$ and concave in $\lambda$. It follows that $b=c=\vec{f}$ is the
optimal solution to the original problem, as required. Since $b=c=\vec{f}$ is
feasible,
$\nabla_{\lambda}L(b,c,\lambda)=\nabla_{\lambda}L(\vec{f},\vec{f},\lambda)=0$.
The Lagrangian is linear in $\lambda$ and hence trivially concave. The
Lagrangian is convex in $\left({{b},{c}}\right)$ if and only if

$\begin{pmatrix}I&-\Gamma\left({\lambda}\right)\\\
-\Gamma\left({\lambda}\right)&I\end{pmatrix}\succeq 0$

or equivalently (by Schur complements) if
$\left\|\Gamma\left({\lambda}\right)\right\|_{\star}\leq 1$ where
$\left\|Q\right\|_{\star}\coloneqq\sup\\{\|Qx\|_{2}:\|x\|_{2}\leq 1\\}$
denotes the operator norm.

#### Construction of primal-dual optimal solution:

Now it only remains to compute $\lambda=\lambda^{*}$ such that
$\nabla_{b,c}L(b,c,\lambda)=0$ for $b=c=\vec{f}$ and
$\left\|\Gamma\left({\lambda}\right)^{2}\right\|_{\star}\leq 1$.

For all $k\in[n]$, we have $\frac{\partial}{\partial
b_{k}}L(b,c,\lambda)=b_{k}-\sum_{i=k}^{n-1}\lambda_{i}c_{i-k}$ and
$\frac{\partial}{\partial
c_{k}}L(b,c,\lambda)=c_{k}-\sum_{i=k}^{n-1}\lambda_{i}b_{i-k}$. Equivalently,
$\nabla_{b}L(b,c,\lambda)=b-\Gamma(\lambda)c$ and
$\nabla_{c}L(b,c,\lambda)=c-\Gamma(\lambda)b$. Thus, in order to ensure
$\nabla_{b,c}L(b,c,\lambda)=0$ when $b=c=\vec{f}$, it suffices to set
$\lambda=\lambda^{\star}$, where $\lambda^{\star}$ is chosen so as to solve

$\displaystyle\vec{f}=\Gamma\left({\lambda^{\star}}\right)\vec{f}$
$\displaystyle\iff P\vec{f}=P\Gamma\left({\lambda^{\star}}\right)\vec{f}$
$\displaystyle\iff\begin{pmatrix}f_{n-1}\\\ f_{n-2}\\\ \vdots\\\
f_{0}\end{pmatrix}=\begin{pmatrix}\lambda^{\star}_{n-1}&0&0&\ldots&0\\\
\lambda^{\star}_{n-2}&\lambda^{\star}_{n-1}&0&\ldots&0\\\
\lambda^{\star}_{n-3}&\lambda^{\star}_{n-2}&\lambda^{\star}_{n-1}&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
\lambda^{\star}_{0}&\lambda^{\star}_{1}&\lambda^{\star}_{2}&\ldots&\lambda^{\star}_{n-1}\end{pmatrix}\begin{pmatrix}f_{0}\\\
f_{1}\\\ \vdots\\\ f_{n-1}\end{pmatrix}$
$\displaystyle\iff\begin{pmatrix}f_{n-1}\\\ f_{n-2}\\\ \vdots\\\
f_{0}\end{pmatrix}=\begin{pmatrix}f_{0}&0&0&\ldots&0\\\
f_{1}&f_{0}&0&\ldots&0\\\ f_{2}&f_{1}&f_{0}&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
f_{n-1}&f_{n-2}&f_{n-3}&\ldots&f_{0}\end{pmatrix}\begin{pmatrix}\lambda^{\star}_{n-1}\\\
\lambda^{\star}_{n-2}\\\ \vdots\\\ \lambda^{\star}_{0}\end{pmatrix},$

where $P$ is the permutation matrix such that $Px$ is the reversal of the
vector $x$. Thus, we can solve for $\lambda^{\star}$ to obtain

$P\lambda^{\star}=\begin{pmatrix}\lambda^{\star}_{n-1}\\\
\lambda^{\star}_{n-2}\\\ \vdots\\\
\lambda^{\star}_{0}\end{pmatrix}=\begin{pmatrix}f_{0}&0&0&\ldots&0\\\
f_{1}-f_{0}&f_{0}&0&\ldots&0\\\ f_{2}-f_{1}&f_{1}-f_{0}&f_{0}&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
f_{n-1}-f_{n-2}&f_{n-2}-f_{n-3}&f_{n-3}-f_{n-4}&\ldots&f_{0}\end{pmatrix}\begin{pmatrix}f_{n-1}\\\
f_{n-2}\\\ \vdots\\\ f_{0}\end{pmatrix},$

where we used the fact that $\vec{f}$ is a feasible solution:

$\displaystyle\begin{pmatrix}f_{0}&0&0&\ldots&0\\\ f_{1}&f_{0}&0&\ldots&0\\\
f_{2}&f_{1}&f_{0}&\ldots&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\
f_{n-1}&f_{n-2}&f_{n-3}&\ldots&f_{0}\end{pmatrix}\begin{pmatrix}f_{0}&0&0&\ldots&0\\\
f_{1}&f_{0}&0&\ldots&0\\\ f_{2}&f_{1}&f_{0}&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
f_{n-1}&f_{n-2}&f_{n-3}&\ldots&f_{0}\end{pmatrix}=\begin{pmatrix}1&0&0&\ldots&0\\\
1&1&0&\ldots&0\\\ 1&1&1&\ldots&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\
1&1&1&\ldots&1\end{pmatrix}$
$\displaystyle\iff{\begin{pmatrix}f_{0}&0&0&\ldots&0\\\
f_{1}&f_{0}&0&\ldots&0\\\ f_{2}&f_{1}&f_{0}&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
f_{n-1}&f_{n-2}&f_{n-3}&\ldots&f_{0}\end{pmatrix}}^{-1}={\begin{pmatrix}1&0&0&\ldots&0\\\
1&1&0&\ldots&0\\\ 1&1&1&\ldots&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\
1&1&1&\ldots&1\end{pmatrix}}^{-1}\begin{pmatrix}f_{0}&0&0&\ldots&0\\\
f_{1}&f_{0}&0&\ldots&0\\\ f_{2}&f_{1}&f_{0}&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
f_{n-1}&f_{n-2}&f_{n-3}&\ldots&f_{0}\end{pmatrix}$
$\displaystyle\iff{\begin{pmatrix}f_{0}&0&0&\ldots&0\\\
f_{1}&f_{0}&0&\ldots&0\\\ f_{2}&f_{1}&f_{0}&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
f_{n-1}&f_{n-2}&f_{n-3}&\ldots&f_{0}\end{pmatrix}}^{-1}={\begin{pmatrix}1&0&0&\ldots&0\\\
-1&1&0&\ldots&0\\\ 0&-1&1&\ldots&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&\ldots&1\end{pmatrix}}\begin{pmatrix}f_{0}&0&0&\ldots&0\\\
f_{1}&f_{0}&0&\ldots&0\\\ f_{2}&f_{1}&f_{0}&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
f_{n-1}&f_{n-2}&f_{n-3}&\ldots&f_{0}\end{pmatrix}.$

From the definition $f_{0}=1$ and $f_{k}=f_{k-1}(1-1/2k)$, for all $k\geq 1$,
we have

$\displaystyle\lambda^{\star}_{n-1-i}$
$\displaystyle=f_{0}f_{n-1-i}-\sum_{j=1}^{i}\left({f_{j}-f_{j-1}}\right)f_{n-1-i+j}$
$\displaystyle=f_{n-1-i}-\sum_{j=1}^{i}f_{j-1}\left({\frac{1}{2j}}\right)f_{n-1-i+j}$
$\displaystyle=f_{n-1-i}\left({1-\sum_{j=1}^{i}f_{j-1}\left({\frac{1}{2j}}\right)\frac{f_{n-1-i+j}}{f_{n-1-i}}}\right).$

Since
$\frac{f_{n-1-i+j}}{f_{n-1-i}}=\prod_{\ell=1}^{j}\left({1-\frac{1}{2(n-i+\ell)}}\right)\leq
1$, we have

$\displaystyle\lambda^{\star}_{n-1-i}\geq
f_{n-1-i}\left({1-\sum_{j=1}^{i}f_{j-1}\left({\frac{1}{2j}}\right)}\right)>f_{n-1-i}\left({1-\sum_{j=1}^{\infty}f_{j-1}\left({\frac{1}{2j}}\right)}\right)=f_{n-1-i}\left({1-\sum_{j=0}^{\infty}\frac{f_{j}}{2(j+1)}}\right).$
(2.23)

Since $\sum_{j=0}^{\infty}f_{j}x^{j}=\frac{1}{\sqrt{1-x}}$ for all
$x\in[0,1)$, integrating both sides between the limits $0$ and $t$, we have

$\sum_{j=0}^{\infty}\frac{f_{j}}{j+1}t^{j+1}=2\left({1-\sqrt{1-t}}\right).$

Taking limits as $t\to 1$, we obtain

$\sum_{j=0}^{\infty}\frac{f_{j}}{j+1}=2.$

Thus, from Eq. 2.23, we can conclude that $\lambda^{\star}_{n-1-i}>0$ for all
$i\in[n]$ and hence $\Gamma\left({\lambda^{\star}}\right)^{2}$ is a matrix
with all entries
$(\Gamma(\lambda^{\star})^{2})_{i,j}=\sum_{k=0}^{n-1-\max\\{i,j\\}}\lambda^{\star}_{i+k}\lambda^{\star}_{j+k}$
strictly positive. Hence, by the Perron-Frobenius theorem, we have that
$\Gamma\left({\lambda^{\star}}\right)^{2}$ has a unique eigenvector with all
coordinates positive, and this corresponds to the maximum eigenvalue. Since
$f$ is an eigenvector of $\Gamma\left({\lambda^{\star}}\right)^{2}$ by
construction, it must be the unique leading eigenvector corresponding to the
eigenvalue $1$. Hence
$\left\|\Gamma\left({\lambda^{\star}}\right)^{2}\right\|_{*}\leq 1$.

Thus, $L\left({{b},{c},\lambda^{\star}}\right)$ is convex in
$\left({{b},{c}}\right)$, since
$\left\|\Gamma\left({\lambda}\right)^{2}\right\|_{\star}\leq
1\iff\left\|\Gamma\left({\lambda}\right)\right\|_{\star}\leq 1$. And
$\lambda^{\star}$ was chosen to ensure $\nabla_{b,c}L(b,c,\lambda^{\star})=0$.
This lets us conclude that $\left({\vec{f},\vec{f},\lambda^{\star}}\right)$ is
a stationary point of the Lagrangian that is convex in ${b},{c}$ and concave
in $\lambda$, hence it is the global optimum.

∎

Proposition 2.2 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:

######  Definition 2.3.

Let $A\in\\{0,1\\}^{n\times n}$ be the lower triangular matrix of all ones.
Recall that $C\in\mathbb{R}^{n\times n}$ is lower triangular Toeplitz if
$C_{i,j}=0$ for all $i<j$ and $C_{i,j}=c_{i-j}$ for all $i\geq j$, where
$c\in\mathbb{R}^{n}$ is a vector. Define

$\displaystyle\mathsf{OptLTToe}(n)$
$\displaystyle\coloneqq\inf\left\\{\operatorname{MaxErr}(B,C):B,C\in\mathbb{R}^{n\times
n}\text{ lower triangular Toeplitz},BC=A\right\\},$ (2.24)
$\displaystyle\mathsf{Opt}(n)$
$\displaystyle\coloneqq\inf\left\\{\operatorname{MaxErr}(B,C):B,C\in\mathbb{R}^{n\times
n},BC=A\right\\},$ (2.25)

where $\operatorname{MaxErr}(B,C)\coloneqq\|B\|_{2\to\infty}\|C\|_{1\to
2}=\sqrt{\max_{i}\sum_{j}B_{i,j}^{2}}\sqrt{\max_{j}\sum_{i}C_{i,j}^{2}}$.

By Proposition 2.2 and Lemma 2.1, for all integers $n\geq 1$, we have

$1+\frac{\gamma-1+\log(n)}{\pi}\leq\mathsf{OptLTToe}(n)=1+\sum_{k=1}^{n-1}\left(2^{-2k}{2k\choose
k}\right)^{2}\leq 1+\frac{\gamma+\log(n)}{\pi}.$ (2.26)

The lower bound of [MNT20] in Eq. 2.9 states the optimal over all
factorizations is

$\mathsf{Opt}(n)=\gamma_{2}(A^{(n)})\geq 1+\frac{\log\left(2n+1\right)}{\pi}.$

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

######  Corollary 2.4.

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

$\mathsf{OptLTToe}(n)-\mathsf{Opt}(n)\leq\left(1+\frac{\gamma+\log(n)}{\pi}\right)-\left(1+\frac{\log(2n+1)}{\pi}\right)\leq\frac{\gamma+\log(2)}{\pi}\leq
0.405.$

Numerically, we have $\mathsf{OptLTToe}(n)-\mathsf{Opt}(n)\leq 0.365$ for all
$n\geq 1$. Figure 2 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=\hat{B}\hat{C}$ are functionally equivalent as long as
$BB^{T}=\hat{B}\hat{B}^{T}$. To see why this fact is true, note that for
$z\leftarrow\mathcal{N}(0,\sigma^{2}I)$, the distribution of $Bz$ is
$\mathcal{N}(0,\sigma^{2}BB^{T})$. So if $BB^{T}=\hat{B}\hat{B}^{T}$, then
$Bz$ and $\hat{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 $\hat{B}$ satisfying
$BB^{T}=\hat{B}\hat{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.

###  2.3 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 [ST13, AS17, KMSTTX21, AFKT23, ZTC22], and in practical deployments
[KMSTTX21, DMRST22, CCMRT22, CCGMMRGTX24, XZACCKMRZ23]. 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 vectors161616There
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) [KMSTTX21, ST13] 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 [XZACCKMRZ23].

For a given set of $\ell_{2}$-Lipschitz loss functions
$\\{f_{1}(\theta),\ldots,f_{n}(\theta)\\}$ (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:

$\texttt{DP-FTRL:}\ \ \ \
\theta_{t}\leftarrow\theta_{0}-\eta\left(\sum_{k=0}^{t-1}\nabla
f_{k}(\theta_{k})+\hat{z}_{k}\right).$ (2.27)

In Eq. 2.27, $\eta$ is the learning rate, and $\\{\hat{z}_{t}\\}$ are the
correlated Gaussian noise added to ensure that the computation of the sequence
of $\\{\theta_{t}\\}$’s are differentially private. Mapping to the MF view in
Eq. 1.1, the adaptively chosen data set is set to $x=\\{\nabla
f_{t}(\theta_{t}):t\in[n]\\}$, and the noise sequence $\\{\hat{z}_{t}\\}$ is
set to $Bz$, where $z\leftarrow\mathcal{N}(0,\sigma^{2}I)$. 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 [SCS13,
BST14, ACGMMTZ16] (another heavily used private learning algorithm) as a
special case [CCGMMRGTX24], and has much better privacy/utility trade-offs in
general [CCDPGST24]. [KMSTTX21, AFKT23] showed that for obtaining the best
regret/SCO guarantees from DP-FTRL style algorithms, it suffices to optimize
for $\operatorname{MaxErr}(B,C)$ to instantiate the noise mechanism in Eq.
2.27.

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 Section 1. 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.)

###  2.4 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.171717The error of our algorithms depends on the number
of iterations $n$. 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

$\forall
x~{}~{}~{}{\underset{\mathcal{M}}{\mathbb{E}}\left[\|\mathcal{M}(x)-Ax\|_{\infty}\right]}={\underset{\mathcal{M}}{\mathbb{E}}\left[\max_{k}|(\mathcal{M}(x)-Ax)_{k}|\right]}\leq
O(\log(n))^{3/2},$ (2.28)

while the maximum expected error is

$\forall
x~{}~{}~{}\max_{k}{\underset{\mathcal{M}}{\mathbb{E}}\left[|(\mathcal{M}(x)-Ax)_{k}|\right]}\leq
O(\log(n)).$ (2.29)

The additional $\sqrt{\log(n)}$ term in Equation 2.28 compared to Equation
2.29 comes from taking a union bound over $n$ Gaussians. Note that, for
simplicity, we suppress the privacy parameters; to achieve
$(\varepsilon,\delta)$-differential privacy, the error above scales with
$\frac{1}{\varepsilon}\log(1/\delta)$; to achieve $\rho$-zCDP [BS16], the
error above scales with $1/\sqrt{\rho}$.

The bound in Eq. 2.29 on the maximum expected error is tight up to constants
in the streaming setting – in fact the lower bound holds for _average_
expected error:

$\forall\mathcal{M}~{}\exists
x~{}~{}~{}\frac{1}{n}\sum_{k}{\underset{\mathcal{M}}{\mathbb{E}}\left[|(\mathcal{M}(x)-Ax)_{k}|\right]}\geq\Omega(\log(n)),$
(2.30)

assuming $\mathcal{M}$ is differentially private [HUU23, Theorem 4]. However,
the bound in Eq. 2.28 on expected maximum error is _not_ tight if we move
beyond additive noise mechanisms [DNRR15]. [DNRR15] show that we can obtain an
improved error guarantee:181818[DNRR15] stated this result for pure
differential privacy, with $\log(|x|)^{2}$ instead of $\log(|x|)^{3/2}$.

$\exists\mathcal{M}~{}\forall
x~{}~{}~{}{\underset{\mathcal{M}}{\mathbb{E}}\left[\|\mathcal{M}(x)-Ax\|_{\infty}\right]}\leq
O(\log(n)+\log(|x|)^{3/2}),$ (2.31)

where $|x|$ is the number of data points.191919In 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 $\mathcal{M}$ 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 $\leq|x|$ updates, rather than $n$ updates.
Thus our improvements for additive noise mechanisms can be applied in this
setting.

Eq. 2.31 is only an improvement over Eq. 2.28 in the sparse setting ($|x|\ll
n$). Machine learning applications are typically dense ($|x|\geq 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 [BNS13, BNSV15, BDRS18,
ALMM19, KLMNS20, GJK21, CLNSS23]. The offline setting permits more
sophisticated algorithms, which achieve asymptotically better error
guarantees. [CLNSS23] give an $(\varepsilon,\delta)$-differentially private
algorithm $\mathcal{M}$ with error

$\forall
x~{}~{}~{}{\underset{\mathcal{M}}{\mathbb{E}}\left[\|\mathcal{M}(x)-Ax\|_{\infty}\right]}\leq\widetilde{O}\left(\frac{1}{\varepsilon}\log^{2}(1/\delta)\log^{*}(n)\right)$
(2.32)

Note that offline algorithms are not applicable in the machine learning
applications discussed in Section 2.3.

#### Pure differential privacy:

The early work on continual counting [DNPR10, CSS11, XWG10] works under pure
$(\varepsilon,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\to\infty}\cdot\|C\|_{1\to 1}$, where $\|C\|_{1\to
1}=\max_{j}\sum_{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$, [Mat93] gives exact expressions for
$\gamma_{2}(A+A^{T})$ and $\gamma_{2}(A+A^{T}+I)$. There is a different line
of work [LMHMR15, MMHM21] 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:

$n$ | Number of steps on which DP estimates are released.
---|---
$d$ | Order of recurrence / degree / number of memory buffers.
$m$ | Dimension of per-step user contributions (e.g., model size).
$[k]$ |  $=\\{0,\dots,k-1\\}$ for $k\geq 1$.
$\mathbb{N}$ |  $=\\{0,1,2,\dots\\}$, the natural numbers including zero.
$f$ | Rational (generating) function where $f(x)=f_{0}+f_{1}x+f_{2}x^{2}\cdots$.
$p$ | Polynomial of degree $d$ with $p(x)=p_{0}+p_{1}x+p_{2}x^{2}+\cdots+p_{d}x^{d}$.
$(f_{k})_{k=0}^{\infty}$ | Sequence $f_{0},f_{1},\dots$ with corresponding generating function $f$.
$\operatorname{M}(f,n)\in\mathbb{R}^{n\times n}$ | Lower-triangular Toeplitz matrices defined by a generating function $f$.
$\cdot$ | Context-dependent multiplication, used where it improves readability.
$f*r$ | Cauchy product / convolution $(f*r)_{n}=\sum_{i=0}^{n}f_{i}\cdot r_{n-i}$.
$A^{T}$ | Transpose of a matrix $A$.
$A^{\star}$ | A matrix $A$ that is “optimal” in a context-dependent sense.
$A^{\dagger}$ | Moore-Penrose pseudoinverse of matrix $A$.
$A_{i,j}$ | The $(i,j)^{\text{th}}$ entry of matrix $A$, zero-indexed.
$A_{i,:}$ and $A_{:,j}$ | The $i^{\text{th}}$ row and $j^{\text{th}}$ column.
$\left\|A\right\|_{F}$ | The Frobenius norm of a matrix $A$.
$A\in\mathbb{R}^{n\times n}$ | Specifically, the lower-triangular matrix of ones to be factorized as $A=BC$.
$\log$ | Natural logarithm, $\log(2.718)\approx 1$.
Table 1: Summary of notation

##  3 Efficiently-Sampleable Factorization via Rational Generating Functions

In this section we describe our main algorithmic tools. First (in §3.1) we
present the view of lower triangular Toeplitz matrices in terms of generating
functions or sequences. This is a convenient mathematical formalism for
analysis [FHU22]. Second (in §3.2) 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 §3.3) we present our algorithm for sampling from lower triangular
Toeplitz matrices with rational generating functions (that is, the
multiplication algorithm for BLT matrices). In subsequent sections (§4, §5) we
instantiate these rational generating functions.

###  3.1 Sequences, Lower Triangular Toeplitz Matrices, & Generating Functions

We begin by describing the generating function view of lower triangular
Toeplitz matrices used by [FHU22]. Our analysis moves fluidly between three
different views of the same mathematical object:

  * •

The sequence $(f_{k})_{k=0}^{\infty}=f_{0},f_{1},f_{2},\dots$ with
$f_{k}\in\mathbb{R}$ for each index $k$.

  * •

The ordinary generating function of the sequence,
$f:\mathbb{C}\rightarrow\mathbb{C}$,

$f(x)=\sum_{k=0}^{\infty}f_{k}x^{k}=\sum_{k=0}^{\infty}\frac{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 ${n\times n}$ lower-triangular Toeplitz matrices
$M(f,n)\in\mathbb{R}^{n\times n}$ for $n\in\mathbb{N}$ that are generated by
$f$ as

$\forall
i,j\in[n]~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}M(f,n)_{i,j}\coloneqq\left\\{\begin{array}[]{cl}f_{i-j}&\text{
if }i\geq j\\\ 0&\text{ if }i<j\end{array}\right\\}$ (3.1)

When the specific $n$ 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})_{k=0}^{\infty}=1,1,1,1,\dots,\quad
g(x)=\frac{1}{1-x}=\sum_{k=0}^{\infty}x^{k},\quad A^{(4)}\coloneqq
M(g,4)=\begin{bmatrix}1&0&0&0\\\ 1&1&0&0\\\ 1&1&1&0\\\ 1&1&1&1\\\
\end{bmatrix}$ (3.2)

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})_{k=0}^{\infty}*(g_{k})_{k=0}^{\infty}=(h_{k})_{k=0}^{\infty}\quad\text{where}\quad
h_{k}\coloneqq\sum_{i=0}^{k}f_{i}\cdot g_{k-i},$ (3.3)

and write $h=f*g$ as a shorthand.

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

######  Lemma 3.1.

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,

$\left(h=f*g\right)\quad\Longleftrightarrow\quad\big{(}h(x)=f(x)g(x)\big{)}\quad\Longleftrightarrow\quad\big{(}\forall
n~{}~{}M(h,n)=M(f,n)\cdot M(g,n)\big{)}.$

Similarly the usual definitions of addition are equivalent:

$\big{(}\forall
i~{}~{}h_{i}=f_{i}+g_{i}\big{)}\quad\Longleftrightarrow\quad\big{(}h(x)=f(x)+g(x)\big{)}\quad\Longleftrightarrow\quad\big{(}\forall
n~{}~{}M(h,n)=M(f,n)+M(g,n)\big{)}.$

###### Proof.

The results for sequences and their ordinary generating functions are well
known, see e.g. [KP11] or their equivalent results based on well-known
properties of the Z-transform; cf. [OWN97]. The results for addition are also
immediate. Hence, it is sufficient to show $h(x)=f(x)g(x)$ iff $\forall
n\in\mathbb{N},M(h,n)=M(f,n)\cdot M(g,n).$

Fix $i,j\in[n]$. If $i<j$, then $M(h,n)_{i,j}=0$ and $(M(f,n)\cdot
M(g,n))_{i,j}=\sum_{k\in[n]}M(f,n)_{i,k}\cdot M(g,n)_{k,j}=\sum_{k}0$, since,
for all $k\in[n]$, either $i<k$, whence $M(f,n)_{i,k}=0$, or $k<j$, whence
$M(g,n)_{k,j}=0$, (or both). Thus we can assume $i\geq j$.

By the product rule and induction, we have the derivative
$h^{(\ell)}(x)=\sum_{k=0}^{\ell}{\ell\choose k}f^{(\ell-k)}(x)\cdot
g^{(k)}(x)$ for all $\ell\geq 0$ and all applicable $x$. Thus

$\displaystyle M(h,n)_{i,j}$ $\displaystyle=\frac{h^{(i-j)}(0)}{(i-j)!}$
$\displaystyle=\frac{1}{(i-j)!}\sum_{k=0}^{i-j}{i-j\choose
k}f^{(i-j-k)}(0)\cdot g^{(k)}(0)$
$\displaystyle=\sum_{k=0}^{i-j}\frac{f^{(i-j-k)}(0)}{(i-j-k)!}\cdot\frac{g^{(k)}(0)}{k!}$
$\displaystyle=\sum_{k^{\prime}=j}^{i}\frac{f^{(i-k^{\prime})}(0)}{(i-k^{\prime})!}\cdot\frac{g^{(k^{\prime}-j)}(0)}{(k^{\prime}-j)!}$
($k^{\prime}=k+j$)
$\displaystyle=\sum_{k^{\prime}=j}^{i}M(f,n)_{i,k^{\prime}}\cdot
M(g,n)_{k^{\prime},j}$
$\displaystyle=\sum_{k^{\prime}\in[n]}M(f,n)_{i,k^{\prime}}\cdot
M(g,n)_{k^{\prime},j}$ ($i<k^{\prime}\implies M(f,n)_{i,k^{\prime}}=0$,
$k^{\prime}<j\implies M(g,n)_{k^{\prime},j}=0$)
$\displaystyle=\left(M(f,n)\cdot M(g,n)\right)_{i,j}.$

∎

Given Eq. 3.2 and the above result for multiplication, in order to factorize
$A=BC$, it is natural to factorize by taking the square root of the generating
function. Let $f(x)=\frac{1}{\sqrt{1-x}}$. Since $f(x)\cdot
f(x)=\frac{1}{1-x}=g(x)$ is the generating function of the all ones sequence,
it follows that $M(f,n)\cdot M(f,n)=M(g,n)=A^{(n)}$. Indeed, this is the
factorization given by [FHU22], which we show to be optimal among the class of
lower triangular Toeplitz matrices in Proposition 2.2.

Given any generating function $r(x)$, we can obtain a valid factorization
$A=BC$ where $B=M(b,n)$ for $b(x)=r(x)/(1-x)$ and $C=M(c,n)$ for
$c(x)=1/r(x)$. Our approach is to choose a generating function such that
$r(x)\approx\sqrt{1-x}$ but which also permits us to design an efficient
algorithm.

###  3.2 Rational Generating Functions and Constant-Recurrent Sequences

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

$r(x)=\frac{p(x)}{q(x)}=\frac{\sum_{i=0}^{d-1}p_{i}x^{i}}{1+\sum_{j=1}^{d}q_{j}x^{j}}.$
(3.4)

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 $\leq d$. Note that we normalize the denominator so that
$q(0)=q_{0}=1$.202020We 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.

Equation 3.4 is a convenient mathematical representation, but for algorithmic
purposes we need the sequence $(r_{k})_{k=0}^{\infty}$ of Taylor series
coefficients

$r(x)=\sum_{k=0}^{\infty}r_{k}x^{k}.$ (3.5)

Thus our first task is to map from Equation 3.4 to Equation 3.5 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})_{k=0}^{\infty}$ being a constant-recurrent sequence (a.k.a.
linear-recursive or C-finite sequence); see for example [KP11]. 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 Lemma 3.2 and
provide a proof for completeness.

######  Lemma 3.2 (Constant-recurrent Taylor series representation of a
rational function).

Let

$r(x)=\frac{p(x)}{q(x)}=\frac{\sum_{i=0}^{d-1}p_{i}x^{i}}{1+\sum_{j=1}^{d}q_{j}x^{j}}$

be a rational function of degree $\leq d$. As in Section 3.1, let

$r(x)=\sum_{k=0}^{\infty}r_{k}x^{k}.$

That is, $r(x)$ generates the sequence $(r_{k})_{k=0}^{\infty}$. Then, for all
$k\geq d$, this sequence satisfies the recurrence

$r_{k}=-\sum_{j=1}^{d}q_{j}r_{k-j}.$ (3.6)

And, for $0\leq k<d$, we have $r_{k}=p_{k}-\sum_{j=1}^{k}q_{j}r_{k-j}$.
Furthermore, for all $k$, we have

$r_{k}=u^{T}W^{k}v,$ (3.7)

where

$u\coloneqq\left(\begin{array}[]{c}1\\\ 0\\\ 0\\\ 0\\\ \vdots\\\ 0\\\
0\end{array}\right)\in\mathbb{R}^{d\times 1},\qquad
v\coloneqq\left(\begin{array}[]{c}p_{0}\\\ p_{1}\\\ p_{2}\\\ p_{3}\\\
\vdots\\\ p_{d-2}\\\ p_{d-1}\end{array}\right)\in\mathbb{R}^{d\times 1},$
(3.8)

and

$W\coloneqq\left(\begin{array}[]{ccccccc}-q_{1}&1&0&0&\cdots&0&0\\\
-q_{2}&0&1&0&\cdots&0&0\\\ -q_{3}&0&0&1&\cdots&0&0\\\
-q_{4}&0&0&0&\cdots&0&0\\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
-q_{d-1}&0&0&0&\cdots&0&1\\\
-q_{d}&0&0&0&\cdots&0&0\end{array}\right)\in\mathbb{R}^{d\times d}.$ (3.9)

###### Proof.

We have

$r(x)=\frac{p(x)}{q(x)}=\frac{\sum_{i=0}^{d-1}p_{i}x^{i}}{1+\sum_{j=1}^{d}q_{j}x^{j}}=\sum_{k=0}^{\infty}r_{k}x^{k},$

which rearranges to

$p(x)=\sum_{i=0}^{d-1}p_{i}x^{i}=q(x)\cdot
r(x)=\left(1+\sum_{j=1}^{d}q_{j}x^{j}\right)\cdot\left(\sum_{k=0}^{\infty}r_{k}x^{k}\right)=\sum_{k=0}^{\infty}\left(r_{k}+\sum_{j=1}^{\min\\{d,k\\}}q_{j}r_{k-j}\right)\cdot
x^{k}.$

Matching coefficients gives $p_{k}=r_{k}+\sum_{j=1}^{k}q_{j}r_{k-j}$ for
$1\leq k\leq d-1$ and $0=r_{k}+\sum_{j=1}^{d}q_{j}r_{k-j}$ for $k\geq 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\leq i\leq d-1$. For all $k\geq 1$, let
${p}_{j}^{(k)}\coloneqq p_{j+1}^{(k-1)}-p_{0}^{(k-1)}q_{j+1}$ for $0\leq j\leq
d-2$ and let ${p}_{d-1}^{(k)}\coloneqq-p_{0}^{(k-1)}q_{d}$. We can check that
these values satisfy

$\frac{\sum_{i=0}^{d-1}p_{i}^{(k-1)}x^{i}}{1+\sum_{j=1}^{d}q_{j}x^{j}}=p_{0}^{(k-1)}+x\cdot\frac{\sum_{i=0}^{d-1}{p}_{i}^{(k)}x^{i}}{1+\sum_{j=1}^{d}q_{j}x^{j}}.$
(3.10)

The left hand side of Equation 3.10 with $k=1$ is simply $r(x)$; from this we
can conclude that $r_{0}=p_{0}^{(0)}=p_{0}$. Expanding Equation 3.10 for $k=1$
and $k=2$ gives $r_{1}={p}_{0}^{(1)}$. Specifically, we have

$r(x)=\frac{\sum_{i=0}^{d-1}p_{i}^{(0)}x^{i}}{1+\sum_{j=1}^{d}q_{j}x^{j}}=p_{0}^{(0)}+x\cdot\frac{\sum_{i=0}^{d-1}{p}_{i}^{(1)}x^{i}}{1+\sum_{j=1}^{d}q_{j}x^{j}}=p_{0}^{(0)}+x\cdot\left(p_{0}^{(1)}+x\cdot\frac{\sum_{i=0}^{d-1}{p}_{i}^{(2)}x^{i}}{1+\sum_{j=1}^{d}q_{j}x^{j}}\right).$

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)=\sum_{\ell=0}^{k-1}p_{0}^{(\ell)}x^{\ell}+x^{k}\cdot\frac{\sum_{i=0}^{d-1}{p}_{i}^{(k)}x^{i}}{1+\sum_{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\geq 1$, we have

$\left(\begin{array}[]{c}{p}_{0}^{(k)}\\\ {p}_{1}^{(k)}\\\ \vdots\\\
{p}_{d-1}^{(k)}\end{array}\right)=W\cdot\left(\begin{array}[]{c}{p}_{0}^{(k-1)}\\\
{p}_{1}^{(k-1)}\\\ \vdots\\\ {p}_{d-1}^{(k-1)}\end{array}\right),$

where $W\in\mathbb{R}^{d\times d}$ is given in Equation 3.9. Given this,
induction shows that, for all $k\geq 0$,

$\left(\begin{array}[]{c}{p}_{0}^{(k)}\\\ {p}_{1}^{(k)}\\\ \vdots\\\
{p}_{d-1}^{(k)}\end{array}\right)=W^{k}\cdot\left(\begin{array}[]{c}{p}_{0}^{(0)}\\\
{p}_{1}^{(0)}\\\ \vdots\\\ {p}_{d-1}^{(0)}\end{array}\right)=W^{k}\cdot v$

and, hence, $r_{k}=p_{0}^{(k)}=u^{T}\cdot W^{k}\cdot v$, as required. ∎

We make some remarks about Lemma 3.2:

  1. 1.

Lemma 3.2 assumes $\deg(p)<d$ and $\deg(q)\leq 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.

  2. 2.

Lemma 3.2 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. [KP11]). However, we do not use the converse of the result.

  3. 3.

The eigenvalues of $W$ can be related to the poles of the corresponding
rational function. Specifically,

$\det(\lambda
I-W)=0\iff\lambda^{d}+\sum_{j=1}^{d}q_{j}\lambda^{d-j}=0\iff\big{(}q(\tfrac{1}{\lambda})=0\text{
or }\lambda=0=q_{d}\big{)}.$ (3.11)

  4. 4.

Furthermore, the eigenvectors of $W$ can be written in terms of the poles:
Suppose
$q(\tfrac{1}{\lambda})=1+\sum_{j=1}^{d}q_{j}(\tfrac{1}{\lambda})^{j}=0$. Let
$u^{T}=(\lambda^{d-1},\lambda^{d-2},\cdots,\lambda,1)\in\mathbb{R}^{1\times
d}$ (or $\lambda=0$ and $q_{d}=0$). Then $u^{T}W=\lambda\cdot u^{T}$ or,
equivalently, $W^{T}u=\lambda\cdot u$.

  5. 5.

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 $\mathbb{R}$). A diagonal matrix $W$ is particularly convenient for
our streaming algorithm (Algorithm 1).

  6. 6.

If the eigenvalues of $W$ are contained in the unit circle (i.e.,
$|\lambda|\leq 1$), then the computation of $W^{k}$ is numerically stable. By
Equation 3.11, this is equivalent to the poles of the rational function having
magnitude $\geq 1$.

  7. 7.

Suppose the rational function $r$ has distinct poles $\frac{1}{\theta_{i}}$,
and can be presented as212121Diagonalizing the matrix $W$ given by Lemma 3.2
corresponds to finding a representation of this form.

$r(x)=\sum_{i\in[d]}\frac{\omega_{i}}{1-\theta_{i}\cdot x}$

via a partial-fraction decomposition (instead of as a ratio of polynomials as
in Lemma 3.2). (Indeed, our rational function in Section 4.2 is presented in
this form in Eq. 4.9.) Then, we have a simpler closed form for the sequence,

$r_{k}=\sum_{i\in[d]}\omega_{i}\theta_{i}^{k}.$

This immediately yields a matrix representation: $W$ is a diagonal matrix with
$W_{i,i}=\theta_{i}$ for all $i$, and the vectors $u$ and $v$ simply need to
satisfy $u_{i}v_{i}=\omega_{i}$ for all $i$. This matrix representation may be
more convenient than the one given by Lemma 3.2; we use this approach
extensively in Section 5.

  8. 8.

The form of the matrix $W$ in Lemma 3.2 is known as a “companion matrix”
[Wik24]. 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 Eq. 3.11).

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

######  Remark 3.3.

The representation in Equation 3.4 assumes the numerator has lower degree than
the denominator, i.e., $\deg(p)<d$ versus $\deg(q)\leq d$. But, in this work,
we typically consider rational functions with the numerator and denominator
both having the same degree, i.e., $\deg(\overline{p})=\deg(q)=d$. This case
can either be handled by padding (i.e., increment $d\mapsto d+1$ and set
$q_{d+1}=0$) or by including an additive constant:

$r(x)=\frac{\overline{p}(x)}{q(x)}=t+\frac{{p}(x)}{q(x)}=\frac{{p}(x)+tq(x)}{q(x)},\quad\text{
where }\quad
t=\frac{\overline{p}_{d}}{q_{d}}\quad\text{and}\quad{p}(x)=\overline{p}(x)-tq(x),$
(3.12)

so that $\deg({p})\leq d-1$. We then have

$r_{k}=u^{T}W^{k}v+t\mathbb{I}[k=0]$ (3.13)

with $(u,W,v)$ as given in Lemma 3.2 applied to $p(x)/q(x)$. Here $\mathbb{I}$
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 $\mathbb{I}[k=0]=0^{k}$
for $k\geq 0$. Thus we can incorporate the $t\mathbb{I}[k=0]$ into the matrix
representation by appending a row/column of $0$s to $W$, appending $1$ to $u$
and $t$ to $v$:

$\left(\begin{array}[]{c}u\\\
1\end{array}\right)^{T}\left(\begin{array}[]{cc}W&\mathbf{0}\\\
\mathbf{0}^{T}&0\end{array}\right)^{k}\left(\begin{array}[]{c}v\\\
t\end{array}\right)=\left(\begin{array}[]{c}u\\\
1\end{array}\right)^{T}\left(\begin{array}[]{cc}W^{k}&\mathbf{0}\\\
\mathbf{0}^{T}&0^{k}\end{array}\right)\left(\begin{array}[]{c}v\\\
t\end{array}\right)=u^{T}W^{k}v+t\cdot\mathbb{I}[k=0].$ (3.14)

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\mathbb{I}[k=0]$ representation rather than the pure
matrix representation in Eq. 3.14.

###  3.3 Efficient Sampling via BLT Multiplication

The generating function view outlined in Section 3.1 gives us a matrix
factorization from any function $r(x)$. Namely, $A^{(n)}=M(b,n)\cdot M(c,n)$
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)\approx\sqrt{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),n)$. This algorithm relies on the representation in
Section 3.2. The final missing ingredient, which we provide in Sections 4 and
5, is to instantiate $r(x)$.

Our algorithmic task is as follows. Fix a lower-triangular Toeplitz matrix
$M(r,n)\in\mathbb{R}^{n\times n}$ with a generating function $r$. We are given
as input a matrix $Z\in\mathbb{R}^{n\times m}$ and must produce as output
$\widehat{Z}\coloneqq M(r,n)\cdot Z\in\mathbb{R}^{n\times m}$. The obvious
algorithm for this task is to use standard matrix multiplication; this would
take $O(n^{2}m)$ time and $O(nm)$ space. Another approach is to exploit the
Toeplitz structure of $M(r,n)$ which makes the matrix multiplication a
convolution; this means it can be accelerated using the fast Fourier transform
to take $O(nm\log n)$ time, but this is not a streaming algorithm and it still
requires $O(nm)$ 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 $nm$
_independent_ samples from a Gaussian and $\widehat{Z}=M(r,n)\cdot 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 $\widehat{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,\cdot}\in\mathbb{R}^{m}$ and must output the next
row $\widehat{Z}_{k,\cdot}\in\mathbb{R}^{m}$. Our goal is to develop an
algorithm that runs in time and space $O(m)$ per iteration (i.e., $O(nm)$ time
in total). Reducing the space usage relies on the fact that we do not need to
store all of $Z$ or $\widehat{Z}$ in the streaming setting. Accomplishing this
relies on the structure of the lower triangular Toeplitz matrix
$M(r,n)\in\mathbb{R}^{d\times d}$ generated by a rational function

$r(x)=t+\frac{\sum_{i=0}^{d-1}p_{i}x^{i}}{1+\sum_{j=1}^{d}q_{j}x^{j}}=r_{0}+r_{1}x+r_{2}x^{2}+r_{3}x^{3}+\cdots.$

Following 3.3 and Lemma 3.2, we use the representation of the sequence
$(r_{k})_{k=0}^{\infty}$ in terms of matrix powers:
$r_{k}=u^{T}W^{k}v+t\mathbb{I}[k=0]$. This is the representation that we use
for our algorithm.

Parameters: Matrix $M(r,n)$ defined following Lemmas 3.2 and 3.3 via column
vectors $u,v\in\mathbb{R}^{d\times 1}$, a matrix $W\in\mathbb{R}^{d\times d}$,
and a scalar $t\in\mathbb{R}$.

Streaming Input: Row vectors
$Z_{0,\cdot},Z_{1,\cdot},\cdots,Z_{n-1,\cdot}\in\mathbb{R}^{1\times m}$.

Streaming Output: Row vectors
$\widehat{Z}_{0,\cdot},\widehat{Z}_{1,\cdot},\cdots,\widehat{Z}_{n-1,\cdot}\in\mathbb{R}^{1\times
m}$ satisfying $\widehat{Z}_{k}=(\operatorname{M}(r)Z)_{k}$.

Goal: $\widehat{Z}_{k,\cdot}=\sum_{i=0}^{k}r_{i}Z_{k-i,\cdot}$, where
$r_{k}=u^{T}W^{k}v+t\mathbb{I}[k=0]$.

Initialize $S_{0}=0\in\mathbb{R}^{d\times m}$. $\triangleright$ State of the
Algorithm

for $k=0,\cdots,n-1$ do

Receive input $Z_{k,\cdot}\in\mathbb{R}^{1\times m}$.

Compute $S_{k+1}=vZ_{k,\cdot}+W\cdot S_{k}\in\mathbb{R}^{d\times m}$.

Return output
$\widehat{Z}_{k,\cdot}=tZ_{k,\cdot}+u^{T}S_{k+1}\in\mathbb{R}^{1\times m}$.

Algorithm 1 Streaming Multiplication by a BLT Matrix

For a BLT given via a rational generating function $r$ in the matrix
representation of Lemma 3.2, Algorithm 1 in fact computes $\hat{Z}=M(r)Z$ in
row-by-row streaming fashion:

######  Lemma 3.4 (Properties of Algorithm 1).

Algorithm 1 taskes as input a stream
$Z_{0,\cdot},Z_{1,\cdot},\cdots,Z_{n-1,\cdot}\in\mathbb{R}^{1\times m}$ and
outputs a stream
$\widehat{Z}_{0,\cdot},\widehat{Z}_{1,\cdot},\cdots,\widehat{Z}_{n-1,\cdot}\in\mathbb{R}^{1\times
m}$. It takes parameters $u,v\in\mathbb{R}^{d\times 1}$ and
$W\in\mathbb{R}^{d\times d}$ and $t\in\mathbb{R}$. At each iteration $k$,
output satisfies

$\widehat{Z}_{k,\cdot}=\sum_{j=0}^{k}r_{k-j}\cdot
Z_{j,\cdot}=tZ_{k,\cdot}+\sum_{j=0}^{k}u^{T}W^{k-j}v\cdot Z_{j,\cdot},$ (3.15)

where $r_{k}=u^{T}W^{k}v+t\mathbb{I}[k=0]$. The space usage is $O(dm+d^{2})$
and the runtime per iteration is dominated by a matrix multiplication $W\cdot
S$, where $S\in\mathbb{R}^{d\times m}$.

###### Proof.

By induction, for all $k$, we have

$S_{k+1}=\sum_{j=0}^{k}W^{k-j}vZ_{j,\cdot}$

and, hence,

$tZ_{k,\cdot}+u^{T}S_{k+1}=tZ_{k,\cdot}+\sum_{j=0}^{k}u^{T}W^{k-j}vZ_{j,\cdot}=\sum_{j=0}^{k}r_{k-j}Z_{j,\cdot}=\widehat{Z}_{k,\cdot},$

where $\widehat{Z}=M(r,n)Z$.

The state of the algorithm after iteration $k$ is given by
$S_{k+1}\in\mathbb{R}^{d\times m}$. The space usage of the algorithm is
dominated by storing this state ($dm$ 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}=vZ_{k,\cdot}+W\cdot S_{k}$. (The output
computation $\widehat{Z}_{k,\cdot}=tZ_{k,\cdot}+u^{T}S_{k+1}$ requires $m$
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 Eq. 3.9, then the update can be performed
“in place” and we only require $m$ extra registers to store the input
$Z_{k,\cdot}$. In general, we can always perform this update using $dm$
additional registers.

The running time of one iteration consists of the update
$S_{k+1}=vZ_{k,\cdot}+W\cdot S_{k}$ and the matrix-vector product
$\widehat{Z}_{k,\cdot}=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(dm)$ time. In general, we can perform this update in
$O(d^{2}m)$ time. ∎

##  4 Factorizations via Rational Function Approximation

In this section we prove Theorem 1.1. We follow the generating functions view
presented in Section 3 and use the corresponding sampling algorithm (Algorithm
1). The proof is split into three steps. In Section 4.1 we recap how a
generating function $r$ yields a valid matrix factorization and prove that, if
$r(x)\approx\sqrt{1-x}$, then this yields an approximately optimal
factorization. In Section 4.2 we construct a low-degree rational function $r$
that appropriately approximates the square root. Finally we assemble the parts
of the proof in Section 4.3.

###  4.1 Reduction to Approximating the Square Root

Proposition 2.2 shows that the optimal matrix factorization is given by
$A^{(n)}=BC$ where $B=C=M(f,n)$ and $f(x)=1/\sqrt{1-x}$. Suppose we have a
generating function $r(x)\approx\sqrt{1-x}$. Let
$b(x)\coloneqq\frac{r(x)}{1-x}\approx\frac{1}{\sqrt{1-x}}=f(x)$ and
$c(x)\coloneqq\frac{1}{r(x)}\approx\frac{1}{\sqrt{1-x}}$. Since $b(x)\cdot
c(x)=\frac{r(x)}{1-x}\cdot\frac{1}{r(x)}=\frac{1}{1-x}=g(x)$, we always obtain
a valid factorization $M(b,n)\cdot M(c,n)=A^{(n)}$. We can bound the objective
value of this factorization in terms of the approximation
$r(x)\approx\sqrt{1-x}$.

######  Proposition 4.1 (Approximating the square root approximates the matrix
factorization objective).

Let $\tau>0$ and $n\in\mathbb{N}$. Let $r:C\to\mathbb{C}$, where
$C\coloneqq\\{z\in\mathbb{C}:|z|<1\\}$ is the open unit disc in the complex
plane centered at zero. Let

$\gamma_{\tau}\coloneqq\sup\left\\{|r(x)-\sqrt{1-x}|:x\in\mathbb{C},|x|=\exp(-\tau)\right\\}.$

Let $b(x)\coloneqq\frac{r(x)}{1-x}$ and $c(x)\coloneqq\frac{1}{r(x)}$ and
$f(x)=\frac{1}{\sqrt{1-x}}$. Then $M(b,n)\cdot M(c,n)=M(c,n)\cdot
M(b,n)=A^{(n)}=M(f,n)^{2}$, where $M(\cdot,n)$ is defined in Equation 3.1 and
$A^{(n)}$ is defined in Equation 2.1. Moreover, if
$\gamma_{\tau}\leq\left(\frac{1-\exp(-2\tau)}{4}\right)^{2}$, then

$\displaystyle\|M(b,n)\|_{2\to\infty}$
$\displaystyle\leq\|M(f,n)\|_{2\to\infty}+\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt{\exp(2\tau)-1}},$ $\displaystyle\|M(c,n)\|_{1\to
2}$ $\displaystyle\leq\|M(f,n)\|_{1\to 2}+\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt[4]{(\exp(2\tau)-1)^{2}-2^{7/2}\gamma_{\tau}\exp(4\tau)}}.$

We will end up setting $\tau=\frac{1}{2n}$. This ensures that $\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt{\exp(2\tau)-1}}\leq\sqrt{n}\cdot\exp(1/2)\cdot\gamma_{\tau}<2\sqrt{n}\gamma_{\tau}$
and, if $n\geq 6$ and $\gamma_{\tau}\leq\frac{1}{32\cdot n^{2}}$, then
$\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt[4]{(\exp(2\tau)-1)^{2}-2^{7/2}\gamma_{\tau}\exp(4\tau)}}\leq\sqrt[4]{2}\cdot\sqrt{n}\cdot\exp(1/2)\cdot\gamma_{\tau}<2\sqrt{n}\gamma_{\tau}$.

###### Proof.

Let $b(x)=\sum_{k=0}^{\infty}b_{k}x^{k}$,
$c(x)=\sum_{k=0}^{\infty}c_{k}x^{k}$, and
$f(x)=\sum_{k=0}^{\infty}f_{k}x^{k}$. By linearity and the triangle
inequality,

$\displaystyle\|M(b,n)\|_{2\to\infty}$
$\displaystyle=\|M(f,n)+M(b-f,n)\|_{2\to\infty}\leq\|M(f,n)\|_{2\to\infty}+\|M(b-f,n)\|_{2\to\infty},$
$\displaystyle\|M(c,n)\|_{1\to 2}$ $\displaystyle=\|M(f,n)+M(c-f,n)\|_{1\to
2}\leq\|M(f,n)\|_{1\to 2}+\|M(c-f,n)\|_{1\to 2}.$

Thus it suffices to bound
$\|M(b-f,n)\|_{2\to\infty}=\sqrt{\sum_{k=0}^{n-1}(b_{k}-f_{k})^{2}}$ and
$\|M(c-f,n)\|_{1\to 2}=\sqrt{\sum_{k=0}^{n-1}(c_{k}-f_{k})^{2}}$. We bound
these using a weighted version of Parseval’s identity.

######  Lemma 4.2 (Weighted Parseval’s Identity).

Let $g,h:C\to\mathbb{C}$ be analytic, where
$C\coloneqq\\{z\in\mathbb{C}:|z|<1\\}$ is the open unit disc in the complex
plane centered at zero. Let $g(x)=\sum_{k=0}^{\infty}g_{k}x^{k}$ and
$h(x)=\sum_{k=0}^{\infty}h_{k}x^{k}$ for $x\in C$, where
$g_{k},h_{k}\in\mathbb{C}$. Let $\tau>0$. Then

$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}|g(\exp(\sqrt{-1}\theta-\tau))-h(\exp(\sqrt{-1}\theta-\tau))|^{2}\mathrm{d}\theta$
$\displaystyle=\sum_{k=0}^{\infty}|g_{k}-h_{k}|^{2}\cdot\exp(-2\tau k)$
$\displaystyle\geq\exp(-2\tau(n-1))\sum_{k=0}^{n-1}|g_{k}-h_{k}|^{2}.$

###### Proof.

Let $w(x)=g(x)-h(x)=\sum_{k=0}^{\infty}w_{k}x^{k}$, where $w_{k}=g_{k}-h_{k}$.
Then

$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}|g(\exp(\sqrt{-1}\theta-\tau))-h(\exp(\sqrt{-1}\theta-\tau))|^{2}\mathrm{d}\theta$
$\displaystyle=\frac{1}{2\pi}\int_{-\pi}^{\pi}|w(\exp(\sqrt{-1}\theta-\tau))|^{2}\mathrm{d}\theta$
$\displaystyle=\frac{1}{2\pi}\int_{-\pi}^{\pi}w(\exp(\sqrt{-1}\theta-\tau))\cdot\overline{w(\exp(\sqrt{-1}\theta-\tau))}\mathrm{d}\theta$
$\displaystyle=\frac{1}{2\pi}\int_{-\pi}^{\pi}\left(\sum_{k=0}^{\infty}w_{k}\cdot(\exp(\sqrt{-1}\theta-\tau))^{k}\right)\cdot\overline{\left(\sum_{\ell=0}^{\infty}w_{\ell}\cdot(\exp(\sqrt{-1}\theta-\tau))^{\ell}\right)}\mathrm{d}\theta$
$\displaystyle=\sum_{k,\ell=0}^{\infty}\frac{1}{2\pi}\int_{-\pi}^{\pi}w_{k}\cdot(\exp(\sqrt{-1}\theta-\tau))^{k}\cdot\overline{w_{\ell}}\cdot(\exp(-\sqrt{-1}\theta-\tau))^{\ell}\mathrm{d}\theta$
$\displaystyle=\sum_{k,\ell=0}^{\infty}w_{k}\cdot\overline{w_{\ell}}\cdot\exp(-(k+\ell)\tau)\cdot\frac{1}{2\pi}\int_{-\pi}^{\pi}\exp((k-\ell)\sqrt{-1}\theta)\mathrm{d}\theta$
$\displaystyle=\sum_{k,\ell=0}^{\infty}w_{k}\cdot\overline{w_{\ell}}\cdot\exp(-(k+\ell)\tau)\cdot\mathbb{I}[k-\ell=0]$
$\displaystyle=\sum_{k=0}^{\infty}|w_{k}|^{2}\cdot\exp(-2k\tau)$
$\displaystyle=\sum_{k=0}^{\infty}|g_{k}-h_{k}|^{2}\cdot\exp(-2\tau k).$

∎

Recall $b(x)\coloneqq\frac{r(x)}{1-x}$ and $f(x)=\frac{1}{\sqrt{1-x}}$, whence
$b(x)-f(x)=\frac{r(x)-\sqrt{1-x}}{1-x}$. Now we have

$\displaystyle\|M(b-f,n)\|_{2\to\infty}$
$\displaystyle=\sqrt{\sum_{k=0}^{n-1}(b_{k}-f_{k})^{2}}$
$\displaystyle\leq\exp(\tau(n-1))\sqrt{\frac{1}{2\pi}\int_{-\pi}^{\pi}|b(\exp(\sqrt{-1}\theta-\tau))-f(\exp(\sqrt{-1}\theta-\tau))|^{2}\mathrm{d}\theta}$
$\displaystyle=\exp(\tau(n-1))\sqrt{\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{\left|r(\exp(\sqrt{-1}\theta-\tau))-\sqrt{1-\exp(\sqrt{-1}\theta-\tau))}\right|^{2}}{\left|1-\exp(\sqrt{-1}\theta-\tau))\right|^{2}}\mathrm{d}\theta}$
$\displaystyle\leq\exp(\tau(n-1))\sqrt{\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{\gamma_{\tau}^{2}}{\left|1-\exp(\sqrt{-1}\theta-\tau))\right|^{2}}\mathrm{d}\theta}$
$\displaystyle=\exp(\tau(n-1))\cdot\gamma_{\tau}\cdot\sqrt{\frac{\exp(\tau)}{4\pi}\int_{-\pi}^{\pi}\frac{1}{\cosh(\tau)+\cos(\theta)}\mathrm{d}\theta}$
$\displaystyle=\exp(\tau(n-1))\cdot\gamma_{\tau}\cdot\sqrt{\frac{\exp(\tau)}{2}\frac{1}{\sqrt{\cosh(\tau)^{2}-1}}}$
$\displaystyle=\frac{\exp(\tau(n-1))\cdot\gamma_{\tau}}{\sqrt{1-\exp(-2\tau)}}=\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt{\exp(2\tau)-1}}\leq\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt{2\tau}}.$

In the above we use the identity

$\displaystyle|1-\exp(\sqrt{-1}\theta-\tau)|^{2}$
$\displaystyle=(1-e^{-\tau}\cos\theta)^{2}+(e^{-\tau}\sin\theta)^{2}=1-2e^{-\tau}\cos\theta+e^{-2\tau}(\cos^{2}\theta+\sin^{2}\theta)$
$\displaystyle=e^{-\tau}\cdot\left(e^{\tau}-2\cos\theta+e^{-\tau}\right)=2e^{-\tau}\cdot\left(\cosh\tau+\cos\theta\right)$

and the integral

$\forall\sigma>1~{}~{}~{}~{}~{}\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{1}{\sigma+\cos(\theta)}\mathrm{d}\theta=\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{1}{\sigma+\sin(\theta)}\mathrm{d}\theta=\frac{1}{\sqrt{\sigma^{2}-1}},$
(4.1)

which follows from the derivative
$\frac{\mathrm{d}}{\mathrm{d}\theta}\left[\frac{2}{\sqrt{\sigma^{2}-1}}\tan^{-1}\left(\frac{1+\sigma\tan(\theta/2)}{\sqrt{\sigma^{2}-1}}\right)\right]=\frac{1}{\sigma+\sin(\theta)}$.

Similarly, recall $c(x)\coloneqq\frac{1}{r(x)}$ and
$f(x)=\frac{1}{\sqrt{1-x}}$, whence
$c(x)-f(x)=\frac{\sqrt{1-x}-r(x)}{r(x)\cdot\sqrt{1-x}}$. Assuming
$\gamma_{\tau}<\frac{(1-\exp(-2\tau))^{2}}{2^{7/2}}$, we have

$\displaystyle\|M(c-f,n)\|_{1\to 2}$
$\displaystyle=\sqrt{\sum_{k=0}^{n-1}(c_{k}-f_{k})^{2}}$
$\displaystyle\leq\exp(\tau(n-1))\sqrt{\frac{1}{2\pi}\int_{-\pi}^{\pi}|c(\exp(\sqrt{-1}\theta-\tau))-f(\exp(\sqrt{-1}\theta-\tau))|^{2}\mathrm{d}\theta}$
$\displaystyle=\exp(\tau(n\\!-\\!1))\\!\\!\sqrt{\\!\frac{1}{2\pi}\\!\int_{-\pi}^{\pi}\\!\frac{\left|r(\exp(\sqrt{-1}\theta-\tau))-\sqrt{1-\exp(\sqrt{-1}\theta-\tau)}\right|^{2}}{\left|r(\exp(\sqrt{-1}\theta\\!-\\!\tau))\\!\cdot\\!\sqrt{1\\!-\\!\exp(\sqrt{-1}\theta\\!-\\!\tau)}\right|^{2}}\\!\mathrm{d}\theta}$
$\displaystyle=\exp(\tau(n\\!-\\!1))\\!\\!\sqrt{\\!\frac{1}{2\pi}\\!\int_{-\pi}^{\pi}\\!\frac{\left|r(\exp(\sqrt{-1}\theta-\tau))-\sqrt{1-\exp(\sqrt{-1}\theta-\tau)}\right|^{2}}{\left|1\\!-\\!\exp(\\!\sqrt{-1}\theta\\!-\\!\tau\\!)\\!+\\!\left(r(\exp(\sqrt{-1}\theta\\!-\\!\tau))\\!-\\!\sqrt{1\\!-\\!\exp(\sqrt{-1}\theta\\!-\\!\tau)}\right)\\!\cdot\\!\sqrt{1\\!-\\!\exp(\sqrt{-1}\theta\\!-\\!\tau)}\right|^{2}}\\!\mathrm{d}\theta}$
$\displaystyle\leq\exp(\tau(n\\!-\\!1))\\!\\!\sqrt{\\!\frac{1}{2\pi}\\!\int_{-\pi}^{\pi}\\!\frac{\gamma_{\tau}^{2}}{\left(\left|1\\!-\\!\exp(\\!\sqrt{-1}\theta\\!-\\!\tau\\!)\right|-\gamma_{\tau}\cdot\sqrt{|1-\exp(\sqrt{-1}\theta-\tau)|}\right)^{2}}\mathrm{d}\theta}$
$\displaystyle\leq\exp(\tau(n\\!-\\!1))\\!\\!\sqrt{\\!\frac{1}{2\pi}\\!\int_{-\pi}^{\pi}\\!\frac{\gamma_{\tau}^{2}}{\left|1\\!-\\!\exp(\\!\sqrt{-1}\theta\\!-\\!\tau\\!)\right|^{2}-2\sqrt{2}\gamma_{\tau}}\mathrm{d}\theta}$
$\displaystyle=\exp(\tau(n\\!-\\!1))\\!\cdot\\!\gamma_{\tau}\\!\cdot\\!\sqrt{\\!\frac{\exp(\tau)}{4\pi}\\!\int_{-\pi}^{\pi}\\!\frac{1}{\cosh(\tau)-\sqrt{2}\exp(\tau)\gamma_{\tau}+\cos(\theta)}\mathrm{d}\theta}$
$\displaystyle=\exp(\tau(n\\!-\\!1))\\!\cdot\\!\gamma_{\tau}\\!\cdot\\!\sqrt{\\!\frac{\exp(\tau)}{2}\frac{1}{\sqrt{(\cosh(\tau)-\sqrt{2}\exp(\tau)\gamma_{\tau})^{2}-1}}}$
$\displaystyle=\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt{2\exp(\tau)\sqrt{(\cosh(\tau)-\sqrt{2}\exp(\tau)\gamma_{\tau})^{2}-1}}}$
$\displaystyle=\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt[4]{{(\exp(2\tau)(1-2\sqrt{2}\gamma_{\tau})+1)^{2}-4\exp(2\tau)}}}$
$\displaystyle=\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt[4]{1-2(1+2^{3/2}\gamma_{\tau})\exp(2\tau)+(1-2^{3/2}\gamma_{\tau})^{2}\exp(4\tau)}}$
$\displaystyle=\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt[4]{(\exp(2\tau)-1)^{2}-2^{5/2}\gamma_{\tau}\exp(2\tau)-2^{5/2}\gamma_{\tau}\exp(4\tau)+2^{3}\gamma_{\tau}^{2}\exp(4\tau)}}$
$\displaystyle\leq\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt[4]{(\exp(2\tau)-1)^{2}-2^{5/2}\gamma_{\tau}\exp(2\tau)(1+\exp(2\tau))}}$
$\displaystyle\leq\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt[4]{(\exp(2\tau)-1)^{2}-2^{7/2}\gamma_{\tau}\exp(4\tau)}}.$

Combining the bounds yields the result. ∎

###  4.2 Rational Approximation of $\sqrt{x}$

Sections 3 and 4.1 reduce the problem to finding a low-degree rational
function that uniformly approximates $\sqrt{1-x}$ for $x\in\mathbb{C}$ with
$|x|\leq 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 $\sqrt{x}$ instead of $\sqrt{1-x}$. To get an
equivalent result, we must uniformly approximate $\sqrt{x}$ for
$x\in\mathbb{C}$ with $|x-1/2|\leq 1/2$.

[New64] gave an explicit rational function that approximates $\sqrt{x}$ for
$x\in[0,1]$. Specifically, let

$p(x)=\prod_{k=0}^{d-1}(x+\exp(-k/\sqrt{d}))~{}~{}~{}\text{ and
}~{}~{}~{}r(x)=\frac{\sqrt{x}\cdot(p(\sqrt{x})-p(-\sqrt{x}))}{p(\sqrt{x})+p(-\sqrt{x})}.$
(4.2)

Then $r(x)$ is a rational function of degree $\lceil d/2\rceil$ –
specifically, $r(x)$ is the ratio where the numerator is a polynomial in $x$
of degree $\lceil d/2\rceil$ and the denominator is a polynomial in $x$ of
degree $\lfloor d/2\rfloor$. Furthermore,

$\sup_{x\in[0,1]}|r(x)-\sqrt{x}|\leq 3\cdot\exp(-\sqrt{d}).$ (4.3)

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 [Sta93] obtained optimal
constants for this approximation.

In order to apply Proposition 4.1, we need to extend Newman’s result into the
complex plane. We provide a proof adapted from that of [GT19].

######  Theorem 4.3.

Let $d\in\mathbb{N}$ with $d\geq 2$. Then there exists a rational function $r$
of degree $d$ with real coefficients and real negative simple poles such that

$|r(x)-\sqrt{x}|\leq\big{(}4+\frac{4}{\pi}\big{)}\cdot\exp\left(-\pi\sqrt{\frac{1}{2}\left\lfloor\frac{d-1}{2}\right\rfloor}\right)\leq
6\cdot\exp\left(-\frac{\pi}{2}\sqrt{d-2}\right)$ (4.4)

for all $x\in\mathbb{C}$ with $\Re(x)\geq 0$ and $|x|\leq 1$. Specifically,

$r(x)=\frac{2h}{\pi}\sum_{k=-d_{-}}^{d_{+}}\frac{x\cdot\exp(hk)}{x+\exp(2hk)},~{}~{}\text{
where
}~{}~{}d_{+}=\left\lfloor\frac{d-1}{2}\right\rfloor,~{}~{}d_{-}=\left\lceil\frac{d-1}{2}\right\rceil,~{}~{}\text{
and }~{}~{}h=\frac{\pi}{\sqrt{2d_{+}}}.$ (4.5)

###### Proof.

We invoke Proposition 4.5. Let $d_{-}=\lceil(d-1)/2\rceil$ and
$d_{+}=\lfloor(d-1)/2\rfloor$ and let $h\in(0,4.9)$ be determined later. Let
$r(x)=r_{d_{+},d_{-},h}(x)$ be as in Equation 4.11. Clearly $r(x)$ is a
rational function of degree $d_{+}+d_{-}+1=d$ with real coefficients and the
poles are given by $-\exp(2hk)\in(-\infty,0)$ for
$k\in\\{-d_{-},-d_{-}+1,\cdots,d_{+}\\}$. It remains to simplify the
approximation guarantee and set $h$. From Equation 4.11, for all
$x\in\mathbb{C}$ with $\Re(x)\geq 0$ and $|x|\leq 1$, we have

$\displaystyle|r(x)-\sqrt{x}|$ $\displaystyle\leq
2\sqrt{|x|}\left(\frac{1}{\exp((1-c)\pi^{2}/h)-1}+\frac{1}{\exp((1+c)\pi^{2}/h)-1}\right)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{2h}{\pi(\exp(h)-1)}\big{(}|x|\cdot{\exp(-hd_{+})}+{\exp(-hd_{-})}\big{)}$
$\displaystyle\leq
2\left(\frac{1}{\exp((1-|c|)\pi^{2}/h)-1}+\frac{1}{\exp((1+|c|)\pi^{2}/h)-1}\right)$
($|x|\leq 1$ & symmetry in $c$)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{4h\exp(-hd_{+})}{\pi(\exp(h)-1)}$
($d_{-}\geq d_{+}$) $\displaystyle\leq
2\left(\frac{1}{\exp(\pi^{2}/2h)-1}+\frac{1}{\exp(\pi^{2}/h)-1}\right)+\frac{4}{\pi}\cdot\exp(-hd_{+})\cdot\frac{h}{\exp(h)-1}$
($|c|=|\frac{1}{\pi}\arg(x)|\leq\frac{1}{2}$) $\displaystyle\leq
4\cdot\exp(-\pi^{2}/2h)+\frac{4}{\pi}\cdot\exp(-hd_{+})$ (4.6)
$\displaystyle=\big{(}4+\frac{4}{\pi}\big{)}\cdot\exp\left(-\pi\sqrt{\frac{d_{+}}{2}}\right)$
(4.7)
$\displaystyle=\big{(}4+\frac{4}{\pi}\big{)}\cdot\exp\left(-\pi\sqrt{\frac{1}{2}\left\lfloor\frac{d-1}{2}\right\rfloor}\right)$
$\displaystyle\leq 6\cdot\exp\left(-\frac{\pi}{2}\sqrt{d-1}\right).$

The fourth inequality Eq. 4.6 follows from two simple facts: First, $h\geq 0$
implies $\frac{h}{\exp(h)-1}\leq 1$.Second,
$\frac{1}{a-1}+\frac{1}{a^{2}-1}\leq\frac{2}{a}$ when $a=\exp(\pi^{2}/2h)\geq
1+\sqrt{3}$ or, equivalently, $h\leq\frac{\pi^{2}}{2\log(1+\sqrt{3})}\approx
4.91$. The penultimate equality Eq. 4.7 follows by setting
$h=\frac{\pi}{\sqrt{2d_{+}}}\leq\frac{\pi}{\sqrt{2}}<2.5$. ∎

######  Corollary 4.4 (Rescaling of Theorem 4.3).

Let $d\geq 2$ be an integer. Then there exists a rational function
$\widetilde{r}$ of degree $d$ with real coefficients such that all poles are
simple and $>1$ and we have

$|\widetilde{r}(x)-\sqrt{1-x}|\leq
8\cdot\exp\left(-\frac{\pi}{2}\sqrt{d-2}\right)$ (4.8)

for all $x\in\mathbb{C}$ with $|x|\leq 1$. Specifically,

$\widetilde{r}(x)=\frac{2h\sqrt{2}}{\pi}\sum_{k=-d_{-}}^{d_{+}}\frac{(1-x)\cdot\exp(hk)}{1-x+2\cdot\exp(2hk)}=\frac{2h\sqrt{2}}{\pi}\sum_{k=-d_{-}}^{d_{+}}\exp(hk)-\frac{2\cdot\exp(3hk)}{1+2\cdot\exp(2hk)-x}$
(4.9)

where $d_{+}=\lfloor(d-1)/2\rfloor$, $d_{-}=\lceil(d-1)/2\rceil$, and
$h=\frac{\pi}{\sqrt{2d_{+}}}$.

###### Proof.

Let $r$ be the rational function of degree $\leq d$ promised by Theorem 4.3.
Define $\widetilde{r}(x)=\sqrt{2}\cdot r(\tfrac{1-x}{2})$. Suppose
$x\in\mathbb{C}$ with $|x|\leq 1$. Then
$\Re(\tfrac{1-x}{2})=\frac{1}{2}-\frac{1}{2}\Re(x)\geq 0$,
$|\tfrac{1-x}{2}|\leq\tfrac{1+|x|}{2}\leq 1$, and

$|\widetilde{r}(x)-\sqrt{1-x}|=\sqrt{2}\cdot\left|r(\tfrac{1-x}{2})-\sqrt{\tfrac{1-x}{2}}\right|\leq\sqrt{2}\cdot\left(4+\frac{4}{\pi}\right)\cdot\exp\left(-\pi\sqrt{\frac{1}{2}\left\lceil\frac{d-1}{2}\right\rceil}\right)\leq
8\exp\left(-\frac{\pi}{2}\sqrt{d-1}\right).$

If $\widetilde{x}$ is a pole of $\widetilde{r}$, then
$\tfrac{1-\widetilde{x}}{2}$ is a pole of $r$; since the poles of $r$ are
negative we have $\tfrac{1-\widetilde{x}}{2}<0$ and, hence, $\widetilde{x}>1$.
∎

The general form of the approximation promised by Theorem 4.3 is given below.

######  Proposition 4.5.

Let $d_{+},d_{-}\in\mathbb{N}$ and $h>0$. For
$x\in\mathbb{C}\setminus[-\exp(2d_{+}h),-\exp(-2d_{-}h)]$, define

$r_{d_{+},d_{-},h}(x)\coloneqq\frac{2h}{\pi}\sum_{k=-d_{-}}^{d_{+}}\frac{x\cdot\exp(hk)}{x+\exp(2hk)}.$
(4.10)

Then, for all $x\in\mathbb{C}$ with $\Re(x)\geq 0$, we have

$\displaystyle|r_{d_{+},d_{-},h}(x)-\sqrt{x}|$ $\displaystyle\leq
2\sqrt{|x|}\left(\frac{1}{\exp((1-c)\pi^{2}/h)-1}+\frac{1}{\exp((1+c)\pi^{2}/h)-1}\right)$
$\displaystyle~{}~{}~{}~{}~{}+\frac{2h}{\pi(\exp(h)-1)}\big{(}|x|\cdot{\exp(-hd_{+})}+{\exp(-hd_{-})}\big{)},$
(4.11)

where $c=\frac{1}{\pi}\arg(x)\in\left[-\frac{1}{2},\frac{1}{2}\right]$.

###### Proof.

Fix $x\in\mathbb{C}$ with $\Re(x)>0$.222222For the proof we assume $\Re(x)>0$.
The case $\Re(x)=0$ follows by continuity. Let
$\arg(x)=\tan^{-1}\left(\frac{\Im(x)}{\Re(x)}\right)=c\pi\in\left(-\frac{\pi}{2},\frac{\pi}{2}\right)$.

We begin with the Cauchy distribution integral

$\int_{0}^{\infty}\frac{1}{1+u^{2}}\mathrm{d}u=\int_{0}^{\infty}\left(\frac{\mathrm{d}}{\mathrm{d}u}\tan^{-1}(u)\right)\mathrm{d}u=\frac{\pi}{2}.$

We perform the variable substitutions $u=v/\sqrt{x}$ and $v=\exp(s)$ and
rearrange to obtain

$\sqrt{x}=\frac{2}{\pi}\int_{0}^{\infty}\frac{x}{x+v^{2}}\mathrm{d}v=\frac{2}{\pi}\int_{-\infty}^{\infty}\frac{x\cdot\exp(s)}{x+\exp(2s)}\mathrm{d}s.$
(4.12)

Now we make the approximation

$\sqrt{x}=\frac{2}{\pi}\int_{-\infty}^{\infty}\frac{x\cdot\exp(s)}{x+\exp(2s)}\mathrm{d}s\approx\frac{2h}{\pi}\sum_{k\in\mathbb{Z}}\frac{x\cdot\exp(hk)}{x+\exp(2hk)}\approx\frac{2h}{\pi}\sum_{k=-d_{-}}^{d_{+}}\frac{x\cdot\exp(hk)}{x+\exp(2hk)}=r_{d_{+},d_{-},h}(x).$
(4.13)

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)\coloneqq\frac{2h}{\pi}\frac{x\cdot\exp(hz)}{x+\exp(2hz)}$. Let

$\displaystyle\widehat{f}(t)$
$\displaystyle\coloneqq\int_{-\infty}^{\infty}f(z)\cdot\exp(-2\pi
zt\sqrt{-1})\mathrm{d}z$
$\displaystyle=\int_{-\infty}^{\infty}\frac{2h}{\pi}\frac{x\cdot\exp(hz-2\pi
zt\sqrt{-1})}{x+\exp(2hz)}\mathrm{d}z$
$\displaystyle=\frac{2h}{\pi}\int_{-\infty}^{\infty}\frac{\exp(-2\pi
zt\sqrt{-1})}{\exp(-hz)+\exp(hz-\log x)}\mathrm{d}z$
$\displaystyle=\frac{2h}{\pi}\int_{-\infty}^{\infty}\frac{\exp(-2\pi(y+\frac{\log
x}{2h})t\sqrt{-1})}{\exp(-hy-\frac{1}{2}\log x)+\exp(hy-\frac{1}{2}\log
x)}\mathrm{d}y$ ($z=y+\frac{\log x}{2h}$)
$\displaystyle=\sqrt{x}\cdot\frac{2h}{\pi}\cdot\exp\left(-\pi t\frac{\log
x}{h}\sqrt{-1}\right)\cdot\int_{-\infty}^{\infty}\frac{\exp(-2\pi
yt\sqrt{-1})}{\exp(-hy)+\exp(hy)}\mathrm{d}y$
$\displaystyle=\sqrt{x}\cdot\frac{2h}{\pi}\cdot\exp\left(-\pi t\frac{\log
x}{h}\sqrt{-1}\right)\cdot\frac{1}{2}\int_{-\infty}^{\infty}\frac{\exp(-2\pi
yt\sqrt{-1})}{\cosh(hy)}\mathrm{d}y$
$\displaystyle=\sqrt{x}\cdot\frac{2h}{\pi}\cdot\exp\left(-\pi t\frac{\log
x}{h}\sqrt{-1}\right)\cdot\frac{\pi}{2h\cdot\cosh(\pi^{2}t/h)}$ (4.14)
$\displaystyle=\sqrt{x}\cdot\exp\left(-\pi t\frac{\log
x}{h}\sqrt{-1}\right)\cdot\frac{1}{\cosh(\pi^{2}t/h)}.$

The penultimate equality (4.14) is the Fourier transform of the hyperbolic
secant [use14]. The Poisson summation formula states that

$\sum_{k\in\mathbb{Z}}f(k)=\sum_{t\in\mathbb{Z}}\widehat{f}(t).$ (4.15)

Per Equation 4.12, we have

$\widehat{f}(0)=\int_{-\infty}^{\infty}f(z)\mathrm{d}z=\int_{-\infty}^{\infty}\frac{2h}{\pi}\frac{x\cdot\exp(hz)}{x+\exp(2hz)}\mathrm{d}z=\int_{-\infty}^{\infty}\frac{2}{\pi}\frac{x\cdot\exp(s)}{x+\exp(2s)}\mathrm{d}s=\sqrt{x}.$

Thus the Poisson summation formula allows us to analyze the first
approximation in Equation 4.13. Namely,

$\displaystyle\left|\sqrt{x}-\frac{2h}{\pi}\sum_{k\in\mathbb{Z}}\frac{x\cdot\exp(hk)}{x+\exp(2hk)}\right|$
$\displaystyle=\left|\widehat{f}(0)-\sum_{k\in\mathbb{Z}}f(k)\right|$
$\displaystyle=\left|\widehat{f}(0)-\sum_{t\in\mathbb{Z}}\widehat{f}(t)\right|$
$\displaystyle\leq\sum_{t=1}^{\infty}|\widehat{f}(t)|+|\widehat{f}(-t)|$
$\displaystyle=\sum_{t=1}^{\infty}|\sqrt{x}|\cdot\left(\left|\exp\left(-\pi
t\frac{\log x}{h}\sqrt{-1}\right)\right|+\left|\exp\left(\pi t\frac{\log
x}{h}\sqrt{-1}\right)\right|\right)\cdot\frac{1}{\cosh(\pi^{2}t/h)}$
$\displaystyle=\sum_{t=1}^{\infty}|\sqrt{x}|\cdot\left(\exp\left(\Im\left(\pi
t\frac{\log x}{h}\right)\right)+\exp\left(\Im\left(-\pi t\frac{\log
x}{h}\right)\right)\right)\cdot\frac{1}{\cosh(\pi^{2}t/h)}$
($|\exp(v\sqrt{-1})|=\exp(\Re(v\sqrt{-1}))=\exp(\Im(-v))$)
$\displaystyle=\sqrt{|x|}\sum_{t=1}^{\infty}2\cosh\left(\pi\frac{t}{h}\Im(\log
x)\right)\cdot\frac{1}{\cosh(\pi^{2}t/h)}$
$\displaystyle=\sqrt{|x|}\sum_{t=1}^{\infty}2\cosh\left(\pi\frac{t}{h}\cdot
c\pi\right)\cdot\frac{1}{\cosh(\pi^{2}t/h)}$ ($\Im(\log x)=\arg(x)=c\pi$)
$\displaystyle=2\sqrt{|x|}\sum_{t=1}^{\infty}\frac{\exp(c\pi^{2}t/h)+\exp(-c\pi^{2}t/h)}{\exp(\pi^{2}t/h)+\exp(-\pi^{2}t/h)}$
$\displaystyle\leq
2\sqrt{|x|}\sum_{t=1}^{\infty}\frac{\exp(c\pi^{2}t/h)+\exp(-c\pi^{2}t/h)}{\exp(\pi^{2}t/h)+0}$
$\displaystyle=2\sqrt{|x|}\sum_{t=1}^{\infty}\exp(-(1-c)\pi^{2}t/h)+\exp(-(1+c)\pi^{2}t/h)$
$\displaystyle=2\sqrt{|x|}\left(\frac{1}{\exp((1-c)\pi^{2}/h)-1}+\frac{1}{\exp((1+c)\pi^{2}/h)-1}\right).$
(4.16)

The final equality (4.16) follows from the usual geometric series formula: If
$a>0$, then
$\sum_{t=1}^{\infty}\exp(-at)=\frac{\exp(-a)}{1-\exp(-a)}=\frac{1}{\exp(a)-1}$.

Next we analyze the second approximation in Equation 4.13, which requires
bounding the tails of $f$. On one side, for $k\in\mathbb{N}$, we have

$\displaystyle|f(k)|$
$\displaystyle=\left|\frac{2h}{\pi}\frac{x\cdot\exp(hk)}{x+\exp(2hk)}\right|$
$\displaystyle=\frac{2h}{\pi}\frac{|x|}{|x\cdot\exp(-hk)+\exp(hk)|}$
$\displaystyle=\frac{2h}{\pi}\frac{|x|}{\sqrt{\Re(x\cdot\exp(-hk)+\exp(hk))^{2}+\Im(x\cdot\exp(-hk)+\exp(hk))^{2}}}$
$\displaystyle\leq\frac{2h}{\pi}\frac{|x|}{|\Re(x\cdot\exp(-hk)+\exp(hk)|}$
$\displaystyle=\frac{2h}{\pi}\frac{|x|}{|\Re(x)\cdot\exp(-hk)+\exp(hk)|}$
($hk\in\mathbb{R}$)
$\displaystyle=\frac{2h}{\pi}\frac{|x|}{\Re(x)\cdot\exp(-hk)+\exp(hk)}$
($\Re(x)>0$) $\displaystyle\leq\frac{2h}{\pi}\frac{|x|}{\exp(hk)}.$ (4.17)

On the other side, for $k\in\mathbb{N}$, we have

$\displaystyle|f(-k)|$
$\displaystyle=\left|\frac{2h}{\pi}\frac{x\cdot\exp(-hk)}{x+\exp(-2hk)}\right|$
$\displaystyle=\frac{2h}{\pi}\frac{\exp(-hk)}{|1+\exp(-2hk)/x|}$
$\displaystyle\leq\frac{2h}{\pi}\frac{\exp(-hk)}{|\Re(1+\exp(-2hk)/x)|}$
$\displaystyle=\frac{2h}{\pi}\frac{\exp(-hk)}{|1+\exp(-2hk)/\Re(x)|}$
($\Re(1/x)=1/\Re(x)$)
$\displaystyle=\frac{2h}{\pi}\frac{\exp(-hk)}{1+\exp(-2hk)/\Re(x)}$
($\Re(x)>0$) $\displaystyle\leq\frac{2h}{\pi}\exp(-hk).$ (4.18)

Combining Equations 4.17 and 4.18, we have

$\displaystyle\left|\frac{2h}{\pi}\sum_{k\in\mathbb{Z}}\frac{x\cdot\exp(hk)}{x+\exp(2hk)}-r_{d_{+},d_{-},h}(x)\right|$
$\displaystyle=\frac{2h}{\pi}\left|\sum_{k\in\mathbb{Z}}\frac{x\cdot\exp(hk)}{x+\exp(2hk)}-\sum_{k=-d_{-}}^{d_{+}}\frac{x\cdot\exp(hk)}{x+\exp(2hk)}\right|$
$\displaystyle=\left|\sum_{k=d_{+}+1}^{\infty}f(k)+\sum_{k=d_{-}+1}^{\infty}f(-k)\right|$
$\displaystyle\leq\sum_{k=d_{+}+1}^{\infty}|f(k)|+\sum_{k=d_{-}+1}^{\infty}|f(-k)|$
$\displaystyle\leq\frac{2h}{\pi}\sum_{k=d_{+}+1}^{\infty}\frac{|x|}{\exp(hk)}+\frac{2h}{\pi}\sum_{k=d_{-}+1}^{\infty}\exp(-hk)$
$\displaystyle=\frac{2h}{\pi}|x|\frac{\exp(-h(d_{+}+1))}{1-\exp(-h)}+\frac{2h}{\pi}\frac{\exp(-h(d_{-}+1))}{1-\exp(-h)}$
(4.19)
$\displaystyle=\frac{2h}{\pi(\exp(h)-1)}\big{(}|x|\cdot{\exp(-hd_{+})}+{\exp(-hd_{-})}\big{)}.$
(4.20)

Combining Equations 4.16 and 4.20 yields the result in Equation 4.11. ∎

###  4.3 Putting Everything Together

Now we can assemble the tools we have developed to prove our main result
(Theorem 1.1). In Section 4.1 we connected near-optimal lower triangular
Toeplitz matrix factorizations to approximating the square root function. In
Sections 3.1 and 3.3 we connected low-degree rational functions to efficient
streaming algorithms. In Section 4.2 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.

######  Theorem 4.6 (Main Theorem - Formal version of Theorem 1.1).

Let $n,d,m\in\mathbb{N}$ satisfy $n\geq 5$ and $d\geq 2+\left(\frac{12+4\log
n}{\pi}\right)^{2}$. There exists a lower triangular Toeplitz matrix
factorization $B,C,A=BC\in\mathbb{R}^{n\times n}$ with the following
properties.

  * •

Near-Optimality: Let $B^{*}=C^{*}=M(f,n)$ for $f(x)=1/\sqrt{1-x}$ be the
optimal lower triangular Toeplitz matrix factorization of size $n$. Then

$\displaystyle\|B\|_{2\to\infty}$
$\displaystyle\leq\|B^{*}\|_{2\to\infty}+16\sqrt{n}\exp\left(-\frac{\pi}{2}\sqrt{d-2}\right),$
$\displaystyle\|C\|_{1\to 2}$ $\displaystyle\leq\|C^{*}\|_{1\to
2}+16\sqrt{n}\exp\left(-\frac{\pi}{2}\sqrt{d-2}\right).$

  * •

Efficiency: There is an algorithm that takes as input a stream
$Z_{1,\cdot},Z_{2,\cdot},\cdots,Z_{n,\cdot}\in\mathbb{R}^{1\times m}$ and
outputs a stream
$\widehat{Z}_{1,\cdot},\widehat{Z}_{2,\cdot},\cdots,\widehat{Z}_{n,\cdot}\in\mathbb{R}^{1\times
m}$ such that $\widehat{Z}=C^{-1}Z$ and the space and time (per iteration)
required by the algorithm is $O(dm)$.

To obtain Theorem 1.1 in the introduction from Theorem 4.6, we simply need to
set $d=O(\log(n/\mu))^{2}$ such that $d\geq 2+\left(\frac{12+4\log
n}{\pi}\right)^{2}$ and
$16\sqrt{n}\exp\left(-\frac{\pi}{2}\sqrt{d-2}\right)\leq\frac{\mu}{4+\log(n)}$.

###### Proof.

Corollary 4.4 gives us a rational function $r$ of degree $d$ with real
coefficients such that

$|{r}(x)-\sqrt{1-x}|\leq 8\cdot\exp\left(-\frac{\pi}{2}\sqrt{d-2}\right)$

for all $x\in\mathbb{C}$ with $|x|\leq 1$.

Let $b(x)\coloneqq\frac{r(x)}{1-x}$ and $c(x)\coloneqq\frac{1}{r(x)}$ and
$f(x)\coloneqq\frac{1}{\sqrt{1-x}}$ and $g(x)\coloneqq 1/(1-x)=b(x)c(x)$. Let
$B\coloneqq M(b,n),C\coloneqq M(c,n)\in\mathbb{R}^{n\times n}$ be as in
Equation 3.1. By Lemma 3.1, this is a valid matrix factorization, as $B\cdot
C=M(b,n)\cdot M(c,n)=M(g,n)$.

Let $\tau>0$ and

$\gamma_{\tau}\coloneqq\sup\left\\{|r(x)-\sqrt{1-x}|:x\in\mathbb{C},|x|=\exp(-\tau)\right\\}\leq
8\cdot\exp\left(-\frac{\pi}{2}\sqrt{d-2}\right).$

By Proposition 4.1, if
$\gamma_{\tau}\leq\left(\frac{1-\exp(-2\tau)}{4}\right)^{2}$, then

$\displaystyle\|M(b,n)\|_{2\to\infty}$
$\displaystyle\leq\|M(f,n)\|_{2\to\infty}+\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt{\exp(2\tau)-1}},$ $\displaystyle\|M(c,n)\|_{1\to
2}$ $\displaystyle\leq\|M(f,n)\|_{1\to 2}+\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt[4]{(\exp(2\tau)-1)^{2}-2^{7/2}\gamma_{\tau}\exp(4\tau)}}.$

Set $\tau=\frac{1}{2n}$. Then $\exp(2\tau)-1\geq 1/n$ and $\exp(\tau
n)=\exp(1/2)$. Thus

$\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt{\exp(2\tau)-1}}\leq\sqrt{n}\cdot\exp(1/2)\cdot\gamma_{\tau}<2\sqrt{n}\gamma_{\tau}$

and, if $n\geq 5$ and $\gamma_{\tau}\leq\left(\frac{1}{6\cdot n}\right)^{2}$,
then

$\frac{\exp(\tau
n)\cdot\gamma_{\tau}}{\sqrt[4]{(\exp(2\tau)-1)^{2}-2^{7/2}\gamma_{\tau}\exp(4\tau)}}\leq\frac{\exp(1/2)\cdot\gamma_{\tau}}{\sqrt[4]{1/n^{2}-2^{7/2}\gamma_{\tau}\exp(2/n)}}\leq\frac{\exp(1/2)\gamma_{\tau}}{\sqrt[4]{\frac{1}{n^{2}}-\frac{2^{7/2}\exp(2/5)}{36n^{2}}}}<2\sqrt{n}\gamma_{\tau}.$

Assuming $n\geq 5$ and $d\geq
2+\left(\frac{2}{\pi}\log(8\cdot(6n)^{2})\right)^{2}$, we have