Differential privacy for aggregates - add dp_laplace and dp_gaussian built-in functions with budget accounting

benchmark/scripts/collect_results.py

@@ -0,0 +1,41 @@
+"""
+Parse per-run accuracy files into a single results.csv.
+
+Output columns: label, epsilon, private, accuracy
+"""
+import pathlib, csv, re
+
+RESULTS = pathlib.Path("benchmark/results")
+
+rows = []
+
+def parse_acc(path: pathlib.Path) -> float:
+    txt = path.read_text().strip()
+    # SystemDS writes a bare float.
+    return float(txt)
+
+# Non-private baseline.
+baseline_path = RESULTS / "acc_baseline.txt"
+if baseline_path.exists():
+    rows.append(dict(label="baseline", epsilon="inf",
+                     private=0, accuracy=parse_acc(baseline_path)))
+
+# DP runs.
+for eps in [0.5, 1, 4, 8]:
+    p = RESULTS / f"acc_eps_{eps}.txt"
+    if p.exists():
+        rows.append(dict(label=f"ε={eps}", epsilon=eps,
+                         private=1, accuracy=parse_acc(p)))
+    else:
+        print(f"Warning: {p} not found — skipping")
+
+out = RESULTS / "results.csv"
+with open(out, "w", newline="") as f:
+    w = csv.DictWriter(f, fieldnames=["label","epsilon","private","accuracy"])
+    w.writeheader()
+    w.writerows(rows)
+
+print(f"Wrote {out}")
+for r in rows:
+    print(f"  {r['label']:12s}  acc={r['accuracy']:.4f}")
+

benchmark/scripts/eval.dml

@@ -0,0 +1,22 @@
+# eval.dml — compute binary classification accuracy on held-out test set.
+# Arguments: data_dir, model_path, out_acc
+data_dir   = $data_dir;
+model_path = $model_path;
+out_acc    = $out_acc;
+
+X_test = read(data_dir + "/X_test.csv",
+              data_type="matrix", value_type="double", format="csv");
+y_test = read(data_dir + "/y_test.csv",
+              data_type="matrix", value_type="double", format="csv");
+w      = read(model_path,
+              data_type="matrix", value_type="double", format="csv");
+
+scores    = X_test %*% w;
+preds     = (scores > 0.0);          # threshold at 0 (log-odds)
+correct   = sum(preds == y_test);
+n_test    = nrow(y_test);
+accuracy  = correct / n_test;
+
+print("Accuracy: " + accuracy);
+write(accuracy, out_acc, format="csv");
+

benchmark/scripts/fedavg_dp.dml

@@ -0,0 +1,130 @@
+# ── fedavg_dp.dml ────────────────────────────────────────────────────────────
+# Arguments (passed via -nvargs):
+#   data_dir   : path to benchmark/data/
+#   epsilon    : DP privacy budget ε  (use 9999 for non-private baseline)
+#   delta      : DP delta             (1e-5 can be used)
+#   clip_norm  : per-example gradient L2-norm clip bound (default 4.0) —
+#                sensitivity = clip_norm / n follows from this, since
+#                clipping each record's gradient contribution to clip_norm
+#                is what makes that bound actually hold.
+#   n_rounds   : number of FedAvg rounds (default 50)
+#   lr         : learning rate          (default 0.1)
+#   w1_rows    : rows in worker 1 shard
+#   w2_rows    : rows in worker 2 shard
+#   w3_rows    : rows in worker 3 shard
+#   w4_rows    : rows in worker 4 shard
+#   n_features : number of features
+#   out        : path to write final weights
+#   private    : 1 = apply DP noise (default), 0 = non-private baseline
+
+data_dir   = $data_dir;
+epsilon    = $epsilon;
+delta      = $delta;
+clip_norm  = ifdef($clip_norm,  4.0);
+n_rounds   = ifdef($n_rounds,   50);
+lr         = ifdef($lr,         0.1);
+private    = ifdef($private,    1);
+out_path   = $out;
+
+# Per-round epsilon: dp_gaussian is called once per round, and
+# DPBudgetAccountant composes the cost of every release against the single
+# total budget set by dp_set_budget(). Spending the full epsilon on every
+# round would exhaust the budget almost immediately, so split it evenly
+# across rounds instead.
+round_epsilon = epsilon / n_rounds;
+
+w1r = $w1_rows;
+w2r = $w2_rows;
+w3r = $w3_rows;
+w4r = $w4_rows;
+d   = $n_features;
+n   = w1r + w2r + w3r + w4r;
+
+# Sensitivity of the released (mean) gradient to a single record changing:
+# with each record's contribution clipped to L2-norm <= clip_norm below,
+# the sum can move by at most clip_norm, so the average moves by clip_norm/n.
+sensitivity = clip_norm / n;
+
+# ── Build federated matrix from 4 local workers ───────────────────────────
+# Row ranges are 0-based [begin, end) for each partition.
+r1s=0;        r1e=w1r;
+r2s=w1r;      r2e=w1r+w2r;
+r3s=w1r+w2r;  r3e=w1r+w2r+w3r;
+r4s=w1r+w2r+w3r; r4e=n;
+
+X = federated(
+    addresses=list(
+        "localhost:8301/" + data_dir + "/worker1/X_train.csv",
+        "localhost:8302/" + data_dir + "/worker2/X_train.csv",
+        "localhost:8303/" + data_dir + "/worker3/X_train.csv",
+        "localhost:8304/" + data_dir + "/worker4/X_train.csv"),
+    ranges=list(
+        list(r1s, 0), list(r1e, d),
+        list(r2s, 0), list(r2e, d),
+        list(r3s, 0), list(r3e, d),
+        list(r4s, 0), list(r4e, d)));
+
+y = federated(
+    addresses=list(
+        "localhost:8301/" + data_dir + "/worker1/y_train.csv",
+        "localhost:8302/" + data_dir + "/worker2/y_train.csv",
+        "localhost:8303/" + data_dir + "/worker3/y_train.csv",
+        "localhost:8304/" + data_dir + "/worker4/y_train.csv"),
+    ranges=list(
+        list(r1s, 0), list(r1e, 1),
+        list(r2s, 0), list(r2e, 1),
+        list(r3s, 0), list(r3e, 1),
+        list(r4s, 0), list(r4e, 1)));
+
+# ── Initialise weights ────────────────────────────────────────────────────
+w = matrix(0, rows=d, cols=1);
+
+if (private == 1) {
+        eps = dp_set_budget($epsilon, $delta);
+}
+
+# ── Training rounds ───────────────────────────────────────────────────────
+for (round in 1:n_rounds) {
+
+    # Forward pass — executes on federated workers.
+    scores = X %*% w;                           # (n × 1), federated
+    probs  = 1.0 / (1.0 + exp(-scores));        # (n × 1), federated
+
+    residuals = probs - y;                       # (n × 1), federated
+
+    # DP noise injection — only when private=1.
+    if (private == 1) {
+        # Per-example L2-norm clipping: each record's raw gradient
+        # contribution is X[i,:] * residual_i, with norm
+        # ||X[i,:]||_2 * |residual_i|. Scaling it down to clip_norm
+        # whenever it exceeds that bound is what makes
+        # sensitivity = clip_norm / n (set above) an actual, provable
+        # bound instead of an arbitrary constant.
+        row_norms         = sqrt(rowSums(X^2));               # (n × 1) ||X[i,:]||_2
+        contrib_norms     = abs(residuals) * row_norms;        # (n × 1) ||X[i,:]*residual_i||_2
+        clip_scale        = clip_norm / pmax(contrib_norms, clip_norm); # (n × 1), in (0,1]
+        clipped_residuals = residuals * clip_scale;            # (n × 1)
+
+        # Gradient aggregation — t(X) %*% clipped_residuals is a (d × 1)
+        # local sum that the coordinator collects in one federated
+        # aggregate instruction.
+        grad = t(X) %*% clipped_residuals / n;                # (d × 1), LOCAL after agg, clipped
+
+        noisy_grad = dp_gaussian(grad,
+                                 "identity",
+                                 sensitivity=sensitivity,
+                                 epsilon=round_epsilon,
+                                 delta=delta);
+        w = w - lr * noisy_grad;
+    } else {
+        # Gradient aggregation — t(X) %*% residuals is an (d × 1) local sum
+        # that the coordinator collects in one federated aggregate instruction.
+        grad = t(X) %*% residuals / n;                        # (d × 1), LOCAL after agg
+        w = w - lr * grad;
+    }
+}
+
+# ── Write model weights ───────────────────────────────────────────────────
+write(w, out_path, format="csv");
+print("FedAvg done. epsilon=" + epsilon + "  rounds=" + n_rounds);
+

benchmark/scripts/plot.py

@@ -0,0 +1,114 @@
+"""
+Read results.csv and produce two figures:
+
+1. accuracy_vs_epsilon.png
+   Line plot: x = ε, y = accuracy.
+   Horizontal dashed line = non-private baseline.
+   Points labelled with accuracy values.
+
+2. privacy_cost.png
+   Bar chart showing accuracy loss relative to baseline (utility cost of DP).
+"""
+import pathlib
+import csv
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import matplotlib.ticker as ticker
+
+RESULTS = pathlib.Path("benchmark/results")
+
+# ── Load ──────────────────────────────────────────────────────────────────
+rows = []
+with open(RESULTS / "results.csv") as f:
+    for r in csv.DictReader(f):
+        rows.append({
+            "label":    r["label"],
+            "epsilon":  float(r["epsilon"]) if r["epsilon"] != "inf" else None,
+            "private":  int(r["private"]),
+            "accuracy": float(r["accuracy"]),
+        })
+
+baseline = next(r for r in rows if r["private"] == 0)
+dp_rows  = sorted([r for r in rows if r["private"] == 1],
+                  key=lambda r: r["epsilon"])
+
+eps_vals = [r["epsilon"] for r in dp_rows]
+acc_vals = [r["accuracy"] for r in dp_rows]
+baseline_acc = baseline["accuracy"]
+
+# ── Figure 1: Accuracy vs ε ───────────────────────────────────────────────
+fig, ax = plt.subplots(figsize=(7, 4.5))
+
+ax.plot(eps_vals, acc_vals, marker="o", linewidth=2,
+        color="#028090", label="DP-FedAvg (Gaussian)")
+ax.axhline(baseline_acc, linestyle="--", color="#1C3A5E",
+           linewidth=1.5, label=f"Non-private baseline ({baseline_acc:.3f})")
+
+# Annotate each DP point.
+for eps, acc in zip(eps_vals, acc_vals):
+    ax.annotate(f"{acc:.3f}", xy=(eps, acc),
+                xytext=(0, 8), textcoords="offset points",
+                ha="center", fontsize=9, color="#028090")
+
+ax.set_xscale("log")
+ax.set_xticks(eps_vals)
+ax.get_xaxis().set_major_formatter(ticker.ScalarFormatter())
+ax.set_xlabel("Privacy budget ε  (smaller = stronger privacy)", fontsize=11)
+ax.set_ylabel("Test accuracy", fontsize=11)
+ax.set_title("Accuracy vs. Privacy Budget — DP-FedAvg on Adult (4 workers)",
+             fontsize=12)
+ax.legend(fontsize=9)
+ax.set_ylim(max(0, min(acc_vals) - 0.05), min(1.0, baseline_acc + 0.05))
+ax.grid(True, which="both", linestyle=":", alpha=0.5)
+
+plt.tight_layout()
+out1 = RESULTS / "accuracy_vs_epsilon.png"
+fig.savefig(out1, dpi=150)
+print(f"Saved {out1}")
+plt.close()
+
+# ── Figure 2: Utility cost (accuracy drop) ────────────────────────────────
+fig, ax = plt.subplots(figsize=(6, 4))
+
+drops = [baseline_acc - acc for acc in acc_vals]
+colors = ["#B91C1C" if d > 0.02 else "#028090" for d in drops]
+# Position bars at their true ε value on a log-scaled x-axis (rather than
+# evenly-spaced categorical slots) so the visual spacing between 0.5→1 and
+# 4→8 reflects the same 2x ratio. Bar widths scale with x so they stay a
+# constant fraction of their slot in log space instead of shrinking/growing.
+widths = [e * 0.4 for e in eps_vals]
+bars = ax.bar(eps_vals, drops, color=colors, width=widths,
+              edgecolor="white")
+ax.set_xscale("log")
+ax.set_xticks(eps_vals)
+ax.get_xaxis().set_major_formatter(ticker.ScalarFormatter())
+
+for bar, drop in zip(bars, drops):
+    ax.text(bar.get_x() + bar.get_width() / 2,
+            bar.get_height() + 0.001,
+            f"{drop:.3f}", ha="center", va="bottom", fontsize=9)
+
+ax.legend(["drop > 0.02", "drop <= 0.02"])
+
+ax.axhline(0, color="black", linewidth=0.8)
+ax.set_xlabel("Privacy budget ε", fontsize=11)
+ax.set_ylabel("Accuracy drop vs. baseline", fontsize=11)
+ax.set_title("Utility Cost of Differential Privacy — DP-FedAvg on Adult",
+             fontsize=11)
+ax.grid(True, axis="y", linestyle=":", alpha=0.5)
+
+plt.tight_layout()
+out2 = RESULTS / "privacy_cost.png"
+fig.savefig(out2, dpi=150)
+print(f"Saved {out2}")
+plt.close()
+
+# ── Console summary table ─────────────────────────────────────────────────
+print()
+print(f"{'ε':>8}  {'accuracy':>10}  {'drop':>8}")
+print("-" * 34)
+print(f"{'baseline':>8}  {baseline_acc:10.4f}  {'—':>8}")
+for eps, acc, drop in zip(eps_vals, acc_vals, drops):
+    print(f"{eps:>8.1f}  {acc:10.4f}  {drop:8.4f}")
+