Computer Science Mathematics
Computer Science-2021.gz:143156
Computer Science
65ca88b4d9c3891c96811c13fd7e46cad20ad789
D3p -- A Python Package for Differentially-Private Probabilistic Programming

We present d3p, a software package designed to help fielding runtime efficient widely-applicable Bayesian inference under differential privacy guarantees. d3p achieves general applicability to a wide range of probabilistic modelling problems by implementing the differentially private variational inference algorithm, allowing users to fit any parametric probabilistic model with a differentiable density function. d3p adopts the probabilistic programming paradigm as a powerful way for the user to flexibly define such models. We demonstrate the use of our software on a hierarchical logistic regression example, showing the expressiveness of the modelling approach as well as the ease of running the parameter inference. We also perform an empirical evaluation of the runtime of the private inference on a complex model and find a $\sim$10 fold speed-up compared to an implementation using TensorFlow Privacy.


Introduction
Probabilistic modelling presents a natural way to model data by describing their (assumed) generative process. The model is then fit to observations by probabilistic inference algorithms. Probabilistic programming aims to make the process easy by allowing the user to only specify the model while the system manages the inference process. Probabilistic programming frameworks such as Stan [6], Pyro [3] and PyMC3 [14] have become popular, but they currently offer no support for privacy-preserving algorithms, which are needed for learning from sensitive data.
Differential privacy (DP) [10] provides a rigorous mathematical framework for addressing privacy concerns and has become the de-facto standard notion for privacy in machine learning. It essentially assures that an algorithm's outputs will not differ significantly whether a specific individual's data record is included in the data set or not. Unfortunately, differentially-private algorithms are usually more complex than their non-private counterparts. Software support for easily performing fast differentially-private inference is therefore a crucial tool to achieve privacy-preserving probabilistic programming. This will greatly simplify applications such as differentially-private data anonymisation using a generative probabilistic model to publish a privacy-preserving synthetic twin of a sensitive data set [19].
Using existing probabilistic programming frameworks with privacy-preserving inference is a highly non-trivial task. Practitioners are forced to come up with their own implementation, either from scratch or by adapting existing privacyenabling libraries, which can be an onerous process and leads to many users having to implement the same (or quite similar) wrapper code. There is therefore a clear need for software solutions that enable privacy-preserving probabilistic programming in a convenient and runtime efficient way to allow for fast prototyping and development involving probabilistic programming under privacy constraints.
We address this gap and extend the tool set for growing adoption of DP by introducing an open-source Python software package called d3p. 1 d3p focuses on providing a reliable high-performance implementation of differentially-private doubly stochastic variational inference (DP-VI) [18] for tabular data, where each record corresponds to a single indi-vidual. d3p extends the NumPyro probabilistic programming framework [32,3], allowing modellers to express and fit a large class of probabilistic models under strict privacy guarantees. Alongside the fundamental DP-VI algorithm, d3p uses a state-of-the-art privacy accounting technique [24] to compute tight bounds on the privacy parameters, allowing it to achieve higher levels of utility than with other commonly employed accountants.
Behind the scenes, d3p relies on the JAX framework [4] to perform computation on GPUs and implements a GPUoptimised minibatch sampling algorithm to further optimise performance. Using d3p we achieve a ∼10 fold speedup for fitting a variational auto-encoder model [23] compared to a similar implementation using TensorFlow Privacy [33] on modern GPUs.
d3p addresses a research audience of probabilistic modelling practitioners working with sensitive data. We aim to provide a helpful tool for experimental modelling under privacy constraints. Our main focus in its design therefore is on usability and runtime performance to enable fast modelling iterations. Due to this, d3p currently does not address technical issues arising from implementing idealised differentially private algorithms on machines with imperfect sources of randomness and finite precision, discussed further in Sec. 2.5; these could theoretically be exploited by an adversary if deployed in a production setting.
In summary, we contribute a versatile and performant off-the-shelf implementation of a privacy-preserving probabilistic programming framework as a solid basis for further research. Additionally, we introduce a highly performant subsampling approach based on a slight modification of the CUDA-Shuffle [29], a recently introduced GPU-optimised shuffling approach, and provide a (probabilistic) runtime analysis for it as a minor contribution in methods.
The remainder of the paper is organised as follows: Section 2 reviews probabilistic programming and differentially private variational inference. Based on that discussion we identify software requirements for our software package to clearly outline our design considerations in the same section. We then demonstrate use of our software on a non-trivial hierarchical logistic regression example, illustrating the expressiveness of the probabilistic programming approach and probabilistic modelling (Section 3). Section 4 highlights some implementation details that are orthogonal to the private inference algorithm that forms the core of our framework but that we consider interesting for the user. This includes the discussion of and establishing of (probabilistic) runtime bounds for a special case of the CUDA-Shuffle algorithm which enables GPU-optimised minibatch sampling in our software. Finally, Section 5 presents an evaluation of the d3p framework, including a runtime comparison to a TensorFlow-based implementation, a demonstration of the model introduced in Section 3 and a replication of an experiment for the DP-VI algorithm in [18].

Differentially Private Probabilistic Programming
In this section, we review the background and techniques for differentially private probabilistic programming that inform the implementation choices of our framework. We start with a broad general introduction of the probabilistic programming paradigm and variational inference (Sec. 2.1) and the definition of differential privacy as our main privacy formalism (Sec. 2.2). Following that we give an outline of the powerful DP-VI private inference algorithm (Sec. 2.3) and review privacy accounting tools (Sec. 2.4). Finally we briefly point out technical difficulties resulting from implementing idealised differentially private algorithms on machines with imperfect sources of randomness and finite precision (Sec. 2.5). Each subsection provides an overview of the topic and allows us to identify major requirements for a software implementation, which are summarised in Table 1. The requirements we identify correspond directly to our overarching goals of providing software that is convenient to use and highly performant. Accordingly, we categorise requirements in Table 1 with the labels usability or performance. We hope that explicitly stating our design goals here will allow the reader to evaluate whether our design goals are suitable for their use case and make an informed decision on whether to use d3p. We also do so to emphasise that the implementation of a software package for general use must consider other factors than a (prototypical) implementation of a newly devised method for a research paper.

Probabilistic Programming
Probabilistic programming is a programming paradigm in which a user programmatically defines a statistical model of data which often depends on a set of parameters θ. In mathematical terms, such a model determines a probabilistic density function p(·|θ). A probabilistic inference algorithm is then used to determine the posterior distribution p(θ|X) of the parameter values given a training data set X. The posterior is given by Bayes' formula as p(θ|X) ∝ p(X|θ)p(θ), where, p(X|θ) is the likelihood of the data under the probabilistic model with given parameter values. p(θ) is a prior distribution encapsulating existing knowledge about plausible parameter values.
There exist a number of different probabilistic programming languages that use different ways of specifying the model. To enable easy user adoption, a software package providing differential privacy for probabilistic programming should  Table 1: Requirements for the differentially private probabilistic programming framework with their corresponding category and the subsection they were discussed in (in order of appearance in the text).
not aim to re-define and re-implement yet another solution but rather integrate with existing solutions by extending them with support for differentially private inference (Table 1, Requirement I).
Defining an example model To illustrate the concept of probabilistic programming, we present as an example the implementation of logistic regression for binary classification in NumPyro.
The simple logistic regression model we consider first is for a data set X of records x i ∈ R D with corresponding labels y i ∈ {0, 1}. We assume that each such record and label corresponds to a single individual.
Mathematically the logistic regression model is formulated as ) denotes the sigmoid function and Bernoulli(·; θ) denotes the Bernoulli distribution with success probability θ. The model uses a weight vector w ∈ R D to express the relationship between the records and labels. Bayesian treatment allows us to formulate a prior on w to express any prior knowledge about plausible parameter values. Here we use a weakly informative, zero-centered Gaussian prior with w 0 = 0 and S 0 = 4I. Doubly Stochastic Variational Inference At the heart of probabilistic programming lies the inference algorithm that is used to determine the posterior distribution of parameters p(θ|X). For complex models computing this posterior exactly is typically intractable. Variational inference [20,40], a class of approximate inference algorithms, therefore approximates it with a simpler, tractable distribution q(θ|ψ). The parameters ψ of q are found by solving the optimisation problem min ψ D(q(θ|ψ)||p(θ|X)) where D(·||·) is a divergence measure for probability distributions.
Doubly stochastic variational inference (DSVI) [38] is a gradient ascent algorithm for non-conjugate models with differentiable (joint) probability densities p(X, θ) and a q(θ|ψ) from which values can be easily sampled algorithmically. While these conditions limit applicability somewhat, they still allow for a large class of models to be fitted.
DSVI minimises the KL-divergence by maximising the so-called evidence lower bound (ELBO), defined as The expectation is approximated stochastically by sampling θ from q(θ|ψ) and using a minibatch of the training data for each gradient step. For details we refer to [38].

Differential Privacy
We rely on (approximate) differential privacy [10,9] as the primary privacy notion for a privacy-preserving variant of the DSVI algorithm. Following [8,Def. 2.4], it is defined as: Definition 1 (Approximate Differential Privacy). A randomised algorithm M satisfies (ε, δ)-differential privacy with ε > 0 and 0 ≤ δ ≤ 1 if, for all neighbouring data sets X ∼ X , and for all S ⊂ im(M), we have The two data sets X and X are considered to be neighbours, denoted X ∼ X , when we can obtain one from the other by adding (resp. removing) a single element. ε and δ are privacy bounds (or privacy parameters) restricting the effect that the presence of any particular record in the input data set has on the output of the algorithm M.
Smaller values for these privacy bounds correspond to stricter privacy, however there is a trade-off between privacy and utility of the algorithm's outputs. Larger values for privacy bounds typically result in higher utility of the outputs, as they allow more information to pass through the algorithm. Choosing the privacy bounds therefore requires careful consideration of this trade-off. This is difficult because for many users, especially those inexperienced with DP, it is not clear how to interpret the privacy bounds in a concrete setting. The software should therefore assist the user in choosing appropriate privacy bounds which we reflect in Requirement II.

Differentially Private Doubly Stochastic Variational Inference
Jälkö et al. introduced a (ε, δ)-DP version of the doubly stochastic variational inference algorithm in [18]. This DP-VI algorithm is derived from the influential DP-SGD [34,2] and the relevant steps of a single iteration (out of T many) can be summarised as 1. Sample a random minibatch of size B from the training data set.
3. For each instance in the minibatch: 1. Compute the gradient of the ELBO. 2. Clip the norm of the gradient to a bound C.
6. Update the model parameters ψ with the perturbed gradient.
The main mechanism by which differential privacy is achieved is the perturbation of the minibatch gradient in step 5 via the Gaussian mechanism [8,Thm. 3.22]. The level of noise, characterised by its variance σ 2 , must be carefully calibrated to provide the desired level of privacy. However, the gradient of any data instance could in theory be arbitrarily large, rendering any fixed noise level ineffective. To remedy this, the DP-VI algorithm enforces an upper bound C on the gradient of each data instance (in step 2.2). The important implication of this is that an implementation of DP-VI needs an efficient way of computing and manipulating the per-instance (often also known as per-example) gradients in a minibatch instead of a single gradient over the entire minibatch (Requirement III).
Another important observation is that the DP-VI algorithm has additional hyperparameters C and σ which govern the privacy vs. accuracy trade-off. Especially σ depends non-trivially on the clipping bound C, desired privacy bounds ε and δ, batch size B and the number of iterations T . The next requirement for the software package is therefore the ability to (automatically) derive appropriate values for the DP-VI hyperparameters from these other hyperparameters (Requirement IV).
A final crucial point is that the algorithm is shown to provide differential privacy only under the assumption that minibatches are independently sampled from the training set. As this needs to occur in every iteration of the algorithm, this routine must be especially fast to not slow down the inference as a whole, making another performance requirement for our software (Requirement V).

Privacy Accounting
As we have seen, the DP-VI algorithm consists of an iterative application of the (subsampled) Gaussian mechanism on gradients wrt. random minibatches of the training data. The overall privacy bounds ε and δ of the DP-VI algorithm then result from the DP composition. In order to achieve good utility it is crucial to compute these overall privacy bounds to be as tight as possible: Looser bounds mean that larger perturbations are required for a desired level of privacy, reducing the information extracted from the data and decreasing utility of the inferred model (cf. [11]).
While loose bounds can be computed using general DP composition theorems in a simple way (cf. [8]), obtaining tight bounds typically requires more complex computation using methods called privacy accountants. Abadi et al.'s Moments accountant [2] was the first of these and significantly improved over traditional DP composition theorems. The tightest privacy bounds are currently achieved by the Fourier accountant [24].
Privacy accountants are typically of the form f P A (C, σ, B, T, δ) = ε, i.e., they take in the algorithm's parameters as well as a target value for δ and compute the corresponding upper bound for ε. They are therefore the primary tool to translate between privacy bounds and inference hyperparameters and instrumental for addressing Requirements II and IV.
With these considerations, providing an implementation of a state-of-the-art privacy accountant is an important aspect for a differentially private probabilistic programming framework and becomes Requirement VI.

Remaining Technical Concerns for the Practical Implementation
The definition of approximate differential privacy given in Section 2.2 provides information-theoretic guarantees: It is impossible for the output probabilities of an algorithm to vary too much no matter what the input is. These can typically not be achieved by a computer system which do not have access to perfect sources of randomness for sampling noise and rely on finite-precision approximation of real numbers. Both of these have the potential to completely void the privacy guarantees of DP algorithms in practical implementations: Predictable randomness can allow an attacker to remove the perturbations [16] and finite-precision floating point numbers can leak information due to approximation errors [28].
We recognise these issues as generally important for production systems. However, as we have already pointed out, we undergo no effort to address these for the current version of d3p, which is primarily intended as a research tool. We consider solving these issues for the DP-VI algorithm as important future work.

Summary
We have seen in this section that there is a large number of desiderata for an implementation of a differentially private probabilistic programming framework. These are summarised in Table 1 and fall into different categories that make some mostly technical considerations for the implementation (e.g., Requirements III, V), while others are important aspects in the design of the user interface (e.g., Requirements II, IV).
In the following sections we first explore how our software package addresses these requirements from a user perspective by implementing an example model. However, some of the technical requirements are not experienced by the user in the programming interface directly and cannot be demonstrated in the example: Whether the implementation is performant (Requirements III, V) has almost no effect on the interface seen by the user but is nevertheless a crucial part of their experience. We therefore briefly discuss some implementation details following the examples to illustrate how the identified requirements were addressed and, following that, provide an empirical evaluation of the runtime performance.

d3p Usage Example
We will now demonstrate our d3p software package on a practical example to show how the previously identified requirements are addressed from a user's perspective. d3p centers around an efficient implementation of the DP-VI algorithm and therefore admits non-conjugate models with differentiable probability densities. It is designed to provide differential privacy guarantees for tabular data, where each individual contributed a single sensitive record and records are assumed conditionally independent. d3p uses NumPyro [32,3] as a modelling language and JAX [4,15] as the underlying computation framework, which offers an API similar to NumPy [17]. We start by fully implementing the simple logistic regression model shown in Section 2.1 to demonstrate the basics of probabilistic programming and show how d3p's DP-VI algorithm is invoked to infer the model's parameters. We then highlight the expressiveness of # variational distribution for w sample('w', Normal(w_loc, w_scale)) Listing 2: Implementation of a simple logistic regression model in NumPyro for d3p.
the probabilistic programming approach by adapting the code to a more complex model that achieves a better fit to the data.

Defining a Model
We recall the simple logistic regression model introduced in Sec. 2.1. The model is for a data set X of records x i ∈ R D with corresponding binary labels y i ∈ {0, 1}. We assume that each such record and label corresponds to a single individual.
As before, the model can be formulated mathematically as where σ(x) = 1/(1 + exp(−x)) denotes the sigmoid function and Bernoulli(·; θ) denotes the Bernoulli distribution with success probability θ. The weight vector w ∈ R D links the records and labels. For a Bayesian treatment we place a prior on w to encode existing knowledge. In this first example we assume that we do not have strong prior knowledge but want to enforce some regularisation, and therefore use a weakly informative, zero-centered Gaussian prior with w 0 = 0 and S 0 = 4I. The model is visualised using plate diagram notation in Figure 1.

Implementation of the Model
The implementation of this model using the NumPyro probabilistic programming framework is reproduced in the model function in the top part of Listing 2. The model is defined as a function taking (a minibatch of) the data as an input and specifies a sampling process in an imperative programming style: Values for w are sampled from the specified Gaussian prior with zero mean and a standard deviation of 4. The values for labels are sampled from the Bernoulli distribution (cf. Eq. 3). However, the sample call for y is conditioned to return the values ys passed into the model function using the obs keyword. This is the mechanism by which the labels are passed into the inference algorithm despite the model being specified from a generative perspective.  NumPyro's plate context manager is used to express the independence assumption for the individual data records. Note that this is an important assumption in the d3p package and must be reflected in the model in this way. 2 This requires the additional argument N to the model, which specifies the total amount of data records in the training data.

Implementation of the Variational Posterior
For inference of the model's parameters in our example, we use independent Gaussian distributions for every data dimension j with parameters µ w,j and σ w,j as the variational approximation to the posterior distribution, i.e., The variational parameters µ w and σ w will be optimised according to the discussion in Section 2.1. The corresponding NumPyro implementation is shown in the guide function, following naming conventions of NumPyro, in the lower portion of Listing 1. We register µ w and σ w as parameters w loc and w scale for the inference algorithm and sample w according to Equation 4 in a vectorised fashion. Note that each σ w,j , the standard deviation of the variational Gaussian, must be a positive number, so we actually register a parameter site named w scale log that we pass through the exponential function to obtain w scale. This allows us to perform the optimisation in an unconstrained space but enforces the positivity constraint for w scale.

Running the Inference
In the previous section we programmatically specified the simple logistic regression model using NumPyro. We now turn to the actual private inference of parameter values using d3p. d3p provides an implementation of the DP-VI algorithm via the DPSVI class, which offers the same interface as NumPyro's implementation of the DSVI algorithm in the SVI class. 3 As discussed in Section 2.3, the DP-VI algorithm must be configured for the desired privacy bounds given the batch size and number of training iterations for the inference (Req. IV) and use independent random minibatches (Req. V). The entire code for running the inference is shown in Listing 3.
We first use d3p's subsample batchify data on the data set which returns a function that efficiently samples and returns independent random minibatches. This function is assigned to get batch in our example code. The additional call to batchify init initialises the internal state of the minibatch sampler.
To instantiate a DPSVI object, the main driver of the inference, we need to supply a value for the privacy noise scale σ. We can obtain a σ appropriate for our desired privacy bounds and training hyperparameters using the approximate sigma function. This function returns an approximate value for σ that is guaranteed to achieve the privacy bounds as measured by the Fourier accountant [24], a state-of-the-art privacy accountant method. We store the result in dp scale. Note that the separate computation of σ is a deliberate choice in d3p. While it would be possible to let the instantiation code of the DPSVI class handle this internally, the current approach allows the user easily provide values for σ different from the ones computed by the approximate sigma function, e.g., for research purposes.
After also instantiating implementations of the ELBO (loss) and an optimiser of our choice (optimiser) using classes provided by NumPyro, we are ready to create the DPSVI object (dpsvi). Similar to the minibatch sampler, the DP-VI algorithm provides an initialisation function that is called to produce a state object (svi state). The state object contains randomness state as well as the current values of the parameters and the state of the optimiser. We can now finally run the inference by repeatedly sampling a batch using the get batch function we obtained before and then calling dpsvi.update. The update method completely encapsulates a single iteration of the DP-VI algorithm, including the performant per-instance gradient computation (Req. III), clipping, perturbation and the update of current parameter estimates by the optimiser.
Obtaining parameter estimates from the inference algorithm (and therefore the approximate posterior distribution q(w) in our example model) completes our example at this point. The user can now use standard NumPyro code to interact with the model and the inferred parameters without additional privacy leakage due to DP's invariance to post-processing.

Switching to a More Complex Model
One of the main benefits of the probabilistic programming approach is the ability to easily tailor the model complexity to the information needs and the available prior knowledge and clearly specifying how the respective components of the model relate to each other. So far we have looked at a very simple example where our model makes the implicit assumption that the data are homogeneous and a single parameter vector w describes their relation to the labels equally well for all records. In reality, however, we often face tasks where data comes from different sources that have different local distributions for records. For example, data records containing information about wealth and income of individual persons from different countries are likely to be heavily influenced by the average level of wealth in the respective country. To address this case, we can use a hierarchical logistic regression model like the one considered in [42].
Extending the Model Specification We now extend the previous model notation by vectors g l ∈ R K of group characteristics for L groups and, for each record x i , an indicator l i assigning it to one of the groups. We assume that the records x i and labels y i are sensitive but the group vectors g l are not. Following [42] we use a separate weight vector w l to model the relationship between data record x i and label y i within each group in the same way we did in the simple logistic regression model. However, we now use a hierarchical Gaussian prior centered at M g l for each weight vector w l . The matrix M ∈ R D×K is a new parameter capturing the relation between group characteristics and their corresponding weight vector w l . The full model is visualised in Figure 2 and defined by For simplicity, we consider the covariance matrix Σ l for the distribution of the w l a fixed model parameter of value Σ l = I. We assume an independent weakly informative Gaussian prior for each element of matrix M with µ 0 = 0 and σ 0 = 4 as before. Model Implementation The required changes to adapt our existing model implementation in function guide are straightforward and essentially follow one to one from the textual description above, as shown in Listing 4. First we sample a value for M according to our prior (M), compute the values for η l (etas) and sample the group-specific weight vectors w l (ws) in a vectorised manner. The remaining code implementing the logistic regression for individual records is almost identical to the implementation of the simple model, except that we use the group indicators l i provided in ls to select the entry in ws that corresponds to the w li of the group each record belongs to.
Variational Posterior For the hierarchical logistic regression model, we are interested in a variational approximation to the posterior of the matrix M . Again we can use independent Gaussian distributions for each dimension of M : where µ kd and σ kd for all 1 ≤ k ≤ K and 1 ≤ d ≤ D are the variational parameters.
The implementation of this is shown in the guide function in Listing 4 and is almost identical to the one for the simple model. Invocation of the inference algorithm for the new model also does not change compared to the simple logistic regression model we have looked at before, except for passing the additional data to the model and guide functions, making it especially easy and fast to refine and tweak models.
Note that we do not specify posterior parameters for the group weight vectors w l . This is because in our model these are determined by M , the group characteristics g l and the known covariance matrix Σ l . If we are interested in obtaining values for the w l , we can use NumPyro routines to sample them from the posterior predictive distribution q(w l |g l ) = p(w l |M g l , Σ l )q(M |µ, σ)dM after inference of the variational parameters for M . In the above we slightly abuse notation to denote by q(M |µ, σ) = Listing 4: Implementation of a hierarchical logistic regression model in NumPyro for d3p.

Discussion
As we have seen we were able to expand our model to the additional structure within our data with a few simple changes in the implementation. The changes correspond directly to the textual description of the new hierarchical model, making it easy for the user to translate theory into implementation. While the initial simple logistic regression model is quite a common choice and specialised implementations for this exist in many software packages, extending the model in a similar fashion to what we did in this section is often more complicated, if at all possible, in those.
Note that we still have made a couple of simplifying assumptions here, one of which is the assumption that the covariance Σ l of the distribution of group weight vectors is a known constant. This was merely for convenience and not because of restrictions of model expressiveness. We can easily formulate a prior and a variational posterior for Σ l to learn it from the data, if this assumption does not hold. Another simplifying modelling assumption is the choice of the prior for M , however, we consider selecting more informative priors outside the scope of this example. NumPyro offers a wide range of distributions from which the user can choose adequate priors according to their modelling needs. These two brief notes serve to show that the probabilistic programming paradigm allows the user to make finegrained decisions on the model expressiveness they need by flexibly either excluding details they are not interested in or incorporating detailed prior knowledge in the model.

d3p Implementation Outline
We now discuss some technical details of our d3p package and how it meets the requirements identified in Section 2.
In line with our discussion of these requirements, we hope that this helps the reader decide whether the choices made in implementing d3p make it suitable for their use. We also focus on the GPU-optimised minibatch sampling algorithm in particular, which we base on a recent GPU-optimised shuffling method but slightly modify to better fit our use.
We organise this section by first briefly motivating the use of NumPyro as the basis for d3p (Section 4.1). We then turn to the implementation of the DP-VI algorithm for NumPyro (Section 4.2) with a focus on how it realises its high performance for per-instance gradient computation (Req. III).
As the second major implementation detail we then discuss our GPU-optimised i.i.d. minibatch sampling routine (Req. V) which has a major impact on the overall performance of the inference and for which we also contribute a runtime analysis (Section 4.3). Finally we finish the discussion on implementation by a brief overview of the remaining usability goals in the implementation (Section 4.4).

Underlying Framework
To meet Requirement I (Integrate with an existing popular probabilistic programming framework), we chose NumPyro as the underlying probabilistic programming framework for d3p. It is a spin-off of the popular Pyro [3] framework, providing a very similar API, but relies on Google's JAX [4,15] for the underlying computational optimisation functionality. The JAX framework uses tracing mechanisms to compile pure (side-effect free) functions directly from Python code to efficient XLA kernels that can be run on either CPU or GPU. Implementations of the DP-SGD algorithm in JAX where found to be consistently faster than competing implementations [36], making it a promising backend for for d3p.
The choice of NumPyro as basis for d3p is therefore a compromise between fulfilling Requirement I and being able to provide a high-performance implementation (Req. III, V). While targeting Pyro instead of NumPyro would arguably have enabled access to a larger established user base, implementing per-instance gradients in the underlying PyTorch would have decreased performance of d3p.

Implementation of DP-VI
As we have seen in the examples, d3p centers around the DPSVI class implementing the DP-VI algorithm. DPSVI offers the same interface as NumPyro's non-private implementation in the SVI class and therefore works as a drop-in replacement for it. This allows for especially easy adoption of privacy-preserving methods with minimal required changes to the code base.
We have identified in Section 2.3 that the DP-VI algorithm requires computation of per-instance gradients. Unfortunately, these are typically not readily available in established machine learning frameworks. A naive but inefficient solution is to set the batch size B to one, which gives per-instance gradients at the cost of losing the performance gains resulting from parallel computation on minibatches that are crucial for efficient machine learning applications.
As the JAX framework is based on the manipulation of side-effect free functions, it also provides a range of composable higher order transformations for these. Crucially, the computation of gradients as well as vectorising a function for parallel execution are examples of these higher order transformations. This enables d3p's implementation in the DPSVI's class to efficiently parallelise computation of gradients and the subsequent clipping over a minibatch using vectorisation, i.e., single-instruction-multiple-data (SIMD) style computing. Leveraging the massively parallel processing capabilities of modern GPUs in this way allows us to negate most of the additional overhead introduced by the per-instance gradient computation.
This vectorisation approach effectively turns computation on a batch of size B into B parallel computations on batches of size 1. In order to prevent this from affecting the relative contributions of prior and variational posterior of global parameters, models must use NumPyro's plate environment to scale likelihood contributions from the batch appropriately to the perceived batch size. This typically means that each data record's contribution to the loss is scaled up. The DPSVI class therefore performs some crucial bookkeeping to ensure that the privacy perturbation in each iteration is also scaled appropriately.
As a convenience feature, the DPSVI class offers methods to obtain tight privacy bounds for its current hyperparameter values which are computed using the Fourier accountant [24]. 4

Performant GPU Batch Subsampling
We have seen in Section 2.3 that the DP-VI algorithm's privacy guarantees rely crucially on minibatches sampled from the data set in a truly i.i.d. fashion (Requirement V). This is an important difference to non-private algorithms that can usually get away with a permute-and-iterate approach to sample minibatches. Often, the data set is permuted once in its entirety and minibatches are then consumed by iterating over the permuted set. Due to being invoked comparatively rarely, the performance of the permutation algorithm does not make a noticeable difference on the overall runtime. For the same reason, it can also run on a different device than the learning algorithm as slow bus transfers are infrequent.
The i.i.d. requirement for DP-VI requires that the minibatch sampling routine is invoked once for each iteration, making the cost of a slow sampling routine prohibitive (or at least, much more noticeable). d3p therefore ships with a parallel and GPU-optimised minibatch sampling routine based on a novel shuffling methodology (CUDA-Shuffle) proposed in [29,35].
Conventional shuffling algorithms, such as the Fisher-Yates shuffle, are ill-suited for GPU-acceleration as they are sequential. In contrast, the main idea in the parallel shuffle algorithm is to use a bijective function f k that for a given key k defines a unique pseudo-random mapping on sets of indices I B → I X . Assuming I B = {0, . . . B − 1} to describe indices in a minibatch and I X = {0, . . . , n − 1} indices in the data set, f k allows sampling a minibatch of elements from the input set in parallel without collisions. Previously, it has been proven that a Feistel network [13] with more than two rounds is a pseudo-random bijective function, provided that it uses a round function that is pseudo-random and the set size is a power of two, i.e., n = 2 b [27].
In [29] this is generalised for arbitrary n by taking the smallest bit-length b such that 2 b > n, applying the bijection on the index set of length 2 b and then removing all values larger than n by an efficient GPU compaction algorithm. The overall runtime of this is in Θ(n).
However, for our application of sampling a minibatch, shuffling the entire index set is inefficient. Instead, we apply the Feistel network repeatedly on values from I B until all outputs are a value in I X . Our generalised Feistel permutation generator can be given as where l is the iteration count and i ∈ I B . Sampling a minibatch then requires lB evaluations of the Feistel network for some factor l. These evaluations consist only of independent, parallel bitwise operations, making this approach very well-suited for GPU acceleration.
The factor l is the number of iterations that are required to converge Equation 8, for which we will now establish a probabilistic upper bound. The probability with which the Feistel network f k returns an output f k (x) < n, given a bit-length b such that where we let r = n − 2 b−1 , i.e., the non-power-of-two residue of n.
Now, let L i be the random variable of the iteration count required to converge x i < n, ∀i ∈ I B and F i = L i − 1 be the penultimate iteration count where some x i ≥ n. F i follows a negative binomial distribution, P r[ Its expected value is and the cumulative distribution is In the worst case with r = 1, resulting in the largest gap to the next power of two, the expected values are Hence with non-trivial data set sizes with b 1, the expected number of iterations required to converge the permutation is E[L i ] ≈ 2, in the worst case. To estimate the maximum number of iterations for a given percentage θ of cases, we can calculate Setting θ = 0.99, b 1 and the worst case r = 1, we obtain Therefore, in 99% of cases we will see no more than six failures, i.e., seven total iterations (and no more than six in 95% of cases). In practice, as long as the batch size is sufficiently small, precisely q = B/n < 1/7, our approach will require less evaluations of the bijection than the one of [29].

Privacy Bound and Hyperparameter Selection
d3p offers an API to compute the perturbation hyperparameter σ for the DP-VI algorithm via the approximate sigma function to satisfy Requirement IV (Determine DP inference algorithm parameters automatically.) This is realised by employing standard black-box optimisation techniques to find a suitable input such that the Fourier accountant arrives at the desired value for ε, given all other hyperparameters.
Finally, Requirement II (Provide assistance for the user in finding adequate privacy bounds) is an important piece of guidance for users inexperienced with differential privacy. Optimal choice of ε and δ depends on a delicate balance between the desired utility and the level of risk of privacy violation that is considered adequate by the user -a choice, therefore, that can only be made by the user but requires knowledge of how ε, δ relate to concrete privacy risks. Unfortunately, this relationship is still an open research question, which makes it difficult to give detailed guidance to the user. Common practice is to require that δ < 1 N , where N is the number of individual records in the data, and ε ≤ 1, however d3p does not enforce this currently to allow for free experimentation.

Evaluation
To demonstrate the performance and flexibility of our framework, we explore a few examples in this section. We first compare d3p's runtime performance to TensorFlow on the implementation of a variational auto-encoder (Sec. 5.1). Afterwards, we show some results for the hierarchical logistic regression model discussed in detail in Sec. 3.3 and use that opportunity to highlight some privacy trade-offs in the regime of small data sets (Sec. 5.3). Finally, we briefly compare results of d3p to the original implementation of the DP-VI algorithm in [18] (Sec. 5.4) on a Gaussian mixture model.

Comparison with TensorFlow Privacy
We first compare the performance of our d3p package to an implementation of an identical model using the TensorFlow Probability framework [7] with a manual implementation of the variational inference algorithm and privacy enabled by the TensorFlow Privacy package [33].
We choose a variational auto-encoder (VAE) model [23] for this purpose. VAEs are generative models that consist of an encoder function, mapping data to parameters of a distribution on latent representations, and a decoder function, converting samples from the latent space to data samples. Generating data consists of drawing a sample from the distribution in the latent space and passing it through the decoder function. These mapping functions are represented by neural networks and therefore typically have a large number of parameters.
VAEs therefore provide an excellent test case for performance and have been previously used for the same purpose [3]. We use slight variations of the same model on different image classification data sets, namely MNIST [26], Fashion-MNIST [43] and CIFAR-10 [25]. For MNIST and Fashion-MNIST datasets, we use feed-forward networks with a single hidden layer encoder and decoder, consisting of 688 884 trainable parameters in total. For CIFAR-10 dataset, we employ networks of 3 convolutional layers followed by a single dense layer, consisting of 640 423 parameters in total.
Of particular interest are the respective runtime of the inference as well as a comparison of the inferred model to verify that d3p is fast and accurate. Table 2 shows the runtime per epoch as well as the loss on the held-out test set after 20 epochs, i.e., passes over the training data set of 60 000 images (50 000 for CIFAR-10). We use a minibatch size of 128. The reported numbers are averages over 20 training processes that were run using a single NVIDIA Tesla V100 32G GPU. All runs used the Adam optimiser [22] for parameter updates after computing gradients. For the DP variants we used σ = 1.5, resulting in ε ≈ 0.5 for δ = 1/60000 (δ = 1/50000 for CIFAR).
From the comparison, we can see that the end losses with MNIST and Fashion-MNIST data sets are comparable across frameworks, and all cases were observed to converge well. However, we found CIFAR-10 data set to be more challenging to train due to large variation in samples, relatively modest number of samples per class and three colour channels. Neither framework was able to learn good representations of CIFAR-10 with DP for our choice of hyperparameters. The aforementioned challenges may also explain the higher discrepancy between the non-DP losses. In terms of performance, d3p consistently outperforms the TensorFlow (TF) implementation by a factor of ∼10 for all data sets. Additionally, the relative performance loss of DP-VI compared to non-private inference in the same framework is lower in d3p (up to ∼ 2.5-fold for CIFAR-10) compared to the TF implementation (up to ∼ 22-fold).
For d3p, we also compare using the Feistel-based GPU-optimised minibatch sampler against using JAX's built-in jax.random.choice method and summarise the results in Table 3. We observe that our optimised sampler consis-   tently yields a ∼ 100 ms speed-up (15% on MNIST) per epoch on all data sets. This similarity is due to the similar sizes of the data sets resulting in similar amounts of total iterations and thus invocations of the subsampling.

Effect of Batch Size
While the vectorised implementation of the per-instance gradient computation and manipulation eliminates much of the time overhead required in the DP-VI algorithm, it does not remove it completely. This is due to the additional steps in DP-VI, such as gradient clipping, summing and perturbing, but may also be an effect of a larger memory footprint due to holding gradient values for each data instance intermittently. This necessitates a larger amount of memory accesses, which can slow down computation. We therefore investigate the effect of minibatch size on the runtime of the DP-VI implementation in d3p for MNIST. In Figure 3 we plot the average runtime per iteration (over 100 iterations) over the size of minibatches for DP-VI in d3p and non-private inference in NumPyro. We see that, as expected, runtime per iteration increases more steeply for DP-VI than the non-private case. However, we also observe that DP-VI runtime increases only linearly with minibatch size, which is in line with our expectations.

Hierarchical Logistic Regression and Small Data
We next evaluate the implementation of the hierarchical logistic regression model presented in Section 3.3 with small data sets. These are generally problematic in privacy-preserving inference because every single data point has a comparatively larger effect on the outcome. Achieving good model performance therefore requires looser privacy bounds than in the case with larger data sets. d3p makes it easy for the user explore the trade-off between privacy and utility via its fast inference and the ease of adjusting the algorithm by simply specifying the privacy parameters from which the perturbation noise is automatically derived.
In the following experiments we use synthetic training data sets of varying size N that follow the hierarchical structure described in Section 3.3. Each data point is five-dimensional and points are split into L = 3 groups that are in turn described by K = 3 variables each. To evaluate the goodness-of-fit of the trained model we evaluate the area under the ROC curve (AUC) on a held-out test set of the same size N as the training set. Figure 4 shows the AUC after training on a data set with N = 500 data points for 100 000 iterations. Results are shown for different levels of privacy as well as non-private variational inference for ten runs with different random seeds. Privacy bound δ was kept fixed to 1 N . Each run took less than 15 seconds on a commodity laptop without GPU acceleration. The red line is the AUC for non-privately fitting a simple (non-hierarchical) logistic regression model   using scikit-learn [31]. The trade-off between privacy and utility is clearly visible: Smaller values of ε corresponding to stricter privacy constraints lead to lower AUC on average and a larger spread of results over different runs. Runs for ε = 4 or larger are close to the AUC of the non-private model on average but exhibit larger spread. Runs for ε = 2 fall short of this but still outperform the simpler baseline. This highlights a strength of the probabilistic modelling approach particularly relevant for privacy-preserving machine learning: By encoding prior knowledge of the generative process underlying the data in a principled way, privacy budget does not need to be spent to learn this a-priori known structure, allowing more capacity to learn the remaining parameters of the model. Figure 5 shows the effect of data set size for the same privacy levels after training for 100 000 iterations. The graphs show the spread of one standard deviation above and below the mean over ten runs for each data set size and privacy level. Smaller data sets result in lower average AUC and larger spread. For data sets of less than 500 records utility deteriorates rapidly.
We finally explore the effects of the number of training iterations. Since the amount of noise added to perturb gradients in the DP-VI algorithm increases with the number of iterations, one could expect that choosing too large a number of iterations will negatively effect the learning. Figure 6 shows the evolution of AUC for N = 500 for different numbers of total training iterations. Note that these are results of separate runs each with a different number of iterations and thus different amounts of perturbation per iteration, not the evolution of results over a single long training run.  Figure 5. The additional horizontal line shows the AUC of the simple non-private non-hierarchical logistic regression baseline.
We observe clearly that due to the privacy perturbations of gradients during training, the DP-VI algorithm takes longer to converge than non-private variational inference. Stricter privacy bounds move convergence to higher iteration counts. Contrary to expectation, despite the larger perturbations required for larger iteration counts, we see a general trend in improved utility for longer training (for N = 100, this trend continues up to 500 000 iterations). In our experiments we observe no negative impact of increasing the iteration count on the final AUC even if the training converges earlier, indicating that the DP-VI algorithm is very robust to the privacy perturbations.

Gaussian Mixture Model
We further demonstrate the ease of specifying expressive models in d3p by replicating an experiment on a Gaussian mixture model from the original DP-VI paper [18]. They used two-dimensional data generated from 5 clusters of spherical Gaussians and trained the model for 1 000 iterations for different levels of privacy. The evaluation is in terms of log-likelihood of a held-out test set on the learned predictive model.  The original code 5 required model specific implementation of the DP-VI algorithm due to the absence of a generic framework for privacy-preserving probabilistic modelling. Using d3p it suffices to simply write out the model and an implementation of log-probability calculations and sampling routines for a Gaussian mixture distribution as an implementation of NumPyro's Distribution class. These are shown in Appendix A. Figure 7 shows the resulting log-likelihood (higher is better) of the test data after training the model using the d3p implementation and the original code of [18]. Shown are the average with standard error over five runs in both cases.
To keep results comparable we use the levels of privacy perturbation (parameter σ) from the original code for d3p.
The results from d3p are very consistent with to slightly better than those from the original implementation but exhibit much less variability due to randomness in the inference algorithm, which is particularly pronounced for ε = 0.1 in the original paper code.

Related Work
Major popular probabilistic programming frameworks for the Python programming language are Edward2 [39], Ten-sorFlow Probability [7], (Num)Pyro [3,32] , PyMC [14], and Stan [6]. Edward2, TensorFlow Probability, Pyro and PyMC allow the user to declare models and run the inference from the Python programming language and differ mainly in the computation framework they rely on to run the inference (Edward2 and TensorFlow Probability use TensorFlow [1], Pyro uses PyTorch [30] and PyMC uses Theano [37]; NumPyro is a direct port of Pyro to JAX). The Stan framework follows a different approach and requires models to be specified in a dedicated domain-specific language, which is then evaluated using the Stan runtime, which can be invoked from Python or other major programming languages. None of these frameworks currently offers support for privacy-preserving inference.
A number of general implementations for differentially private machine learning exist for the popular frameworks. Notable are TensorFlow Privacy [33] for TensorFlow, Opacus [12] and PyVacy [41] for PyTorch. These generally provide implementations of the DP-SGD algorithm as alternative optimisers for the computational framework. In principle they could be combined with the dominant probabilistic programming framework for the respective backend, however this integration is usually not as seamless as one would desire. These implementations of DP also often suffer from poor performance in the implementation that can usually be traced back to inefficient computation of per-instance gradients [36]. With d3p we aim to provide better integration with high performance by directly targeting the NumPyro probabilistic programming framework. We are not aware of any general library of DP-SGD for the JAX framework that could be used with NumPyro to achieve the same goal. As far as we are aware, none of these libraries takes active measures to address the previously discussed technical issues of implementing differential privacy.

Conclusion
We have presented our d3p package which extends the NumPyro probabilistic programming framework with runtime efficient differentially private inference. We demonstrated the use of our framework and the expressiveness of the probabilistic programming approach on an extensive example and highlighted the requirements and corresponding implementation choices for our software. Our goal is to provide a helpful tool that encourages use of probabilistic programming as a viable approach to modelling data for privacy practitioners, as well as lowers the threshold for adoption of privacy-preserving methods for probabilistic modelling experts.