解析高斯机制（Analytic Gaussian Mechanism）

解析高斯机制（Analytic Gaussian Mechanism）

本文不保证正确性，任何错误都是有可能的

---

前言

需要注意，Analytic Gaussian Mechanism在文献[1]（没看过原文）中被提出，但本文对Analytic Gaussian Mechanism定义和理解来自于文献[2]。
[1] Borja Balle and Yu-Xiang Wang. 2018. Improving the Gaussian Mechanism for Diferential Privacy: Analytical Calibration and Optimal Denoising. In ICML. 403–412.
[2] Ergute Bao, Yin Yang, Xiaokui Xiao, and Bolin Ding. 2021. CGM: an enhanced mechanism for streaming data collection with local differential privacy. Proc. VLDB Endow. 14, 11 (July 2021), 2258–2270.

Github上有关于Analytic Gaussian Mechanism的实现代码（https://github.com/BorjaBalle/analytic-gaussian-mechanism/blob/master/agm-example.py），但是和文献[2]的定义对不上，所以这里就没用上，读者可自行参考。

---

一、定义

Analytic Gaussian Mechanism是对Gaussian Mechanism的一种改进，其本质依旧是向数据中注入满足$(ϵ,\delta )$-差分隐私的高斯噪声，其噪声是从高斯分布$\mathcal{N}(0,{\sigma }^2)$中的取样。高斯机制的更多解释可参考文章【我们为什么用高斯机制？】。

Analytic Gaussian Mechanism的定义如下（文献[2]第三页Lemma2.4）

上述定义中的$S(F)$表示$L_2$灵敏度，对于单个数据而言，灵敏度计算方式为（最大值-最小值）。$\mathbf{I}$表示单位矩阵，因为文中假设需要扰动的数据是$d$-维向量，所以每个数据需要$d$个取样。本文谈论的是单个数据的扰动，即给定参数$ϵ$和$\delta$，计算出注入数据$x$的噪声$\mathcal{N}(0,{\sigma }^2)$。

二、实现

高斯噪声的计算过程如下：
1.根据给定参数$ϵ$和$\delta$，计算$\mathcal{X}$
2.根据参数$\mathcal{X}$，计算$\sigma$
3.根据参数$\sigma$，取样

1.计算$\mathcal{X}$

定义误差函数（Eerror Function, ERF）
$$erf(x)=\frac{2}{\sqrt{\pi }}{\int }_0^{\infty }{e}^{-t^2}dt.$$
则互补误差函数（Error Function Complement, erfc）写为
$$erfc(x)=1-erf(x).$$

根据定义公式，$\mathcal{X}$的计算方式为（注意，$exp(x)=e^x$）
$$\begin{array}{rl}erfc(\mathcal{X})-e^{ϵ}\cdot erfc(\sqrt{{\mathcal{X}}^2+ϵ})& =2\delta \\ 1-erf(\mathcal{X})-e^{ϵ}\cdot (1-erf(\sqrt{{\mathcal{X}}^2+ϵ}))& =2\delta \\ 1-e^{ϵ}+e^{ϵ}\cdot erf(\sqrt{{\mathcal{X}}^2+ϵ})-erf(\mathcal{X})& =2\delta \\ e^{ϵ}\cdot erf(\sqrt{{\mathcal{X}}^2+ϵ})-erf(\mathcal{X})& =2\delta -1+e^{ϵ}\end{array}$$

将误差函数代入得到
$$\begin{array}{rl}{e}^{ϵ}\cdot \frac{2}{\sqrt{\pi }}{\int }_0^{\sqrt{{\mathcal{X}}^2+ϵ}}{e}^{-t^2}dt-\frac{2}{\sqrt{\pi }}{\int }_0^{\mathcal{X}}{e}^{-t^2}dt& =2\delta -1+e^{ϵ}\\ e^{ϵ}\cdot {\int }_0^{\sqrt{{\mathcal{X}}^2+ϵ}}{e}^{-t^2}dt-{\int }_0^{\mathcal{X}}{e}^{-t^2}dt& =\frac{\sqrt{\pi }}{2}\cdot (2\delta -1+e^{ϵ})\end{array}$$

因为${\int }_0^x{e}^{-t^2}dt$ 没有初等函数的解析表达式，我们没有办法直接解出$\mathcal{X}$，转而使用数值方法来解。在$ϵ$和$\delta$已知的情况下，等式右边可以视为一个常数$C$，定义函数$F(x)={\int }_0^x{e}^{-t^2}dt$，则上述式子可以写为
$$e^{ϵ}\cdot F(\sqrt{x^2+ϵ})-F(x)=C$$

然后，定义函数$G(x)=e^{ϵ}\cdot F(\sqrt{x^2+ϵ})-F(x)$，我们需要计算$G(x)=C$。根据 链式法则 和 积分函数的导数等于被积函数 的原则，我们求出$G^{\prime }(x)$为
$$\begin{array}{rl}{G}^{\prime }(x)& =e^{ϵ}\cdot F^{\prime }(\sqrt{x^2+ϵ})-F^{\prime }(x)\\ & =e^{ϵ}\cdot e^{-x^2-1}\cdot \frac{x}{\sqrt{x^2+1}}-e^{2x^2}\end{array}$$

之后，可以使用牛顿迭代法来计算：
1.选取初始值$x_0$
2.计算$G(x_n)$和$G^{\prime }(x_n)$
3.更新$x_{n+1}=x_n-\frac{G(x_n)-C}{G^{\prime }(x_n)}$
4.检查$|G(x_{n+1})-C|$是否收敛（小于选定的阈值），不收敛则返回步骤2

或者直接使用Python中的fsolve函数求解。此外，在计算$G(x)$时还是不可避免的需要计算${\int }_0^x{e}^{-t^2}dt$，详细的解法可以参考【高斯误差函数erf的数值计算方法】。简单的做法是直接使用Python中的quad函数，最后得到的结果$x$可视为$\mathcal{X}$。

2.计算$\sigma$和噪声

有了$\mathcal{X}$，接下来便是计算标准差$\sigma$，定义中给出的条件公式为不等式小于等于，我们直接考虑等于的情况，即
$$\sigma =\frac{S(F)}{\sqrt{2}(\sqrt{{\mathcal{X}}^2+ϵ}-\mathcal{X})}$$
然后便可从$\mathcal{X}(0,{\sigma }^2)$中取样得到噪声。

代码如下（示例）：

import math

import numpy as np
import numpy.random
from scipy.integrate import quad
from scipy.optimize import fsolve

epsilon = 1
delta = 10 ** -5
sensitivity = math.sqrt(2)


# 定义 F(x)
def F(x):
    return quad(lambda t: np.exp(-t**2), 0, x)[0]


# 定义 G(x)
def G(x, C):
    return (math.e ** epsilon) * F(np.sqrt(x**2 + 1)) - F(x) - C


# 设定 C 的值
C = (2 * delta + math.e ** epsilon - 1) * math.sqrt(math.pi) / 2

# 使用 fsolve 求解
X = fsolve(lambda x: G(x, C), np.array([2.0]))[0]

sigma = sensitivity / (math.sqrt(2) * (math.sqrt(X ** 2 + epsilon) - X))

noise = numpy.random.normal(0, sigma ** 2)
print(X, sigma, noise)