1 INTRODUCTION↩︎

When applying Bayesian inference on sensitive data, one needs to consider the privacy risks the released results might pose. Multiple methods combining the state-of-the-art privacy paradigm differential privacy (DP) [1] with Bayesian inference have been proposed in the past. These include methods that are based on, e.g., Markov chain Monte Carlo (MCMC) [2], [3] and variational inference (VI) [4], [5].

Many of these methods do not consider how the additional noise due to DP affects the learned posteriors, which might lead to poor uncertainty quantification. As a solution, multiple works trying to explicitly model the DP noise have since been proposed [6]–[8], but these are limited to simple or small problems.

In this work, we propose a highly scalable noise-aware DP Bayesian inference approach that is applicable to a broad range of models. Our method uses DP variational inference (DPVI) [4], [5] to obtain noisy gradients from an approximate Bayesian inference problem. Modelling the DP noise in the gradients, we form a probabilistic model connecting the noisy gradients to an optimal solution of the VI problem. Finally, learning the noise-aware posterior for the optimum, we can form a noise-aware posterior for the original Bayesian inference problem.

Related work: Several methods have been proposed to take the DP induced noise into account in private Bayesian inference. Earlier works include sufficient statistic based models such as noise-aware exponential family models [6] and Bayesian linear regression [7]. In these works, the data is released through noisy sufficient statistics, and the perturbation noise is included as a part of the data generation process. Using approximate sufficient statistics, [9] extended the noise-awareness for Bayesian generalized linear models [10]. [11] proposed a general summary statistics–based approach using approximate Bayesian computation (ABC). [8] proposed a Gibbs sampling based approach, where the latent confidential data is augmented into the inference model. Their approach is applicable for a wide range of models; however, the cost of the inference process scales by the number of samples in the data. [12] proposed a method which fits a normalizing flow as a surrogate for the true posterior, by iteratively drawing the model parameters from the surrogate, drawing the latent confidential data conditioned on the proposed parameter, and finally updating the surrogate model based on the observed noisy summary statistics. As an important application of DP probabilistic modelling, [13] proposed a synthetic data generation method, which is based on learning a noise-aware posterior for a discrete marginal query based summary statistic.

Contributions:

2 BACKGROUND↩︎

2.1 Bayesian Inference↩︎

Assume we have a model \(p\left(\mathbf{D}\mid \pmb{\theta}\right)\), where \(\pmb{\theta}\in \Theta \subseteq \mathbb{R}^{n}\) denotes the unobserved model parameters; and \(\mathbf{D}\in \mathcal{D}\) denotes the data. Given a prior \(p(\pmb{\theta})\) for \(\pmb{\theta}\), Bayes’ Theorem states that the posterior can be written as: \[\label{bayes-Theorem} p(\pmb{\theta} \mid \mathbf{D}) = \frac{p\left(\mathbf{D}\mid \pmb{\theta}\right)p(\pmb{\theta})}{\int_{\Theta} p\left(\mathbf{D}\mid \pmb{\theta}\right)p(\pmb{\theta}) \mathrm{d}\pmb{\theta}}.\tag{1}\] In many cases, the denominator in 1 is intractable. For such cases, the true posterior can only be approximated with some other distribution \(\widetilde{p}(\pmb{\theta}\mid \mathbf{D})\) using methods such as Markov chain Monte Carlo (MCMC) [15] or Variational inference (VI) [16].

2.2 Validating Approximate Bayesian Inference↩︎

After obtaining \(\widetilde{p}\), we want to test how well \(\widetilde{p}\) approximates \(p\). In order to devise such a test, let us first assume that \(\widetilde{p}\) has support anywhere \(p\) has support. Next, we define a \((1-\alpha)\) credible region for \(\tilde{p}\) as a mapping \(\widetilde{\mathcal{R}}_{\alpha}: \mathcal{D} \rightarrow \mathcal{P}(\Theta)\), where \(\mathcal{P}\) denotes a power set, such that for any \(\mathbf{D}\in \mathcal{D}\) \[\label{approximate-posterior-coverage} \int_{\Theta} \text{\usefont{U}{bbold}{m}{n}1}_{\widetilde{\mathcal{R}}_{\alpha}\left(\mathbf{D}\right)}\left( \pmb{\theta}\right) \widetilde{p}\left(\pmb{\theta}\mid \mathbf{D}\right) \mathrm{d}\pmb{\theta} = 1 - \alpha.\tag{2}\] Following [14], we define the Expected Coverage Probability (ECP) of \(\widetilde{\mathcal{R}}_{\alpha}\left(\mathbf{D}\right)\) as \[\begin{align} \label{ecp-definition-original} \text{ECP} \left[\widetilde{\mathcal{R}}_{\alpha}\left(\mathbf{D}\right)\right] &= \underset{\pmb{\theta}, \mathbf{D}\sim p\left(\pmb{\theta}, \mathbf{D}\right)}{\mathbb{E}}\left[\text{\usefont{U}{bbold}{m}{n}1}_{\widetilde{\mathcal{R}}_{\alpha}\left(\mathbf{D}\right)}\left( \pmb{\theta}\right) \right]. \end{align}\tag{3}\] [14] showed, that if \(\text{ECP} \left[\widetilde{\mathcal{R}}_{\alpha}\left(\mathbf{D}\right)\right] = 1-\alpha\) for every mapping \(\widetilde{\mathcal{R}}_{\alpha}\) that satisfied 2 and every \(\alpha \in (0, 1)\), we have \(\tilde{p}(\pmb{\theta}\mid \mathbf{D}) = p(\pmb{\theta}\mid \mathbf{D})\). Unfortunately, testing this is not practically feasible, as we would need to enumerate over every possible credible region. In practice, we need to perform the test over a finite selection of credible regions, and see if the expected coverage over this set is close to \(1-\alpha\). However, the way the set of credible regions is constructed is critical as shown by [14]. For example, we can obtain perfect \(1-\alpha\) coverage by choosing \(\widetilde{\mathcal{R}}_{\alpha}\left(\mathbf{D}\right)\) to be the (\(1 - \alpha\)) Highest Posterior Density (HPD) region and then setting \(\widetilde{p}\left(\pmb{\theta}\mid \mathbf{D}\right) = p(\pmb{\theta})\), i.e.using the prior as the posterior approximation. To solve this problem, [14] introduced the concept of a positionable credible region.

Intuitively, this means that we can position credible regions around any point in \(\Theta\). In addition, we want to choose \(\pmb{\theta}_{\text{ref}}\) as a function of \(\mathbf{D}\); that is, \(\pmb{\theta}_{\text{ref}}: \mathcal{D} \rightarrow \Theta\). Based on Theorem 3 in [14], if \[\text{ECP}\left[\widetilde{\mathcal{C}}_{\alpha}\left(\mathbf{D}, \pmb{\theta}_{\text{ref}}\left(\mathbf{D}\right)\right)\right] = 1 - \alpha,\] for all \(\alpha \in (0, 1)\) and any function \(\pmb{\theta}_{\text{ref}}: \mathcal{D} \rightarrow \Theta\), then \(\widetilde{p}\left(\pmb{\theta}\mid \mathbf{D}\right) = p\left(\pmb{\theta}\mid \mathbf{D}\right)\) for all \(\pmb{\theta}\in \Theta\) and \(\mathbf{D}\in \mathcal{D}\). This is an important property of the ECP with the positionable-credible regions, as it establishes the equality of \(\widetilde{p}\) and \(p\) if and only if \(\widetilde{p}\) has perfect coverages. For the rest of the paper we will abbreviate \(\widetilde{\mathcal{C}}_{\alpha}\left(\mathbf{D}, \pmb{\theta}_{\text{ref}}\left(\mathbf{D}\right)\right)\) as \(\widetilde{\mathcal{C}}_{\alpha}\left(\mathbf{D}\right)\).

The resulting coverage test by [14] called TARP can be implemented with the following steps.

2.3 Differential Privacy (DP) and DP Stochastic Gradient Descent↩︎

DP [1] is the standard definition for privacy in modern computer science. DP is defined over so called neighbouring data sets, i.e. data sets that differ in only single element. The formal definition is given as:

DP has a number of appealing theoretical properties, including compositionality, whereby repeated accesses to the data weaken the guarantees in a predictable manner, and post-processing immunity, whereby the guarantees can never be weakened by post-processing.

DP Stochastic Gradient Descent (DP-SGD) [17]–[19] modifies the traditional SGD algorithm to satisfy DP. In order to apply DP-SGD, we need to assume that the loss function, parametrized by \(\pmb{\phi}\), can be represented as a sum over individuals \(\pmb{x}_i \in \mathbf{D}\): \[\label{dp-sgd-loss-definition} \mathcal{L}\left(\pmb{\phi}; \mathbf{D}\right) = \sum_{i= 1}^{N} \ell(\pmb{\phi}; \pmb{x}_i).\tag{4}\] DP-SGD makes three modifications to the SGD algorithm:

The DP-SGD update rule can be written as \[\begin{align} \label{dp-sgd-definition} \begin{aligned} \pmb{g}_{t + 1} &= \sum_{i \in \mathcal{B}_{t + 1}} \mathrm{clip}\left(\nabla_{\pmb{\phi}}\ell(\pmb{\phi}_t; \pmb{x}_i), C\right), \\ \pmb{\phi}_{t + 1} &= \pmb{\phi}_{t} - \lambda_t \left[ \pmb{g}_{t + 1} + \sigma_{\text{DP}} C \pmb{\eta}_{t + 1} \right], \end{aligned} \end{align}\tag{5}\] where \(\lambda_t\) is the learning rate, \(\mathcal{B}_t\) is a Poisson subsampled minibatch and \(\pmb{\eta}_{t+1} \sim \mathcal{N}(\mathbf{0}, \mathbf{I}_d)\). The Gaussian noise-level \(\sigma_{\text{DP}}\) is chosen so that the \((\epsilon, \delta)\)-DP guarantees hold for \(T\) iterations of the DP-SGD algorithm. The general privacy proofs for DP-SGD allow releasing all of the intermediate steps and the noisy gradients of the algorithm under the same privacy guarantee.

2.4 DP Variational Inference↩︎

Variational inference [16] is a method used to approximate a posterior \(p\left(\pmb{\theta}\mid \mathbf{D}\right)\) using another distribution called the variational distribution \(q_{\text{VI}}\left(\pmb{\theta}; \pmb{\phi}\right)\) parameterized by a vector \(\pmb{\phi}\in \Phi \subseteq \mathbb{R}^{d}\). The idea is to find \(\pmb{\phi}\) that maximizes the Evidence Lower Bound (ELBO): \[\begin{align} \label{elbo-definition} \mathcal{L}\left(\pmb{\phi}; \mathbf{D}\right) &= \mathbb{E}_{q_{\text{VI}}(\pmb{\theta}; \pmb{\phi})} \left[ \log \frac{p(\pmb{\theta}, \mathbf{D})}{q_{\text{VI}}(\pmb{\theta}; \pmb{\phi})} \right]. \end{align}\tag{6}\] Furthermore, when the observations are i.i.d. we can decompose the ELBO as \[\begin{align} \begin{aligned} \mathcal{L}\left(\pmb{\phi}; \mathbf{D}\right) &= \mathbb{E}_{q_{\text{VI}}(\pmb{\theta}; \pmb{\phi})} \left[ \log \frac{p(\pmb{\theta})}{q_{\text{VI}}(\pmb{\theta}; \pmb{\phi})} \right] \\ &\phantom{aaa} + \sum_{i=1}^N \mathbb{E}_{q_{\text{VI}}(\pmb{\theta}; \pmb{\phi})} \left[ \log p(\pmb{x}_i \mid \pmb{\theta}) \right]. \end{aligned} \end{align}\] The expectations in ELBO are typically intractable. To solve this, [20] proposed a solution called Stochastic Gradient Variational Bayes (SGVB), which approximates the expectations with Monte-Carlo estimates and parametrizes \(q_{\text{VI}}\) s.t. the Monte-Carlo approximator of the ELBO becomes a differentiable function that can be optimized with gradient based methods. For more details on the Monte-Carlo estimation and the differentiable parametrization of \(q_{\text{VI}}\) see 13.

[4] proposed an DP variant of the SGVB method called DPVI, by replacing the gradient optimizer with DP-SGD. Later [5] improved the DPVI by using certain structural knowledge of the gradients in order to improve the convergence.

3 NOISE-AWARE INFERENCE↩︎

3.1 Formalism↩︎

The distributions in this paper are assumed to be continuous unless explicitly stated otherwise. Optimizing the ELBO using DP-SGD introduces uncertainty into the final VI approximation. Our goal is to integrate the noise from the DP mechanism into the posterior, this motivates the following definition.

Similar to \(\mathbf{D}\), the vector \(\pmb{\xi}\) can be thought of as either one point or a list of points concatenated into one large vector. For instance, DP-SGD provides the iterates \((\pmb{\phi}_0, \ldots, \pmb{\phi}_T)\) as a list of vectors. ?? is obtained by applying Bayes’ theorem using the joint distribution in ?? with additionally marginalizing out the data \(\mathbf{D}\), as it cannot be observed in the DP setting.

Computing ?? analytically is often intractable, and instead it is approximated by another distribution. Given any approximation of ?? , say \(\widetilde{p}\left(\pmb{\theta}\mid \pmb{\xi}\right)\), we want to evaluate this approximation using the ECP. However, compared to the non-private Bayesian inference, we now have three random vectors \(\pmb{\theta}\), \(\mathbf{D}\), and \(\pmb{\xi}\) in the data generating process. Thus, we propose a slight modification to the \(\text{ECP}\) as follows. Note that \(p\left(\pmb{\theta}, \pmb{\xi}\right) = p\left(\pmb{\xi}\mid \pmb{\theta}\right) p\left(\pmb{\theta}\right)\), so we can define the ECP equivalently using the joint distribution in ?? as, \[\begin{align} \label{ecp-definition-alternate} \begin{aligned} \text{ECP} \left[\widetilde{\mathcal{C}}_{\alpha}\left(\pmb{\xi}\right)\right] &= \underset{\pmb{\theta}, \mathbf{D}, \pmb{\xi}\sim p\left(\pmb{\theta}, \mathbf{D}, \pmb{\xi}\right)}{\mathbb{E}}\left[\text{\usefont{U}{bbold}{m}{n}1}_{\widetilde{\mathcal{C}}_{\alpha}\left(\pmb{\xi}\right)}\left( \pmb{\theta}\right) \right]. \end{aligned} \end{align}\tag{7}\] The reason why 7 and 3 are equivalent is that the indicator function \(\text{\usefont{U}{bbold}{m}{n}1}_{\widetilde{\mathcal{C}}_{\alpha}\left(\pmb{\xi}\right)}\left( \pmb{\theta}\right)\) does not depend on \(\mathbf{D}\), so by applying Fubini’s Theorem [21], the expectation in 7 becomes \[\begin{align} \underset{\pmb{\theta}, \mathbf{D}, \pmb{\xi}\sim p\left(\pmb{\theta}, \mathbf{D}, \pmb{\xi}\right)}{\mathbb{E}}\left[\text{\usefont{U}{bbold}{m}{n}1}_{\widetilde{\mathcal{C}}_{\alpha}\left(\pmb{\xi}\right)}\left( \pmb{\theta}\right) \right] = \underset{\pmb{\theta}, \pmb{\xi}\sim p\left(\pmb{\theta}, \pmb{\xi}\right)}{\mathbb{E}} \left[\text{\usefont{U}{bbold}{m}{n}1}_{\widetilde{\mathcal{C}}_{\alpha}\left(\pmb{\xi}\right)}(\pmb{\theta}) \right]. \end{align}\] Similarly, the ECP in 7 with the positionable-credible regions establishes the equality of \(\widetilde{p}\) and \(p\) if and only if \(\widetilde{p}\) has perfect coverages.

To evaluate the approximate noise-aware posterior based on 7 , the coverage test in [alg:non-private-ecp] needs to be modified to include the data generation process for \(\pmb{\xi}\). We first draw samples from the joint distribution ?? , where the samples are drawn in the following order, \[\pmb{\theta}_i \sim p\left(\pmb{\theta}\right) \boldsymbol{\rightarrow} \mathbf{D}_i \sim p\left(\mathbf{D}\mid \pmb{\theta}_{i} \right) \boldsymbol{\rightarrow} \pmb{\xi}_i \sim p\left(\pmb{\xi}_{i} \mid \mathbf{D}_{i}, \pmb{\theta}_{i} \right).\] Second, we replace \(I_k\) in [alg:non-private-ecp] with \(I_k = \text{\usefont{U}{bbold}{m}{n}1}_{\widetilde{\mathcal{C}}_{\alpha}\left(\pmb{\xi}_i\right)}\left(\pmb{\theta}_i\right)\).

3.2 Approximate Inference with DPVI↩︎

We build our noise-aware approximate Bayesian inference method on top of DPVI. Recall that DPVI is based on optimizing the ELBO with DP-SGD. The final posterior approximation returned by DPVI is the \(q(\pmb{\theta}; \pmb{\phi}_T)\), where \(\pmb{\phi}_T\) denote the last parameters returned by the DP-SGD optimization. The DPVI treats these parameters as the true optima, and hence disregards all the stochasticity that was introduces by the Gaussian perturbation. Therefore, the DPVI as such is completely unaware of the noise.

However, due to the privacy accounting of DP-SGD, we do not need to limit our considerations to the last iterate. In fact, we can use the full parameter and noisy gradient traces, respectively \(\mathcal{T} = \{ \pmb{\phi}_{t} \}_{t=1}^{T}\) and \(\widetilde{\mathcal{G}} = \{ \widetilde{\pmb{g}}_{t} \}_{t=1}^T\), for arbitrary post-processing. In the next Section, we will introduce our post-processing model, which allows us to model the DP-induced noise, and make the approximate inference noise-aware.

3.3 Post-processing Model↩︎

Assume for now that our gradients norms are bounded by \(C\). Hence \(\text{clip}(\nabla_{\pmb{\phi}} \ell (\pmb{\phi}_t; x_i), C) = \nabla_{\pmb{\phi}} \ell (\pmb{\phi}_t; x_i)\). From the DP-SGD update equations, we have that the noisy gradient at iteration \(t\) is \[\begin{align} \label{eq:dp-sgd-grad} \widetilde{\pmb{g}}_{t + 1} = \sum_{i \in \mathcal{B}_{t + 1}} \nabla_{\pmb{\phi}} \ell (\pmb{\phi}_t; x_i) + \sigma_{\text{DP}} C \boldsymbol{\eta}, \end{align}\tag{8}\]

where \(\pmb{\eta} \sim \mathcal{N}\left(0, \mathbf{I}_{d}\right)\) with \(\mathbf{I}_{d}\) denoting a \(d\)-dimensional identity matrix. Next, we assume that the subsampled gradients could be approximated with a Gaussian distribution through the central limit theorem. This is possible since, in DP-SGD, a large sampling rate or batch size is usually used (e.g., \(\kappa = 0.1\)). Given the full-data gradient at iteration \(t\), \(\nabla \mathcal{L}(\pmb{\phi}_t; \mathbf{D})\), we have \[\begin{align} \label{eq:subsampled-true-grad} \sum_{i \in \mathcal{B}_{t + 1}} \nabla_{\pmb{\phi}} \ell (\pmb{\phi}_t; x_i) \sim \mathcal{N}\left( \kappa \nabla \mathcal{L}(\pmb{\phi}_t; \mathbf{D}), \Sigma_{\text{sub}}(\pmb{\phi}_t) \right). \end{align}\tag{9}\] Combining 8 9 yields \[\begin{align} \widetilde{\pmb{g}}_{t + 1} \sim \mathcal{N}(\kappa \nabla \mathcal{L}(\pmb{\phi}_t; \mathbf{D}), \sigma_{\text{DP}}^2 C^2 \mathbf{I}_d + \Sigma_{\text{sub}}(\pmb{\phi}_t)). \end{align}\] Since the true gradient, as a data-dependent quantity is clearly unknown, we need to approximate it. We will approximate the true gradient as a linear function parameterized with a matrix \(\mathbf{A}\) and a vector \(\pmb{\phi}^{*}\) through the second-order Taylor approximation of \(\mathcal{L}\): \[\begin{align} \nabla \mathcal{L}(\pmb{\phi}_t; \mathbf{D}) \approx \mathbf{A} (\pmb{\phi}_t - \pmb{\phi}^{*}). \end{align}\] This parametrization is convenient, as it allows us to interpret the \(\pmb{\phi}^{*}\) as the true optimum for the VI problem. If we further assume that \(\Sigma_{\text{sub}}\) is approximately constant around \(\pmb{\phi}^{*}\), then denoting the matrix \(\Sigma_{\text{total}} = \sigma_{\text{DP}}^2 C^2 \mathbf{I}_d + \Sigma_{\text{sub}}\), we can write \[\begin{align} \label{gradient-model-trace-approximate} \begin{aligned} &\widetilde{\pmb{g}}_{t + 1} \mid \pmb{\phi}_{t}, \mathbf{A}, \pmb{\phi}^{*}, \Sigma_{\text{sub}} \sim \mathcal{N}\left(\kappa \mathbf{A} \left(\pmb{\phi}_t - \pmb{\phi}^{*}\right), \Sigma_{\text{total}}\right). \end{aligned} \end{align}\tag{10}\] 2 shows an example of this model based on the exponential distribution [expfam-1], with the noisy gradients \(\widetilde{\pmb{g}}_t\) and a MAP estimate for \(\pmb{\phi}^{*}\) together with the fitted gradient. We now have a model for the noisy gradients, and can perform Bayesian inference on the unknown variables \(\mathbf{A}, \pmb{\phi}^{*}\) and \(\Sigma_{\text{sub}}\). Note that the parameters \(\mathbf{A}\) and \(\pmb{\phi}^{*}\) encapsulate all the information about the latent sensitive data \(\mathbf{D}\). This avoids the costly marginalization of individual samples over the latent data, and instead our models scales for arbitrary number of data samples.

After obtaining the posterior for \(\mathbf{A}, \pmb{\phi}^{*}\) and \(\Sigma_{\text{sub}}\), we can form the final noise-aware posterior approximation. We start by marginalizing out \(\mathbf{A}\) and \(\Sigma_{\text{sub}}\) from the post-processing model, \[\label{post-processing-model-marginalized} p\left(\pmb{\phi}^{*}\mid \mathcal{T}\right) = \int p\left(\mathbf{A}, \pmb{\phi}^{*}, \Sigma_{\text{sub}} \mid \mathcal{T}\right) \mathrm{d}\mathbf{A} \mathrm{d}\Sigma_{\text{sub}}.\tag{11}\] Then by observing that the joint distribution of \(\pmb{\theta}, \pmb{\phi}^{*}\) conditioned on \(\mathcal{T}\) is given by, \[\begin{align} \label{post-processing-joint-model} \begin{aligned} p(\pmb{\theta}, \pmb{\phi}^{*}\mid \mathcal{T}) &= p(\pmb{\theta}\mid \pmb{\phi}^{*}, \mathcal{T}) p(\pmb{\phi}^{*}\mid \mathcal{T}) \\ &= q_{\text{VI}}(\pmb{\theta}; \pmb{\phi}^{*}) p(\pmb{\phi}^{*}\mid \mathcal{T}), \end{aligned} \end{align}\tag{12}\] we obtain the final noise-aware approximate posterior \(\widetilde{p}\left(\pmb{\theta}\mid \mathcal{T}\right)\) by marginalizing out \(\pmb{\phi}^{*}\) in 12 , \[\label{approximate-noise-aware-posterior} \widetilde{p}(\pmb{\theta}\mid \mathcal{T}) = \int q_{\text{VI}}(\pmb{\theta}; \pmb{\phi}^{*}) p(\pmb{\phi}^{*}\mid \mathcal{T}) \mathrm{d}\pmb{\phi}^{*}.\tag{13}\] To provide theoretical justification for our gradient-based model, we first need to state the following assumptions.

1 is a reasonable assumption for the VI problem to be solvable via DP-SGD.

Motivation for the second-order Taylor approximation can be found in 6 which implies that \(\mathbf{A}\) is the Hessian. We can rigorously formalize 2 through the Fréchet derivative of \(\nabla_{\pmb{\phi}} \mathcal{L}\left(\pmb{\phi}; \mathbf{D}\right)\) at \(\pmb{\phi}^{*}\). For some tolerance \(e_{\text{tay}} > 0\), that is problem-dependent, there exists \(1 < T^{*} < T\), and an open ball \(\mathrm{B}_{r^{*}}(\pmb{\phi}^{*})\) around \(\pmb{\phi}^{*}\) with radius \(r^{*}\) that contains \(\pmb{\phi}_{t}\) for all \(t \geq T^{*}\), such that, for all \(\pmb{\phi}\in \mathrm{B}_{r^{*}}(\pmb{\phi}^{*})\) \[\label{frechet-derivative-taylor} \frac{\lVert \nabla_{\pmb{\phi}} \mathcal{L}\left(\pmb{\phi}; \mathbf{D}\right) - \nabla_{\pmb{\phi}}^2 \mathcal{L}\left(\pmb{\phi}^{*};\mathbf{D}\right)(\pmb{\phi}- \pmb{\phi}^{*}) \rVert}{\lVert \pmb{\phi}- \pmb{\phi}^{*}\rVert} < e_{\text{tay}}.\tag{14}\]

Assumption 3 is required to make the clipping operation redundant, allowing us to write the noisy gradient models as in 10 .

We verify this assumption experimentally through the application of our method. We can approximate \(\Sigma_{\text{sub}}\) in principle; however, most of our experiments show that this is not necessary due to the relatively large magnitude of the DP noise, so we can set \(\Sigma_{\text{sub}} = \pmb{0}\). Now we will introduce the following theorem which theoretically justifies our post-processing model.

The proof of 1 is found in 7. If some of these assumptions do not hold in practice, such as 3, then it is sufficient that the approximation given by 10 holds for our method to work. These assumptions mainly aid the theoretical analysis and the derivation of our method. We verify that the approximation in 10 is reasonable in the following section experimentally, by applying the method to various models. Furthermore, we find that the method is sensitive to the learning rate of DP-SGD and the priors. Thus, we analyze this further in 9. We obtain a heuristic for the learning rate based on that analysis, which works well across our experiments.

4 EXPERIMENTS↩︎

4.1 Setup↩︎

In our experiments, we use a diagonal Gaussian distribution as the variational distribution, i.e. \(q_\text{VI}(\pmb{\theta}; \pmb{\phi}) = \mathcal{N}(\pmb{\theta}; \pmb{\mu}, \mathrm{diag}(\pmb{s}))\), where \(\pmb{\mu}\) are the means and \(\pmb{s}\) are the variances. Instead of directly optimizing the variances \(\pmb{s} \in \mathbb{R}_{+}^n\), we parametrize the variational distribution with \(\pmb{u} \in \mathbb{R}^n\) such that \(\pmb{s}\) is an element-wise softplus transformation of \(\pmb{u}\) as \(s_j = \log (1 + e^{u_j})\). Now the variational distribution is parameterized with \(\pmb{\phi}= (\pmb{\mu}, \pmb{u})\) leading to \(d=2n\) variational parameters to optimize in total.

[5] demonstrated that the gradients w.r.t. \(\pmb{u}\) are usually a lot smaller in magnitude than the gradients w.r.t. \(\pmb{\mu}\). Since DP-SGD adds identically distributed noise to each dimension of the gradient, the gradients w.r.t. \(\pmb{u}\) will be affected disproportionally by the noise. To avoid this, we scale up gradients w.r.t. \(\pmb{u}\) before clipping and noise addition, and revert this as a post-processing before we take the update step. For more details about this preconditioning see 13.

We infer the post-processing model in 12 through two methods, the No-U-Turn Sampler (NUTS) [22] and using Laplace’s approximation which we obtain using the Adam optimizer [23]. More details about these methods can be found in 8. Additionally we compare our method against the last iterate DPVI.

We approximate \(\mathbf{A}\) with a diagonal matrix, that is \(\mathbf{A} = \mathrm{diag}(a_1, \ldots, a_d)\). Further, in this section we interpret \(\mathcal{L}\) as the optimization loss minimized by DP-SGD, i.e. the negative ELBO. Hence, around a local optimum, the Hessian is positive-definite, and we constrain \(\mathbf{A}\) to be element-wise positive by parameterizing the post-processing model with \(\pmb{v} \in \mathbb{R}^d\) and mapping the \(v_i\) to \(a_i\) with the softplus transformation, i.e. \(a_i = \text{softplus}(v_i)\). The algorithm blocks for NA-DPVI are provided in 13.7.

We observed empirically that non-informative priors worked poorly for our post-processing model which motivated us to do further analysis (see 9). Instead, we chose the priors for \(\pmb{v}\) and \(\pmb{\phi}^{*}\) based on \(\mathcal{T}\) and \(\widetilde{\mathcal{G}}\) which we discuss in 13.3.

For evaluating the coverages, we apply a slightly modified version of the TARP algorithm ([alg:non-private-ecp]) according to 3.1 which we elaborate on further in 13.4. We repeat TARP \(20\) times in each experiment with \(K=500\) to obtain error estimates for the TARP. We estimate the error as the difference between the coverages at each credible level and the perfect coverages (\(C(\alpha) - (1 - \alpha)\)). We summarize these coverage errors with a single number using the RMSE (Root Mean Square Error) metric for \(N_{\alpha}\) values of \(\alpha\), \[\label{rmse-definition} \text{RMSE} = \sqrt{\frac{1}{N_{\alpha}} \sum_{i = 1}^{N_{\alpha}} \left(C\left(\alpha_i\right) - \left(1 - \alpha_i\right) \right)^2}.\tag{15}\] Across all the coverage experiments we use \(\delta = 10^{-5}\), \(T = 10^4\), \(N = 5000\), \(\kappa = 0.1\), and we use the NUTS method with \(1000\) warmup iterations and then run it for \(4000\) iterations.

4.2 Exponential Families↩︎

We compare our approach to the method of [6] which only works for exponential families, and with DPVIm [5], to show that our method works well for similar models. We implement the following conjugate models,

1. : \(\pmb{\theta} \sim \text{Gamma}\left(\alpha, \beta\right)\), \(\pmb{x} \mid \pmb{\theta} \sim \text{Exp}(\pmb{\theta})\),

2. : \(\pmb{\theta} \sim \text{Beta}\left(\alpha, \beta\right)\), \(\pmb{x} \mid \pmb{\theta} \sim \text{Ber}(\pmb{\theta})\),

3. : \(\pmb{\theta} \sim \text{Dir}\left(\alpha_1, \alpha_2, \alpha_3\right)\), \(\pmb{x} \mid \pmb{\theta} \sim \text{Cat}(\pmb{\theta})\).

For details on the specific parameter transformations we use for these models, see 13.5. The coverages for both our method and Bernstein & Sheldon’s are shown in 3 including the coverages for the naive baseline \(q_{\text{VI}}(\pmb{\theta}\mid \pmb{\phi}_{T})\). We use \(\epsilon = 0.1\) for all the models. We also show the errors between the coverages and the perfect coverages in 3. 1 summarizes the errors using the RMSE in 15 . From the coverages and coverage errors, we can see that both methods work similarly compared to the naive baseline. Our method performs better for [expfam-1] and [expfam-2]; however, Bernstein & Sheldon’s method works better for [expfam-3]. The marginal coverages and the marginal coverage errors for [expfam-3] are shown in 13.5.

4.3 10D Bayesian linear regression↩︎

Noise-aware Bayesian linear regression was done before by [7]; however their approach does not scale to many dimensions. The baselines we provide are Gibbs-SS-Noisy [7], DPVIm [5], and by using the last iterate \(\pmb{\phi}_T\), i.e. constructing credible regions from the distribution \(q_{\text{VI}}(\pmb{\theta}; \pmb{\phi}_T)\). We apply our method and compare it with the baselines over the following model: \[\begin{align} \begin{aligned} &\left(\theta_i, \theta_{\sigma}^2\right) \sim \text{N-}\Gamma^{-1}\left(0, \frac{1}{4}, 20, \frac{1}{2}\right), \\ &x_i \sim \mathcal{N}\left(0, 1\right), \quad y \mid \pmb{x}, \pmb{\theta} \sim \mathcal{N}\left(\pmb{w}^\top \pmb{x}, \theta_{\sigma}^2 \right), \end{aligned} \end{align}\] where the parameters are \(\pmb{\theta} = \left(\theta_1, \ldots, \theta_{10}, \theta_{11}, \theta_{\sigma}^2\right)\) and \(\text{N-}\Gamma^{-1}\) is the normal-inverse-gamma distribution. We define \(\pmb{w} = (\theta_1, \ldots, \theta_{11})\) and \(\pmb{x} = (x_1, \ldots, x_{10}, 1)\). The bias term is \(\theta_{11}\). For the constrained optimization problem, we only condition \(\theta_{\sigma}^2\) and use the same function as 60 .

We compute the coverages for our method using both NUTS and Laplace’s approximation for \(\epsilon \in \{0.1, 0.3, 1.0 \}\). The coverages for our method are shown in 4 including the coverages for the naive baseline \(q_{\text{VI}}(\pmb{\theta}\mid \pmb{\phi}_{T})\). We also show the errors between the coverages and the perfect coverages in 4. [table:linear-regression-10d] summarizes the errors using the RMSE 15 . From the coverages and coverage errors, we can see that both methods work similarly compared to the naive baseline. However, NUTS performs better compared to Laplace’s approximation.

4.4 UCI Adult Bayesian Logistic Regression↩︎

We perform an experiment on the UCI Adult dataset [24] with the Bayesian logistic regression model with standard normal prior. No baseline method can scale to this problem other than the naive (last iterate) baseline and DPVIm [5]. After learning the posterior for \(\pmb{\theta}\), we test the posterior predictive distribution given as \[\label{noise-aware-posterior-predictive} \widetilde{p}(y = 1 \mid \mathcal{T}, \pmb{x}) = \int p(y = 1 \mid \pmb{\theta}, \pmb{x}) \widetilde{p}(\pmb{\theta}\mid \mathcal{T}) \mathrm{d}\pmb{\theta},\tag{16}\] where we replace the integral with its Monte-Carlo estimator. We do \(20\) repeats for each experiment and we plot the calibration curves for NA-DPVI (NUTS & Laplace) and for the other baselines with \(\epsilon \in \{0.1, 0.3, 1.0 \}\). The calibration curves and calibration errors are found in 5. We also compute the RMSE for the calibration errors in 3, this is the same as the square root of the Brier score [25]. From the calibrations and calibration errors, we can see that both NUTS and Laplace outperform the naive baseline and DPVIm. However, NUTS performs slightly better than Laplace’s approximation. For more details about the experiment, see 13.6.

5 DISCUSSION↩︎

Bayesian inference should be a very natural companion to DP because it naturally deals with uncertain information and should thus be able to easily tolerate the noise added for DP. The fact that noise-aware inference is so difficult reminds us that commonly used general-purpose inference methods such as MCMC are not really Bayesian since they do not incorporate a mechanism for reasoning about the accuracy of their results. In this sense, algorithms required for noise-aware inference share a similarity with probabilistic numerics [26].

We cite DPVIm [5] in our work regarding DPVI, which they extensively investigate; however, the noise-awareness concept in their work is not exactly equivalent to ours which we generalized from [7]. Although our proposed NA-DPVI greatly expands the domain of noise-aware inference, it is still restricted by its reliance on potentially inaccurate VI to capture the data-modelling uncertainty. Additionally, our approximations only work under certain conditions, which may not always hold. An explicit privacy-utility-computation trade-off is very difficult to obtain for NA-DPVI, because there are several factors that are challenging to account for, theoretically speaking. These factors include DP-SGD, VI, the post-processing model, and the final approximate posterior obtained either in closed form or through an approximate inference method (MCMC, Laplace).

Another limitation in this work is that we did not account for the privacy leakage from hyper-parameter selection, which was done manually due to the high computational cost. Currently, many papers that employ DP-SGD ignore the privacy accounting of hyper-parameter tuning. For example, see Section 5.2 of [27]. We applied heuristics to choose some of the hyper-parameters (e.g., the learning rate), and other hyper-parameters were shared among most experiments which reduces the privacy leakage. The state of privacy accounting for hyper-parameter tuning is still in its infancy, see Section 5.4 of [28], and this presents an area for future research and improvement, especially for noise-aware inference.

The search for the ultimate noise-aware inference algorithm that could accurately capture both data-modelling uncertainty as well as uncertainty due to DP inference for arbitrary models remains an important goal for future research.

Acknowledgments↩︎

This work was supported by the Research Council of Finland (Flagship programme: Finnish Center for Artificial Intelligence, FCAI, Grant 356499 and Grant 359111), the Strategic Research Council at the Research Council of Finland (Grant 358247) as well as the European Union (Project 101070617). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the granting authority can be held responsible for them. This work has been performed using resources provided by the CSC – IT Center for Science, Finland (Project 2003275).

6 Second-Order Taylor Approximation of the Loss Function↩︎

For any function \(\mathcal{L}(\pmb{\phi})\), if \(\nabla \mathcal{L}\left(\pmb{\phi}^{*}\right) = 0\), then by using the second-order Taylor approximation around \(\pmb{\phi}^{*}\), we can write \(\mathcal{L}(\pmb{\phi})\) as: \[\begin{align} \label{second-order-taylor-approximation} \begin{aligned} \mathcal{L} \left(\pmb{\phi} \right) &\approx \mathcal{L} \left(\pmb{\phi}^{*} \right) + \nabla\mathcal{L} \left(\pmb{\phi}^{*} \right)^T (\pmb{\phi} - \pmb{\phi}^*) + \frac{1}{2}(\pmb{\phi} - \pmb{\phi}^*)^T\nabla^2 \mathcal{L} \left(\pmb{\phi}^{*} \right)(\pmb{\phi} - \pmb{\phi}^*) \\ &\approx \mathcal{L} \left(\pmb{\phi}^{*} \right) + \frac{1}{2}(\pmb{\phi} - \pmb{\phi}^*)^T\nabla^2 \mathcal{L} \left(\pmb{\phi}^{*} \right)(\pmb{\phi} - \pmb{\phi}^*), \end{aligned} \end{align}\tag{17}\] where \(\nabla^2 \mathcal{L} \left(\pmb{\phi}^{*} \right)\) is the Hessian of \(\mathcal{L}\) at \(\pmb{\phi}^{*}\). Taking the gradient of both sides of (17 ) with respect to \(\pmb{\phi}\), we obtain, \[\nabla\mathcal{L}\left(\pmb{\phi} \right) \approx \nabla^2\mathcal{L}\left(\pmb{\phi}^{*} \right) \left(\pmb{\phi} - \pmb{\phi}^*\right).\]

7 Proof of Theorem 1↩︎

8 Noise-Aware Posterior Approximation↩︎

Deriving the closed form of \(p(\pmb{\phi}^{*}, \mathbf{A} \mid \mathbf{T})\) might not always be feasible so we obtain the following approximations. The first is Laplace’s approximation. Since the Hessian is positive-definite, we want \(\mathbf{A}\) to be positive-definite as well, so we need a bijective differentiable function (diffeomorphism) \(F\), such that \(F(\mathbf{V}) = \mathbf{A}\) where \(\mathbf{V}\) is an unconstrained version of \(\mathbf{A}\). Sometimes, \(F\) is referred to as a conditioning transformation [29]. Therefore, \[p\left(\mathbf{V},\pmb{\phi}^{*}\mid \mathcal{T} \right) = p\left(F(\mathbf{V}),\pmb{\phi}^{*}\mid \mathcal{T} \right) \left\lvert J_{F}\left(\mathbf{V}\right) \right\rvert,\] where \(J_{F}\) is the Jacobian of \(F\). To obtain Laplace’s approximation first we need the MAP (Maximum a posteriori) estimate for both \(\mathbf{V}\) and \(\pmb{\phi}^{*}\), \[\label{laplace-map-estimates} \left(\widehat{\mathbf{V}}_{\text{map}}, \widehat{\pmb{\phi}}^{*}_{\text{map}}\right) = \underset{\left(\pmb{\phi}^{*}, \mathbf{V}\right)}{\text{argmax}}\; \log p\left(\mathbf{V},\pmb{\phi}^{*}\mid \mathcal{T} \right).\tag{18}\] Let \(p_{\text{gauss}}\left(\mathbf{V}, \pmb{\phi}^{*}\right)\) be the distribution \[\label{laplace-gaussian-distribution} \mathcal{N}\left( \left( \text{vec}\left(\widehat{\mathbf{V}}_{\text{map}}\right), \widehat{\pmb{\phi}}^{*}_{\text{map}}\right), \Sigma_{\text{lap}} \right),\tag{19}\] where \(\mathrm{vec}(\cdot)\) is the vector representation of a matrix. We can easily marginalize out \(\widehat{\pmb{\phi}}^{*}_{\text{map}}\) from \(p_{\text{gauss}}\left(\mathbf{V}, \pmb{\phi}^{*}\right)\) to approximate (11 ) since it is a Gaussian distribution. We then use this marginalized distribution to approximate (13 ) and obtain a new distribution that we will denote \(\widetilde{p}_{\text{lap}}\left(\pmb{\theta}\mid \mathcal{T} \right)\). The second method to approximate (13 ) is by using any Markov Chain Monte Carlo (MCMC) method to approximately sample from the post-processing model, \[\label{mcmc-samples} \left(\mathbf{A}_{i}, \pmb{\phi}^{*}_{i}\right) \sim p\left( \mathbf{A}, \pmb{\phi}^{*}\mid \mathcal{T} \right),\; 1 \leq i \leq N_{\text{mc}}.\tag{20}\] We can use these samples to approximate (13 ) using the following mixture model, \[\label{mcmc-mixture-model} \widetilde{p}(\pmb{\theta}\mid \mathcal{T}) \approx \frac{1}{N_{\text{mc}}} \sum_{i = 1}^{N_{\text{mc}}} q_{\text{VI}}\left(\pmb{\theta}; \pmb{\phi}^{*}_i \right) \equiv_{\text{def}} \widetilde{p}_{\text{mc}}(\pmb{\theta}\mid \mathcal{T}).\tag{21}\]

9 Estimation of the Hessian Matrix and the Influence of the Learning Rate↩︎

From our experiments, we noticed that the uncertainty estimates of the approximate noise-aware posterior are sensitive to the estimation of the Hessian matrix \(\mathbf{A}\). Through observation, we found that the accuracy of the estimation of \(\mathbf{A}\) is highly correlated with the choice of the constant learning rate \(\lambda\). This motivated us to conduct further analysis which we present in this section.

There is also an equivalent model to (10 ) which we call parameter-based model of the trace that is useful in this analysis, \[\label{parameter-model-trace-approximate} \pmb{\phi}_{t + 1} \mid \pmb{\phi}_{t}, \mathbf{A}, \pmb{\phi}^{*}\sim \mathcal{N}\left(\pmb{\phi}_{t} - \lambda \kappa\mathbf{A} \left(\pmb{\phi}_t - \pmb{\phi}^{*}\right), \lambda^2 \Sigma_{total}\right),\tag{22}\] where we assume that \(\Sigma_{\text{sub}} = \pmb{0}\).

For Assumption 2 to hold, the part of the trace \(\mathcal{T}\) that is stationary around the optimum \(\pmb{\phi}^{*}\) (after \(t = T^{*}\)) needs to be close enough to \(\pmb{\phi}^{*}\). On the other hand, the trace should have enough variability around \(\pmb{\phi}^{*}\) for the inference to be accurate. Intuitively speaking, we want to sample parameters \(\pmb{\phi}_t\) and gradients \(\pmb{g}_t\) around \(\pmb{\phi}^{*}\) to use that information for the post-processing model; and in particular, to estimate the Hessian that describes the curvature of the loss function around \(\pmb{\phi}^{*}\). For a rigorous analysis, we present the following Theorem with the full proof in Supplement 10.

The analysis in 2 can be done without the assumption that \(\mathbf{A}\) is diagonal; however, this assumption makes the analysis easier and the relation more clear. We also perform our experiments with a diagonal \(\mathbf{A}\) matrix so it makes sense to analyze this particular case. Inequality (?? ) entails that there is a trade-off between \(\max_i \sqrt{\text{Var}\left[ \hat{a}_i \right]}\) and \(r^{*}\) (the radius of the open ball that contains the iterates in \(\mathcal{T}^{*}\)), and this puts a limit on how well we can estimate the Hessian with a non-informative prior even if we assume perfect knowledge about \(\pmb{\phi}^{*}\).

Also, according to (22 ), other than \(\Sigma_{\text{DP}}\), it is clear that the value of \(\lambda\) affects the variance of the trace around \(\pmb{\phi}^{*}\). We show the relation between \(\lambda\) and \(\mathbb{E}\left[ \lVert \pmb{\phi}_{t} - \pmb{\phi}^{*} \rVert^2 \right]\) through the Ornstein–Uhlenbeck (OU) SDE which is a known tool for analyzing variants of SGD [30]–[32]. This is usually done by modeling the noise in the discrete-time SGD process through the Wiener process and then approximating it by a continuous-time process using the Itô integral [33].

What 3 entails is that \[\lambda \sigma_{\text{DP}}^2 C^2 \propto \mathbb{E}\left[ \lVert \pmb{\phi}_{t} - \pmb{\phi}^{*} \rVert^2 \right] < (r^{*})^2,\] so if we make \(\lambda\) large enough, then \(\mathbb{E}\left[ \lVert \pmb{\phi}_{t} - \pmb{\phi}^{*} \rVert^2 \right]\) will increase proportional to \(\lambda\) which will invalidate 2 (the second-order Taylor approximation of the loss function) by contraposition. The full proof of 3 and a further discussion is found in 11.

To find a good learning rate, we analyze the convergence of DP-SGD with some additional assumptions and then find the learning rate that optimizes the convergence bound given in the following theorem. Also, we want to write the optimization problem as a minimization problem. If the optimization problem is initially a maximization problem, then we change the direction of optimization by redefining \(\mathcal{L}\) as \(-\mathcal{L}\).

We obtain a heuristic based on finding \(\lambda\) that optimizes the right-hand side of (?? ) by taking the derivative, setting it to zero, and then solving for \(\lambda\) from which we obtain a quantity proportional to the following, \[\lambda_{\text{heur}} = \frac{\sqrt{2} \lambda_{c}}{ \sigma_{\text{DP}} C \sqrt{T d} },\] where \(\lambda_{c} > 0\) is a hyper-parameter. See supplement 12 for the full proof and other details. For our experiments, we mostly use \(\lambda_{c} = 1\) or \(\lambda_{c} = \sqrt{\frac{d}{2}}\). This heuristic makes it easier to tune the learning-rate because it sets \(\lambda\) to the appropriate scale that is required by convergence results. We show experimentally through our experiments that this is a good choice for the learning rate.

10 Proof of Theorem 2↩︎

As said before, an intuitive explanation why the variance of the trace around \(\pmb{\phi}^{*}\) affects the estimation of the Hessian is that assuming no to little variability, i.e. \(\pmb{\phi}_t \approx \pmb{\phi}^{*}\) for \(t = T^{*}, \ldots, T\), then \(\mathbf{A}\left(\pmb{\phi}_t - \pmb{\phi}^{*}\right) \approx \nabla \mathcal{L}\left(\pmb{\phi}_t; \mathbf{D}\right)\) becomes \(\mathbf{A}\left(\pmb{\phi}^{*}- \pmb{\phi}^{*}\right) \approx \nabla \mathcal{L}\left(\pmb{\phi}^{*}; \mathbf{D}\right) = 0\). This implies that \(\mathbf{A} \times 0 \approx 0\), and thus, in this case, it is impossible to estimate \(\mathbf{A}\) from this equation. The following is the proof of 2.

11 Proof of Theorem 3↩︎

For the second-order Taylor approximation to be accurate, certain values of the learning rate are valid. To show this, we will first approximate the discrete DP-SGD process using a stochastic differential equation (SDE). We provide the following proof for 3.

Neither inequality (41 ) nor inequality (39 ) helps us choose the learning rate since the relation between them and inequality (28 ) is not a logical equivalence and we typically don’t have any information about \(\mathbf{A}\) in advance. The reason why we approximated the DP-SGD process using the OU SDE is that it is easier to work with. That said, it is possible to do a similar analysis with the discrete process but with some additional assumptions. To see this, observe that (31 ) can be written as (when \(t = T - 1\)): \[\pmb{\phi}_{T} = \left(\mathbf{I}_{d} - \lambda \kappa \mathbf{A}\right)\pmb{\phi}_{T - 1} + \lambda \kappa \mathbf{A} \pmb{\phi}^{*}+ \lambda \Sigma_{\text{DP}}^{\frac{1}{2}} \pmb{\eta}_{T},\] and that this could be expanded as: \[\pmb{\phi}_{T} = \left(\mathbf{I}_{d} - \lambda \kappa \mathbf{A}\right)^{T}\pmb{\phi}_{0} + \lambda \kappa\sum_{t = 0}^{T - 1} \left(\mathbf{I}_{d} - \lambda \kappa \mathbf{A}\right)^{t} \mathbf{A} \pmb{\phi}^{*}+ \lambda \sum_{t = 0}^{T - 1} \left(\mathbf{I}_{d} - \lambda \kappa\mathbf{A}\right)^{t}\Sigma_{\text{DP}}^{\frac{1}{2}}\pmb{\eta}_{t + 1}.\] Since \(\mathbf{A}\) is a symmetric matrix and thus is diagonalizable, there is a diagonal matrix \(\mathbf{L}\) and an orthogonal matrix \(\mathbf{U}\) such that \(\mathbf{U} \mathbf{U}^\top = \mathbf{I}_{d}\) and \(\mathbf{A} = \mathbf{U} \mathbf{L} \mathbf{U}^\top\), then \(\mathbf{I}_{d} - \lambda \mathbf{L}\) is also a diagonal matrix, Let \(l\) be the absolute value of the largest entry in the diagonal of \(\mathbf{L}\), if \(\lambda < \frac{1}{l\kappa}\), then \[\label{alternative-discrete-analysis-1} \left(\mathbf{I}_{d} - \lambda\kappa \mathbf{A}\right)^{T} \pmb{\phi}_{0} = \left(\mathbf{I}_{d} - \lambda\kappa \mathbf{U}\mathbf{L}\mathbf{U}^{\top}\right)^{T} \pmb{\phi}_{0} =\left( \mathbf{U}\left(\mathbf{I}_{d} - \lambda\kappa \mathbf{L}\right)\mathbf{U}^{\top}\right)^{T} \pmb{\phi}_{0} = \mathbf{U}\left(\mathbf{I}_{d} - \lambda \kappa \mathbf{L}\right)^{T} \mathbf{U}^\top \pmb{\phi}_{0} \overset{T \rightarrow \infty}{\longrightarrow} \pmb{0}.\tag{42}\] Further assume that the diagonal entries of \(\mathbf{L}\) are all non-zero so that the matrix \(\mathbf{A}\) is invertible, let \(\mathbf{Q} = \mathbf{I}_{d} - \lambda \kappa\mathbf{A}\) then, applying the geometric series sum for matrices, \[\begin{align} \label{alternative-discrete-analysis-2} \begin{aligned} \lambda\kappa \sum_{t = 0}^{T - 1} \mathbf{Q}^{t} \mathbf{A} \pmb{\phi}^{*}&=\lambda\kappa\left(\sum_{t = 0}^{T - 1} \mathbf{Q}^{t} \right) \mathbf{A} \pmb{\phi}^{*}\\ &= \lambda\kappa\left(\left(\mathbf{I}_{d} - \mathbf{Q}^{T} \right)\left(\mathbf{I}_{d} - \mathbf{I}_{d} + \lambda \kappa\mathbf{A}\right)^{-1} \right) \mathbf{A} \pmb{\phi}^{*}\\ &=\lambda\kappa \left(\left(\mathbf{I}_{d} - \mathbf{Q}^{T} \right)\left(\lambda \kappa \mathbf{A}\right)^{-1} \right) \mathbf{A} \pmb{\phi}^{*}\\ &= \frac{\lambda\kappa}{\lambda\kappa} \left(\mathbf{I}_{d} - \mathbf{Q}^{T} \right) \pmb{\phi}^{*}\\ &= \left(\mathbf{I}_{d}- \left(\mathbf{I}_{d} - \lambda\kappa\mathbf{U} \mathbf{L}\mathbf{U}^{\top}\right)^{T} \right) \pmb{\phi}^{*}\\ &= \left(\mathbf{I}_{d}- \left(\mathbf{U}\left(\mathbf{I}_{d} - \lambda\kappa \mathbf{L}\right)\mathbf{U}^{\top}\right)^{T} \right) \pmb{\phi}^{*}\\ &= \left(\mathbf{I}_{d}- \mathbf{U} \left(\mathbf{I}_{d} - \lambda\kappa\mathbf{L}\right)^{T} \mathbf{U}^{\top} \right) \pmb{\phi}^{*}\\ &= \mathbf{U}\left(\mathbf{I}_{d}- \left(\mathbf{I}_{d} - \lambda\kappa \mathbf{L}\right)^{T} \right) \mathbf{U}^\top \pmb{\phi}^{*}\overset{T \rightarrow \infty}{\longrightarrow} \pmb{\phi}^{*}. \end{aligned} \end{align}\tag{43}\] Since \(\Sigma_{\text{DP}}\) is \(\sigma_{\text{DP}}^2 C^2 \mathbf{I}_d\), it commutes with \(\mathbf{Q}^{t}\), thus from the fact that the noises \(\pmb{\eta}_{t}\) are i.i.d., \[\begin{align} \label{alternative-discrete-analysis-3} \begin{aligned} \lambda \sum_{t = 0}^{T - 1} \mathbf{Q}^{t}\Sigma_{\text{DP}}^{\frac{1}{2}}\pmb{\eta}_{t + 1} \Bigm| \pmb{\phi}_{0} &\sim \mathcal{N}\left(\pmb{0}, \lambda^2\sigma_{\text{DP}}^2 C^2 \sum_{t = 0}^{T - 1} \mathbf{Q}^{2t}\right) \\ &\sim \mathcal{N}\left(\pmb{0}, \lambda^2\sigma_{\text{DP}}^2 C^2 \left(\mathbf{I}_{d} - \mathbf{Q}^{2T} \right)\left(\mathbf{I}_{d} - \mathbf{Q}^2\right)^{-1} \right). \end{aligned} \end{align}\tag{44}\] Hence, from (42 , (43 ), and (44 ), \[\mathbb{E}\left[ \pmb{\phi}_{T}\right] \overset{T \rightarrow \infty}{\longrightarrow} \pmb{\phi}^{*},\] \[\text{Cov}\left[ \pmb{\phi}_{T}\right] \overset{T \rightarrow \infty}{\longrightarrow} \lambda^2 \sigma_{\text{DP}}^2 C^2 \left(\mathbf{I}_{d} - \mathbf{Q}^2\right)^{-1}.\] Therefore, \[\begin{align} \mathbb{E}\left[ \lVert \pmb{\phi}_{T} - \pmb{\phi}^{*}\rVert^2\right] \overset{T \rightarrow \infty}{\longrightarrow} & \lambda^2 \sigma_{\text{DP}}^2 C^2 \text{Tr}\left(\mathbf{U}^{\top} \left(\mathbf{I}_{d} - \left(\mathbf{I}_{d} - \lambda \kappa\mathbf{L}\right)^2\right)^{-1}\mathbf{U} \right) \\ &= \lambda^2 \sigma_{\text{DP}}^2 C^2 \text{Tr}\left(\left(\mathbf{I}_{d} - \left(\mathbf{I}_{d} - \lambda \kappa\mathbf{L}\right)^2 \right)^{-1}\right) \\ &= \lambda^2 \sigma_{\text{DP}}^2 C^2 \sum_{i = 1}^{d} \frac{1}{1 - (1 - \lambda\kappa l_i)^2} \\ &= \lambda^2 \sigma_{\text{DP}}^2 C^2 \sum_{i = 1}^{d} \frac{1}{1 - 1 + 2 \lambda \kappa l_i - \lambda^2 \kappa^2 l_i^2} \\ &= \frac{\lambda \sigma_{\text{DP}}^2 C^2}{\kappa} \sum_{i = 1}^{d} \frac{1}{2 l_i - \lambda l_i^2} \end{align}\] where \(l_i\) (\(1 \leq i \leq d\)) are the diagonal entries of \(\mathbf{L}\).

12 Proof of Theorem 4↩︎

We state the additional assumptions of 4 again,

Since the per-example losses are assumed to be Lipschitz continuous, then \[\lVert\nabla_{\pmb{\phi}} \mathcal{L}\left(\pmb{\phi}; \mathbf{D}\right) \rVert \leq \sum_{i = 1}^{N} \lVert \nabla_{\pmb{\phi}} \ell\left(\pmb{\phi}; \pmb{x}_{i} \right)\rVert \leq N C,\] so \(\mathcal{L}\) is Lipschitz continuous, also \[\sup_{\pmb{\phi}} \lVert\nabla_{\pmb{\phi}} \mathcal{L}\left(\pmb{\phi}; \mathbf{D}\right) \rVert = G \leq N C.\]

We will denote \(\mathcal{L}\left(\pmb{\phi}; \mathbf{D}\right)\) as \(\mathcal{L}\left(\pmb{\phi}\right)\) for short. Also, we want to write the optimization problem as a minimization problem. If the optimization problem is initially a maximization problem, then we change the direction of optimization by redefining \(\mathcal{L}\) as \(-\mathcal{L}\). The following lemma is a well-known result; however we prove it within this context.

Now the following is the proof of 4:

Define \(\lambda_{\text{opt}}\) as: \[\lambda_{\text{opt}} \equiv_{\text{def}} \frac{\sqrt{2 \left(\mathcal{L}\left(\pmb{\phi}_{0}\right) - M\right)}}{\sqrt{H T \left( 2e_{\text{sub}}^2 + 2\kappa^2 G^2 + \sigma_{\text{DP}}^2 C^2 d \right)}}.\] If \(e_{\text{sub}}^2\) is insignificant compared to \(\kappa^2 G^2 + \sigma_{\text{DP}}^2 C^2 d\), then we can clearly see that \[\label{lambda-optimal-exact} \lambda_{\text{opt}} \propto \frac{\sqrt{2}}{\sqrt{H T \left( \kappa^2 G^2 + \sigma_{\text{DP}}^2 C^2 d \right)}}.\tag{54}\] We choose our \(\lambda\) according to (54 ); however, if we don’t know about \(H\) and \(G\), then we can either estimate them or simply choose \(\lambda_{\text{heur}}\) such that: \[\label{lambda-heuristic} \lambda_{\text{heur}} = \frac{\sqrt{2} \lambda_c}{ \sigma_{\text{DP}} C \sqrt{T d} },\tag{55}\] where \(\lambda_c > 0\) is a hyper-parameter.

13 Experiments Additional Details↩︎

13.1 Per-Example Loss Function Approximation↩︎

The per-example losses based on the ELBO might be difficult to compute in closed form for some problems, so generally we approximate \(\ell\) by sampling from the variational distribution \(q_{\text{VI}}\) and then replacing the expected values with averages. So let \(\left\{\pmb{\theta}_{i}\right\}_{i = 1}^{N_{\text{VI}}}\), \(N_{\text{VI}} > 0\) be samples from \(q_{\text{VI}}\), then \[\begin{align} \label{approximate-per-example-loss} \ell\left(\pmb{\phi}; \pmb{x} \right) \approx \frac{1}{N_{\text{VI}}}\sum_{i = 1}^{N_{\text{VI}}} \left[\log p\left(\pmb{x} \mid \pmb{\theta}_{i} \right)\right] - \frac{1}{N_{\text{VI}} N} \sum_{i = 1}^{N_{\text{VI}}}\left[\log q_{\text{VI}}(\pmb{\theta}_i \mid \pmb{\phi}) - \log p\left(\pmb{\theta}_i \right)\right]. \end{align}\tag{56}\] Throughout our experiments we use \(N_{\text{VI}} = 10\).

If \(p\left(\mathbf{D} \mid \pmb{\theta}_{\text{con}}\right)\) and \(p(\pmb{\theta}_{\text{con}})\) require some constraints on \(\pmb{\theta}\), then we would need to transform \(\pmb{\theta}_{\text{con}}\) to an unconstrained domain to work with. From [34], we can define any diffeomorphism \(\mathcal{U}: \Theta \rightarrow \mathbb{R}^{n}\) such that it transforms \(\pmb{\theta}_{\text{con}}\) from the constrained domain to the unconstrained domain \(\mathbb{R}^{n}\). Because \(\pmb{\theta}\) is unconstrained with respect to the variational distribution, this requires an additional adjustment to (56 ), \[\begin{align} \label{approximate-per-example-loss-constrained} \begin{aligned} \ell\left(\pmb{\phi}; \pmb{x} \right) \approx&\frac{1}{N_{\text{VI}}}\sum_{i = 1}^{N_{\text{VI}}} \left[\log p\left(\pmb{x} \mid \mathcal{U}^{-1}\left(\pmb{\theta}_{i}\right) \right)\right] \\ &- \frac{1}{N_{\text{VI}} N} \sum_{i = 1}^{N_{\text{VI}}}\left[\log q_{\text{VI}}(\pmb{\theta}_i \mid \pmb{\phi}) - \log p\left(\mathcal{U}^{-1}\left(\pmb{\theta}_i\right) \right)\right] \\ & + \frac{1}{N_{\text{VI}} N} \sum_{i = 1}^{N_{\text{VI}}} \left\lvert \det J_{\mathcal{U}^{-1}} \left(\pmb{\theta}_{i}\right)\right \rvert. \end{aligned} \end{align}\tag{57}\] where \(J_{\mathcal{U}^{-1}}\left(\pmb{\theta}_i\right)\) is the Jacobian matrix of \(\mathcal{U}^{-1}\) evaluated at \(\pmb{\theta}_{i}\) and \(\det\) denotes the determinant.

The estimation of \(\pmb{\phi}^{*}\) and \(\mathbf{A}\) is affected by the choice of the learning-rate, especially the estimation of \(\mathbf{A}\) which we fully expand on in 9. We also establish a heuristic for the learning-rate in 10 which we use extensively in our experiments, \[\lambda_{\text{heur}} = \frac{\sqrt{2} \lambda_{c}}{ \sigma_{\text{DP}} C \sqrt{T d} }.\] It makes it easier to get the learning rate to the right scale and the only decision to be made is the choice of \(\lambda_c\) which is a hyper-parameter. Usually, the default value (\(\lambda_c = 1\)) yields good results. Another value for \(\lambda_c\) that we found to yield good results for other models is \(\lambda_c = \sqrt{\frac{d}{2}}\).

13.2 Gradients Preconditioning↩︎

We apply a simple preconditioning technique by modifying (5 ) so that for a vector \(\pmb{\beta} = \left(\pmb{\beta}_{\pmb{\mu}}, \pmb{\beta}_{\pmb{u}}\right)\) which is the concatenation of two vectors \(\pmb{\beta}_{\pmb{\mu}}\) and \(\pmb{\beta}_{\pmb{u}}\) each of dimension \(n\), the update rule becomes, \[\begin{align} \label{modified-dp-sgd-align} \begin{aligned} \pmb{g}_{t + 1} &= \sum_{i \in \mathcal{B}_{t + 1}} \mathrm{clip}\left(\pmb{\beta} \odot \nabla_{\pmb{\phi}}\ell(\pmb{\phi}_t; \pmb{x}_i), C\right), \\ \widetilde{\pmb{g}}_{t + 1} &= \frac{1}{\pmb{\beta}} \odot \left[ \pmb{g}_{t + 1} + \sigma_{\text{DP}} C \pmb{\eta}_{t + 1} \right], \\ \pmb{\phi}_{t + 1} &= \pmb{\phi}_{t} - \lambda \widetilde{\pmb{g}}_{t + 1}, \end{aligned} \end{align}\tag{58}\] where \(\odot\) is the element-wise multiplication operation, and \(\frac{1}{\pmb{\beta}}\) denotes the vector obtained from the reciprocals of the components of \(\pmb{\beta}\). In our experiments, \(\pmb{\beta}_{\pmb{\mu}}\) is filled with \(1\)s and we only choose specific values for \(\pmb{\beta}_{\pmb{u}}\). We also need to modify the gradient-based (10 ) model so that \[\label{modified-gradient-based-model} \widetilde{\pmb{g}}_{t + 1} \mid \pmb{\phi}_{t}, \mathbf{A}, \pmb{\phi}^{*}\sim \mathcal{N}\left(\kappa \mathbf{A} \left(\pmb{\phi}_t - \pmb{\phi}^{*}\right), \frac{1}{\pmb{\beta}} \odot \left( \sigma_{\text{DP}}^2 C^2 \mathbf{I}_d + \Sigma_{\text{sub}} \right)\right),\tag{59}\] where \(\frac{1}{\pmb{\beta}}\odot\Sigma_{\text{DP}}\) is the element-wise multiplication of each row-vector of \(\Sigma_{\text{DP}}\) by \(\frac{1}{\pmb{\beta}}\). Regarding the learning rate \(\lambda\), we use our heuristic but multiply the heuristic with \(\pmb{\beta}\) to ensure also that all the parameters converge at a similar rate, \[\lambda = \lambda_{\text{heur}} \pmb{\beta},\] and this makes lambda into a vector so that the product in the last equation in (58 ), i.e. \(\lambda \widetilde{g}_{t + 1}\) becomes an element-wise product \(\lambda \odot \widetilde{g}_{t + 1}\).

13.3 Setting the Priors↩︎

We choose Gaussian distributions as priors for both \(\pmb{\phi}^{*}\) and each \(v_i\). For \(\pmb{\phi}^{*}\) we set, \[\pmb{\phi}^{*}\sim \mathcal{N}\left(\frac{1}{T - T^{*}} \sum_{t = T^{*}}^{T - 1}\pmb{\phi}_{t}, \mathbf{I}_{d}\right).\] The prior for each \(v_i\) is set based on the MLE estimate for \(\mathbf{A} = \mathrm{diag}(v_1, \ldots, v_d)\) in the proof of 2, \[\begin{align} v_i &\sim \mathcal{N}\left(\mu_{v_i}, \sigma_{v_i}^2\right), \\ \mu_{v_i} &= \frac{\left\lvert \sum_{t=T^{*}}^{T - 1} \widetilde{\pmb{g}}_{t + 1}^{(i)} \left(\pmb{\phi}_t^{(i)} - \overline{\pmb{\phi}}^{(i)}\right) \right \rvert}{\kappa \sum_{t=T^{*}}^{T - 1} \left(\pmb{\phi}_t^{(i)} - \overline{\pmb{\phi}}^{(i)}\right)^2}, \\ \sigma_{v_i} &= \frac{\sigma_{\text{DP}}^2 C^2}{\kappa^2 \left(\pmb{\beta}^{(i)}\right)^2 \sum_{t = T^{*}}^{T - 1} \left(\pmb{\phi}_{t}^{(i)} - \overline{\pmb{\phi}}^{(i)} \right)^2}. \end{align}\] where \(\overline{\pmb{\phi}}\) is the average of the trace after the burn-in index \(\pmb{\phi}^{*}\), \[\overline{\pmb{\phi}} = \frac{1}{T - T^{*}} \sum_{t = T^{*}}^{T - 1}\pmb{\phi}_{t}.\]

13.4 TARP Evaluation Method↩︎

We can further approximate the average coverages: \[C(\alpha) = \frac{1}{K}\sum_{i = 1}^{K} \text{\usefont{U}{bbold}{m}{n}1}_{i, \alpha}\left(\pmb{\theta}_i\right),\] by applying the TARP algorithm from [14]. This is done by sampling from the approximate posterior \(\widetilde{p}\left(\pmb{\theta}\mid \pmb{\xi}\right)\) and assuming that the credible regions are balls centered around \(\pmb{\theta}_{\text{ref}}\left(\pmb{\xi}_i\right)\). Denote \(\pmb{\theta}_{\text{ref}}\left(\pmb{\xi}_i\right)\) as \(\pmb{\theta}_{\text{ref}(i)}\) and let \(\left\{ \widetilde{\pmb{\theta}}_{j, i} : 1 \leq j \leq N_{\text{tarp}}\right\}\) be samples from the approximate posterior \(\widetilde{p}\left(\pmb{\theta}\mid \pmb{\xi}_{i} \right)\) for each \(i\), then \[C(\alpha) \approx \frac{1}{K}\sum_{i = 1}^{K} \text{\usefont{U}{bbold}{m}{n}1}\left[f_i < 1 - \alpha\right],\] where \[f_i = \frac{1}{N_{\text{tarp}}} \sum_{j = 1}^{N_{\text{tarp}}} \text{\usefont{U}{bbold}{m}{n}1}\left[d\left( \widetilde{\pmb{\theta}}_{j, i}, \pmb{\theta}_{\text{ref}(i)}\right) < d\left( \pmb{\theta}_{i}, \pmb{\theta}_{\text{ref}(i)}\right) \right].\] The motivation for why this is a valid approximation is that: \[\pmb{\theta}_i \in \widetilde{\mathcal{C}}_{\alpha}\left(\pmb{\xi}_i, \pmb{\theta}_{\text{ref}(i)}\right) \iff d\left( \pmb{\theta}_{i}, \pmb{\theta}_{\text{ref}(i)}\right) < r\left(\alpha, \pmb{\xi}_i\right),\] where \(r\left(\alpha, \pmb{\xi}_i\right)\) is the radius of the credible region. Also, this is equivalent to \[\text{Ball}\left(\pmb{\theta}_{\text{ref}(i)}, d\left( \pmb{\theta}_{i}, \pmb{\theta}_{\text{ref}(i)}\right)\right) \subset \widetilde{\mathcal{C}}_{\alpha}\left(\pmb{\xi}_i, \pmb{\theta}_{\text{ref}(i)}\right),\] i.e. the ball centered around \(\pmb{\theta}_{\text{ref}(i)}\) with radius \(d\left( \pmb{\theta}_{i}, \pmb{\theta}_{\text{ref}(i)}\right)\) is contained in the credible region. Therefore, it is also equivalent to \[\int_{\Theta} \text{\usefont{U}{bbold}{m}{n}1}\left[ \widetilde{\pmb{\theta}} \in \text{Ball}\left(\pmb{\theta}_{\text{ref}(i)}, d\left( \pmb{\theta}_{i}, \pmb{\theta}_{\text{ref}(i)}\right)\right)\right] \widetilde{p}\left(\widetilde{\pmb{\theta}} \mid \pmb{\xi}_{i} \right) < 1 - \alpha.\] This integral can be approximated by sampling from \(\widetilde{p}\left(\widetilde{\pmb{\theta}} \mid \pmb{\xi}_{i} \right)\) and then using these samples to calculate the average of \(\text{\usefont{U}{bbold}{m}{n}1}\left[ \widetilde{\pmb{\theta}} \in \text{Ball}\left(\pmb{\theta}_{\text{ref}(i)}, d\left( \pmb{\theta}_{i}, \pmb{\theta}_{\text{ref}(i)}\right)\right)\right]\). Also, \[\begin{align} & \widetilde{\pmb{\theta}} \in \text{Ball}\left(\pmb{\theta}_{\text{ref}(i)}, d\left( \pmb{\theta}_{i}, \pmb{\theta}_{\text{ref}(i)}\right)\right) \\ \iff & \text{\usefont{U}{bbold}{m}{n}1}\left[ d\left(\widetilde{\pmb{\theta}} ,\pmb{\theta}_{\text{ref}(i)} \right) < d\left( \pmb{\theta}_{i}, \pmb{\theta}_{\text{ref}(i)}\right) \right], \end{align}\] and thus the approximation holds. One important thing to consider here, is that if \(\pmb{\theta} \in \Theta\) is constrained, we must have a diffeomorphism \(\mathcal{U}: \Theta \rightarrow \mathbb{R}^{n}\) that transforms the prior samples \(\pmb{\theta}_{i}\) to the unconstrained space \(\mathcal{U}\left(\pmb{\theta}_{i}\right)\). Since \(\mathcal{U}\) is a diffeomorphism and thus both injective and surjective, it must preserve the relation \(\in\). Hence, \[\pmb{\theta} \in X \iff \mathcal{U}\left(\pmb{\theta}\right) \in \mathcal{U}\left(X\right)\] similarly if \(\pmb{\theta}' \in \mathbb{R}^{n}\), then \[\pmb{\theta}' \in Y \iff \mathcal{U}^{-1}\left(\pmb{\theta}'\right) \in \mathcal{U}^{-1}\left(Y\right)\]

According to Theorem 3 of [14], if \(\widetilde{p}\left(\pmb{\theta} \mid \mathcal{T} \right)\) had the correct coverages, then the marginals \(\widetilde{p}\left(\pmb{\theta}^{(i)} \mid \mathcal{T} \right)\) should also have correct coverages. We can also use the TARP algorithm to calculate the coverages for each dimension of the parameters vector \(\pmb{\theta}\), i.e., the marginal coverages \(\widetilde{p}\left(\pmb{\theta}^{(i)} \mid \mathcal{T}\right)\) where \(\pmb{v}^{(i)}\) denotes the \(i\)th component of a vector \(\pmb{v}\).

13.5 Exponential Families Experiment Details↩︎

For the constrained optimization problem of [expfam-1], \(\pmb{\theta}_{\text{con}} > 0\), so we define \[\label{softplus-conditioning} \mathcal{U}_1\left(\pmb{\theta}_{\text{con}}\right) = \log(e^{\pmb{\theta}_{\text{con}}} - 1) \implies \mathcal{U}_1^{-1}\left(\pmb{\theta}\right) = \text{softplus}\left(\pmb{\theta}\right).\tag{60}\] For Model [expfam-2], \(\pmb{\theta}_{\text{con}} \in (0, 1)\), so we define \[\mathcal{U}_2\left(\pmb{\theta}_{\text{con}}\right) = \log \left(\frac{\pmb{\theta}_{\text{con}}}{1 - \pmb{\theta}_{\text{con}}}\right) \implies \mathcal{U}_2^{-1}\left(\pmb{\theta}\right) = \frac{1}{1 + e^{-\pmb{\theta}}}.\] Finally, for Model [expfam-3], let \(\pmb{\theta}_{\text{con}} = \left(\theta_1, \theta_2, \theta_3\right)\), then \(\theta_1 + \theta_2 + \theta_3 = 1\) so we only have two degrees of freedom and can let \(\theta_3 = 1 - \theta_1 - \theta_2\). We define, \[\mathcal{U}_3\left(\pmb{\theta}_{\text{con}}\right) = \left(\log\left(\frac{\theta_1}{\theta_3}\right), \log\left(\frac{\theta_2}{\theta_3}\right), 0\right).\] For the inverse, let \(\pmb{\theta} = \left(\theta_1', \theta_2'\right)\), then \[\mathcal{U}_3^{-1}\left(\pmb{\theta}\right) = \text{softmax}\left(\theta_1', \theta_2', 0\right).\]

13.6 UCI Adult Experiment Details↩︎

For data pre-processing, we removed the columns “education-num”, “native-country”, and “relationship”, and converted categorical values to one-hot encoded vectors. The continuous values were normalized to be within \((0, 1)\). Moreover, we evaluated our method for \(\epsilon \in \{0.1, 0.3, 1.0\}\). Across these different values of \(\epsilon\), we used the same number of iterations \(T = 10^{4}\) and a sampling rate of \(\kappa = 0.1\).

13.7 NA-DPVI Algorithm↩︎

[SUBALG]IndentEndIndent

The variables that are used in all the algorithms are found in 6. The DPVI algorithm that we used can be found in 7 from which the parameter and gradient traces (\(\mathcal{T}\), \(\widetilde{\mathcal{G}}\)) are obtained. Finally, our NA-DPVI algorithm can be found in 8.

The time complexity of DPVI (7) is \(\mathcal{O}(T\times B \times d)\) where \(T\) is the number of steps of DP-SGD, \(B = \kappa N\) is the expected batch size (\(\kappa\) is the subsampling ratio and \(N\) is the number of samples), and \(d\) is dimensionality of \(\pmb{\phi}\) (i.e. parameters).

The time complexity of NA-DPVI (8) has an additional cost of sampling \(M\) from the posterior \(p\left(\pmb{\phi}^{*}, \mathbf{A}, \Sigma_{\text{sub}} \mid \widetilde{\mathcal{G}}, \mathcal{T}\right)\), which depends on the method of approximating the noise-aware posterior 13 (e.g., NUTS and Laplace). The complexity of the noise-aware posterior approximation method is a function \(M\), \(d\), and \(T\) since the approximation method takes as its input the trace \(\mathcal{T}\) and perturbed gradients \(\widetilde{\mathcal{G}}\).

13.8 Computational Resources↩︎

We used the computational resources offered by CSC – IT Center for Science, Finland. In particular, the Puhti supercomputer was used to run experiments on CPU nodes in parallel. Each Puhti node is equipped with two Intel Xeon processors, code name Cascade Lake, with 20 cores each running at 2.1 GHz.

References↩︎