Differential Privacy Theory

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Sample from Laplace(0, b) using inverse CDF
5
// If U ~ Uniform(0,1), then X = b * sign(U-0.5) * ln(1 - 2|U-0.5|)
6
// We implement the equivalent branch form for numerical stability.
7
double sample_laplace(double b, mt19937_64 &rng) {
8
    uniform_real_distribution<double> U(0.0, 1.0);
9
    double u = U(rng);
10
    if (u < 0.5) {
11
        return b * log(2.0 * u); // negative value
12
    } else {
13
        return -b * log(2.0 * (1.0 - u)); // positive value
14
    }
15
}
16
17
// Private count query with Laplace noise. Sensitivity Δ1 = 1 for counts.
18
int dp_count_with_laplace(const vector<int> &binary_data, double epsilon, mt19937_64 &rng) {
19
    // Enforce bounds: treat any nonzero as 1. This keeps sensitivity at 1.
20
    int true_count = 0;
21
    for (int x : binary_data) true_count += (x != 0);
22
23
    double b = 1.0 / max(epsilon, 1e-12); // scale = Δ1/ε, with guard
24
    double noise = sample_laplace(b, rng);
25
    // Post-process: round to nearest int and clamp to valid range if desired
26
    int noisy = (int) llround((double)true_count + noise);
27
    noisy = max(0, min((int)binary_data.size(), noisy));
28
    return noisy;
29
}
30
31
int main(){
32
    random_device rd; mt19937_64 rng(rd());
33
34
    // Example dataset of 0/1 participation flags
35
    vector<int> data = {1,0,1,1,1,0,0,1,1,0,1,0,1};
36
37
    double epsilon = 0.5; // stronger privacy when smaller
38
    int noisy_count = dp_count_with_laplace(data, epsilon, rng);
39
40
    cout << "Noisy count (ε=" << epsilon << ") = " << noisy_count << "\n";
41
    return 0;
42
}
43

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Numerically stable sampling from probabilities proportional to exp(alpha * u_i)
5
int sample_exponential_mechanism(const vector<double> &utilities, double alpha, mt19937_64 &rng) {
6
    // alpha = ε / (2Δu). Ensure Δu is known and bounded by design.
7
    double m = *max_element(utilities.begin(), utilities.end());
8
    vector<double> weights;
9
    weights.reserve(utilities.size());
10
    for (double u : utilities) weights.push_back(exp(alpha * (u - m))); // log-sum-exp trick
11
    double Z = accumulate(weights.begin(), weights.end(), 0.0);
12
    uniform_real_distribution<double> U(0.0, Z);
13
    double r = U(rng);
14
    double c = 0.0;
15
    for (int i = 0; i < (int)weights.size(); ++i) {
16
        c += weights[i];
17
        if (r <= c) return i;
18
    }
19
    return (int)weights.size() - 1; // fallback
20
}
21
22
int main(){
23
    random_device rd; mt19937_64 rng(rd());
24
25
    // Suppose we want to privately choose the most popular category among m categories.
26
    // Utility u(D, r) is the count of category r in the dataset (Δu = 1).
27
    vector<int> data = {0,1,2,2,2,1,0,2,0,1,1,2,0,1,2};
28
    int m = 3; // categories 0,1,2
29
30
    vector<double> counts(m, 0.0);
31
    for (int x : data) {
32
        if (0 <= x && x < m) counts[x] += 1.0; // enforce bounds
33
    }
34
35
    double epsilon = 0.8; // privacy budget for this selection
36
    double delta_u = 1.0; // sensitivity of count utility
37
    double alpha = epsilon / (2.0 * delta_u);
38
39
    int chosen = sample_exponential_mechanism(counts, alpha, rng);
40
41
    cout << "Privately chosen category (ε=" << epsilon << ") = " << chosen << "\n";
42
    return 0;
43
}
44

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct Example { double x, y; };
5
6
// Clip a 2D gradient [gw, gb] to L2 norm C
7
array<double,2> clip_grad(array<double,2> g, double C) {
8
    double norm = sqrt(g[0]*g[0] + g[1]*g[1]);
9
    if (norm <= C || norm == 0.0) return g;
10
    double scale = C / norm;
11
    return {g[0]*scale, g[1]*scale};
12
}
13
14
int main(){
15
    random_device rd; mt19937_64 rng(rd());
16
17
    // Tiny synthetic dataset for y ≈ w*x + b
18
    vector<Example> data = {{-2, -3.8}, {-1, -2.0}, {0, 0.2}, {1, 2.1}, {2, 3.9}};
19
20
    // Model parameters
21
    double w = 0.0, b = 0.0;
22
    double lr = 0.1;          // learning rate
23
    double C = 1.0;           // clipping norm
24
    double noise_multiplier = 1.0; // z in σ = z * C
25
26
    // Build a minibatch by subsampling each example with prob q
27
    double q = 0.6; // subsampling rate (amplifies privacy)
28
    bernoulli_distribution take(q);
29
30
    vector<array<double,2>> per_example_grads;
31
    for (auto &ex : data) {
32
        if (!take(rng)) continue;
33
        // Per-example gradient of squared error L = (w*x + b - y)^2
34
        double err = (w * ex.x + b) - ex.y;
35
        array<double,2> g = {2.0 * err * ex.x, 2.0 * err}; // [dL/dw, dL/db]
36
        per_example_grads.push_back(clip_grad(g, C));
37
    }
38
39
    if (per_example_grads.empty()) {
40
        cout << "Empty minibatch; no update this step.\n";
41
        return 0;
42
    }
43
44
    // Aggregate clipped gradients
45
    array<double,2> sumg = {0.0, 0.0};
46
    for (auto &g : per_example_grads) { sumg[0] += g[0]; sumg[1] += g[1]; }
47
48
    // Add Gaussian noise with stddev σ*C to each coordinate (σ = noise_multiplier)
49
    normal_distribution<double> N(0.0, noise_multiplier * C);
50
    sumg[0] += N(rng);
51
    sumg[1] += N(rng);
52
53
    // Average by batch size
54
    double B = (double)per_example_grads.size();
55
    array<double,2> avg = {sumg[0] / B, sumg[1] / B};
56
57
    // Parameter update
58
    w -= lr * avg[0];
59
    b -= lr * avg[1];
60
61
    cout << fixed << setprecision(4);
62
    cout << "After one DP-SGD step: w=" << w << ", b=" << b << " (B=" << B << ")\n";
63
    cout << "Note: Computing final (ε, δ) requires a privacy accountant (e.g., RDP).\n";
64
    return 0;
65
}
66