DML syntax (post-aggregate form): + *
+ * result = dp_laplace(aggregate, sensitivity=1.0, epsilon=0.5) + * result = dp_gaussian(aggregate, sensitivity=1.0, epsilon=0.5, delta=1e-5) + *+ * + *
Sensitivity norm: {@code sensitivity} is not interchangeable + * between the two builtins. {@code dp_laplace} calibrates its noise scale + * to the L1 sensitivity of {@code aggregate} to a single-record + * change; {@code dp_gaussian} calibrates its σ to the L2 sensitivity. + * For a scalar aggregate (e.g. a single sum or mean) the two norms coincide, + * but for a vector-valued aggregate (e.g. column means of a multi-column + * matrix) they generally differ (L2 ≤ L1 ≤ √d·L2 for d entries) — the caller + * is responsible for supplying the norm matching the builtin invoked (see + * {@link #sensitivityOf}). + * + *
The instruction receives a materialised matrix (the aggregate result), + * injects calibrated noise element-wise, records the release with the + * session-scoped {@link DPBudgetAccountant}, and returns the noisy matrix. + * + *
Noise is generated in Java and added via a {@code MatrixBlock} binary + * operation so that the output allocation path is identical to every other + * CP instruction (no special memory-management required). + * + *
The {@link #sensitivityOf} method is deliberately separated from the
+ * noise-scale computation. In Phase 1 it returns the caller-supplied
+ * constant. In the future HOP-level rewrite pass (Phase 2) the body of this
+ * single method is replaced with a static analysis that reads the
+ * sensitivity bound computed by the compiler; every other line in this class
+ * stays unchanged.
+ */
+public class DPBuiltinCPInstruction extends ComputationCPInstruction {
+
+ // -----------------------------------------------------------------------
+ // Constants
+ // -----------------------------------------------------------------------
+
+ /** Opcode registered in Builtins and CPInstructionParser. */
+ public static final String OPCODE_LAPLACE = "dp_laplace";
+ public static final String OPCODE_GAUSSIAN = "dp_gaussian";
+
+ // -----------------------------------------------------------------------
+ // Fields
+ // -----------------------------------------------------------------------
+
+ /**
+ * Named parameters extracted from the serialised instruction string.
+ * Keys: "target", "sensitivity", "epsilon", "delta" (Gaussian only).
+ *
+ * Using the same LinkedHashMap Expected format (OPERAND_DELIM = '\u00b0'):
+ * Phase 1 (now): returns the caller-supplied literal from the
+ * DML script as-is, with no norm conversion or validation — the DML
+ * author must compute the sensitivity in the correct norm for the
+ * builtin they call. Sensitivity analysis is the caller's responsibility.
+ *
+ * Phase 2 (HOP-level rewrite pass): replace this body with a
+ * call that inspects the HOP node that produced {@code aggregate}, reads
+ * the {@code sensitivityBound} field computed during compilation (in the
+ * norm matching {@code instOpcode}), and returns it. No other line in
+ * this class changes.
+ *
+ * @param aggregate the already-computed aggregate block (ignored in
+ * Phase 1; used in Phase 2 to look up lineage)
+ * @return caller-supplied sensitivity constant, expected to already be
+ * in the L1 norm (Laplace) or L2 norm (Gaussian)
+ */
+ private double sensitivityOf(MatrixBlock aggregate) {
+ // Phase 1: unwrap the literal or variable value from the param map.
+ // In Phase 2, replace the body below with HOP-annotation lookup.
+ return parsePositiveDouble("sensitivity");
+ }
+
+ // -----------------------------------------------------------------------
+ // Noise generation
+ // -----------------------------------------------------------------------
+
+ /**
+ * Generates a noise {@link MatrixBlock} of the same shape as
+ * {@code aggregate}, filled with samples from the mechanism-appropriate
+ * distribution calibrated to ({@code sensitivity}, {@code epsilon},
+ * {@code delta}).
+ *
+ * Both mechanisms produce a dense block. Sparsity exploitation is
+ * left for future work; for the aggregate outputs targeted here (e.g.
+ * column means, row sums) the aggregate is already dense.
+ */
+ private MatrixBlock generateNoise(
+ MatrixBlock aggregate,
+ double sensitivity,
+ double epsilon,
+ double delta) {
+
+ int rows = aggregate.getNumRows();
+ int cols = aggregate.getNumColumns();
+ MatrixBlock noise = new MatrixBlock(rows, cols, false); // dense
+ noise.allocateDenseBlock();
+
+ if (instOpcode.equals(OPCODE_LAPLACE)) {
+ // Laplace mechanism
+ // For a given epsilon, noise is drawn from the Laplace distribution at
+ // scale b = sensitivity / epsilon
+ fillLaplaceNoise(noise, sensitivity / epsilon);
+ } else {
+ // Gaussian mechanism: calibrate sigma for (epsilon, delta)-DP.
+ // For a given epsilon, noise is drawn from the normal distribution at
+ // sigma^2 = 2 * sensitivity^2 * log(1.25/delta) / epsilon^2
+ double sigma = sensitivity
+ * Math.sqrt(2.0 * Math.log(1.25 / delta))
+ / epsilon;
+ fillGaussianNoise(noise, sigma);
+ }
+
+ noise.recomputeNonZeros();
+ return noise;
+ }
+
+ /**
+ * Fills {@code block} with i.i.d. Laplace(0, scale) samples using the
+ * inverse-CDF method.
+ *
+ * For u ~ Uniform(0, 1): X = -scale * sign(u - 0.5) * ln(1 - 2|u - 0.5|)
+ */
+ private static void fillLaplaceNoise(MatrixBlock block, double scale) {
+ ThreadLocalRandom rng = ThreadLocalRandom.current();
+ int rows = block.getNumRows();
+ int cols = block.getNumColumns();
+ for (int r = 0; r < rows; r++) {
+ for (int c = 0; c < cols; c++) {
+ double u = rng.nextDouble(); // u in (0, 1)
+ double v = u - 0.5;
+ // Guard against the degenerate u == 0.5 case (ln(0) = -inf).
+ if (v == 0.0) v = 1e-15;
+ double sample = -scale * Math.signum(v) * Math.log(1.0 - 2.0 * Math.abs(v));
+ block.set(r, c, sample);
+ }
+ }
+ }
+
+ /**
+ * Fills {@code block} with i.i.d. N(0, sigma²) samples.
+ *
+ * Uses {@link ThreadLocalRandom#nextGaussian()} which is thread-safe
+ * and does not require external libraries.
+ */
+ private static void fillGaussianNoise(MatrixBlock block, double sigma) {
+ ThreadLocalRandom rng = ThreadLocalRandom.current();
+ int rows = block.getNumRows();
+ int cols = block.getNumColumns();
+ for (int r = 0; r < rows; r++) {
+ for (int c = 0; c < cols; c++) {
+ block.set(r, c, sigma * rng.nextGaussian());
+ }
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Helpers
+ // -----------------------------------------------------------------------
+
+ /**
+ * Parses a parameter value as a positive {@code double}.
+ *
+ * @throws DMLRuntimeException if the key is absent, unparseable, or
+ * non-positive
+ */
+ private double parsePositiveDouble(String key) {
+ String raw = _params.get(key);
+ if (raw == null)
+ throw new DMLRuntimeException(
+ instOpcode + ": parameter '" + key + "' is missing");
+ double v;
+ try {
+ v = Double.parseDouble(raw);
+ } catch (NumberFormatException e) {
+ throw new DMLRuntimeException(
+ instOpcode + ": parameter '" + key
+ + "' is not a valid number: " + raw);
+ }
+ if (!(v > 0.0))
+ throw new DMLRuntimeException(
+ instOpcode + ": parameter '" + key
+ + "' must be strictly positive, got " + v);
+ return v;
+ }
+}
diff --git a/src/main/java/org/apache/sysds/runtime/privacy/dp/DPBudgetAccountant.java b/src/main/java/org/apache/sysds/runtime/privacy/dp/DPBudgetAccountant.java
new file mode 100644
index 00000000000..9dd011fcb1e
--- /dev/null
+++ b/src/main/java/org/apache/sysds/runtime/privacy/dp/DPBudgetAccountant.java
@@ -0,0 +1,283 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied. See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+package org.apache.sysds.runtime.privacy.dp;
+
+import org.apache.sysds.runtime.DMLRuntimeException;
+import org.apache.sysds.runtime.instructions.cp.DPBuiltinCPInstruction;
+
+/**
+ * Session-scoped differential privacy budget accountant.
+ *
+ * When both mechanisms are used in the same script the total cost is:
+ * Laplace gives pure ε-DP (no δ). Basic composition is exact and
+ * tighter than the RDP conversion path for Laplace (which would introduce
+ * an unnecessary δ and produce a looser bound). Each Laplace release adds
+ * its ε here; the total is added directly in {@link #totalEpsilonSpent()}.
+ */
+ private double _pureEpsilonSum = 0.0;
+
+ /** Total privacy budget (ε) for the script execution. */
+ private final double _epsilonBudget;
+
+ /** δ used for the Gaussian RDP-to-(ε,δ) conversion. */
+ private final double _delta;
+
+ /** Number of releases recorded so far (for error messages). */
+ private int _releaseCount = 0;
+
+ // -----------------------------------------------------------------------
+ // Constructors
+ // -----------------------------------------------------------------------
+
+ /**
+ * Creates an accountant with the given global budget.
+ *
+ * Typical usage: the DML script sets the budget once at the top
+ * (future work: a {@code dp_set_budget(epsilon, delta)} built-in),
+ * or the accountant is created with defaults and the budget is checked
+ * after each release.
+ *
+ * @param epsilonBudget total ε budget for the script execution (must be > 0)
+ * @param delta δ used for the Gaussian RDP-to-(ε,δ) conversion (must be in (0,1))
+ */
+ public DPBudgetAccountant(double epsilonBudget, double delta) {
+ if (!(epsilonBudget > 0))
+ throw new DMLRuntimeException(
+ "DPBudgetAccountant: epsilonBudget must be > 0, got " + epsilonBudget);
+ if (!(delta > 0 && delta < 1))
+ throw new DMLRuntimeException(
+ "DPBudgetAccountant: delta must be in (0,1), got " + delta);
+ _epsilonBudget = epsilonBudget;
+ _delta = delta;
+ }
+
+ /**
+ * Convenience constructor using a liberal default δ = 1e-5.
+ * Suitable when the calling script does not specify δ explicitly.
+ */
+ public DPBudgetAccountant(double epsilonBudget) {
+ this(epsilonBudget, 1e-5);
+ }
+
+ // -----------------------------------------------------------------------
+ // Core API
+ // -----------------------------------------------------------------------
+
+ /**
+ * Records one DP release and checks the budget.
+ *
+ * This method must be called before the result is written to
+ * the variable table. If the budget is exhausted it throws and the
+ * caller's result is discarded, preventing an unaccounted release.
+ *
+ * Mechanism selection (see class-level Javadoc for details):
+ * Total = Laplace pure-ε sum + Gaussian RDP-converted ε (clamped to
+ * zero when no Gaussian releases have been recorded).
+ */
+ public double totalEpsilonSpent() {
+ // Take min_α(ε_α) as the current total privacy cost
+ double gaussianEps = Double.MAX_VALUE;
+ for (int i = 0; i < ORDERS.length; i++) {
+ double alpha = ORDERS[i];
+ double eps = _rdpSum[i]
+ + Math.log(1.0 - 1.0 / alpha)
+ - Math.log(_delta * (alpha - 1.0)) / alpha;
+ if (eps < gaussianEps)
+ gaussianEps = eps;
+ }
+ // Clamp: with no Gaussian releases the RDP sum is 0 and the log-delta
+ // term alone drives gaussianEps to a small positive value; clamp to 0
+ // so Laplace-only scripts are not penalised by δ they never requested.
+ if (gaussianEps < 0) gaussianEps = 0.0;
+ return _pureEpsilonSum + gaussianEps;
+ }
+
+ /** Returns the remaining ε budget (negative if the budget is exceeded). */
+ public double remainingBudget() {
+ return _epsilonBudget - totalEpsilonSpent();
+ }
+
+ /** Returns the number of DP releases recorded so far. */
+ public int releaseCount() {
+ return _releaseCount;
+ }
+
+ // -----------------------------------------------------------------------
+ // Private helpers
+ // -----------------------------------------------------------------------
+
+ /**
+ * Rényi divergence of order α for the Gaussian mechanism (Mironov 2017,
+ * Proposition 3, example 2):
+ * The tests are grouped into three levels:
+ * The DML integration tests require a built SystemDS jar and are separated
+ * into a companion class {@link org.apache.sysds.test.functions.privacy.dp.DPBuiltinDMLTest}.
+ */
+public class DPBuiltinCPInstructionTest {
+
+ private static final double EPS = 1e-9;
+
+ // =======================================================================
+ // 1. DPBudgetAccountant unit tests
+ // =======================================================================
+
+ @Test
+ public void testAccountantInitialisesAtZeroCost() {
+ DPBudgetAccountant acc = new DPBudgetAccountant(1.0, 1e-5);
+ // No releases yet: total cost should be a large negative number
+ // (conversion formula gives -∞ when rdpSum = 0 for all orders),
+ // so remainingBudget() should exceed the budget.
+ assertTrue("No releases should leave budget intact",
+ acc.remainingBudget() > 0);
+ assertEquals(0, acc.releaseCount());
+ }
+
+ @Test
+ public void testSingleLaplaceReleaseDoesNotExceedBudget() {
+ // epsilon=0.5, budget=1.0: one release should consume < budget.
+ DPBudgetAccountant acc = new DPBudgetAccountant(1.0, 1e-5);
+ acc.compose(0.5, 0.0, 1.0); // Laplace, sensitivity=1
+ assertEquals(1, acc.releaseCount());
+ assertTrue("Single release within budget",
+ acc.totalEpsilonSpent() <= 1.0);
+ }
+
+ @Test
+ public void testSingleGaussianReleaseDoesNotExceedBudget() {
+ DPBudgetAccountant acc = new DPBudgetAccountant(1.0, 1e-5);
+ acc.compose(0.5, 1e-5, 1.0); // Gaussian
+ assertEquals(1, acc.releaseCount());
+ assertTrue("Single Gaussian release within budget",
+ acc.totalEpsilonSpent() <= 1.0);
+ }
+
+ @Test(expected = DMLRuntimeException.class)
+ public void testBudgetExhaustionThrows() {
+ // Budget = 0.1, but we try to make 10 releases at epsilon=0.5 each.
+ // After enough releases the budget must be exceeded.
+ DPBudgetAccountant acc = new DPBudgetAccountant(0.1, 1e-5);
+ for (int i = 0; i < 10; i++) {
+ acc.compose(0.5, 0.0, 1.0); // will throw before the 10th
+ }
+ }
+
+ @Test
+ public void testCompositionIsMonotonicallyIncreasing() {
+ DPBudgetAccountant acc = new DPBudgetAccountant(100.0, 1e-5); // large budget
+ double prev = acc.totalEpsilonSpent();
+ for (int i = 0; i < 5; i++) {
+ acc.compose(0.3, 1e-5, 1.0);
+ double current = acc.totalEpsilonSpent();
+ assertTrue("Epsilon spent must increase with each release",
+ current > prev);
+ prev = current;
+ }
+ }
+
+ @Test
+ public void testGaussianTighterThanLaplaceForSameEpsilon() {
+ // For the same nominal (ε, δ), Gaussian uses RDP composition which
+ // is tighter than Laplace with basic composition. After 5 releases:
+ // Laplace (basic, worst-case): 5ε
+ // Gaussian (RDP) : something < 5ε
+ double eps = 0.5;
+ double delta = 1e-5;
+
+ DPBudgetAccountant gaussian = new DPBudgetAccountant(100.0, delta);
+ DPBudgetAccountant laplace = new DPBudgetAccountant(100.0, delta);
+
+ for (int i = 0; i < 5; i++) {
+ gaussian.compose(eps, delta, 1.0);
+ laplace.compose(eps, 0.0, 1.0);
+ }
+
+ // After 5 releases, Gaussian RDP bound should be tighter.
+ // (Both may be < 5*eps; the point is Gaussian <= Laplace.)
+ assertTrue("Gaussian RDP bound should be <= Laplace bound after 5 releases",
+ gaussian.totalEpsilonSpent() <= laplace.totalEpsilonSpent() + 1e-6);
+ }
+
+ @Test
+ public void testRemainingBudgetDecreasesMonotonically() {
+ DPBudgetAccountant acc = new DPBudgetAccountant(2.0, 1e-5);
+ double prev = acc.remainingBudget();
+ for (int i = 0; i < 3; i++) {
+ acc.compose(0.2, 1e-5, 1.0);
+ double current = acc.remainingBudget();
+ assertTrue("Remaining budget must decrease", current < prev);
+ prev = current;
+ }
+ }
+
+ @Test
+ public void testHigherEpsilonCostMoreForLaplace() {
+ // For Laplace, the accountant uses basic (pure ε-DP) composition: cost = epsilon.
+ // Sensitivity determines noise scale but NOT the budget consumed — that is set
+ // entirely by the caller's epsilon parameter.
+ // A release at epsilon=1.0 costs more budget than one at epsilon=0.5.
+ DPBudgetAccountant acc1 = new DPBudgetAccountant(100.0, 1e-5);
+ DPBudgetAccountant acc2 = new DPBudgetAccountant(100.0, 1e-5);
+ acc1.compose(0.5, 0.0, 1.0); // epsilon=0.5, Laplace
+ acc2.compose(1.0, 0.0, 1.0); // epsilon=1.0, same sensitivity
+
+ assertTrue("Higher epsilon costs more budget (Laplace basic composition)",
+ acc1.totalEpsilonSpent() < acc2.totalEpsilonSpent());
+ }
+
+ // --- Constructor error paths ------------------------------------
+
+ @Test(expected = DMLRuntimeException.class)
+ public void testConstructorRejectsZeroEpsilonBudget() {
+ new DPBudgetAccountant(0.0, 1e-5);
+ }
+
+ @Test(expected = DMLRuntimeException.class)
+ public void testConstructorRejectsNegativeEpsilonBudget() {
+ new DPBudgetAccountant(-0.5, 1e-5);
+ }
+
+ @Test(expected = DMLRuntimeException.class)
+ public void testConstructorRejectsDeltaZero() {
+ new DPBudgetAccountant(1.0, 0.0);
+ }
+
+ @Test(expected = DMLRuntimeException.class)
+ public void testConstructorRejectsDeltaOne() {
+ new DPBudgetAccountant(1.0, 1.0);
+ }
+
+ // --- Single-argument convenience constructor -------------------
+
+ @Test
+ public void testConvenienceConstructorDefaultsDeltaTo1e5() {
+ // The one-arg form delegates to (epsilonBudget, 1e-5). A Gaussian
+ // release whose per-release delta matches that default must produce
+ // identical totalEpsilonSpent() from both construction paths.
+ DPBudgetAccountant oneArg = new DPBudgetAccountant(10.0);
+ DPBudgetAccountant twoArg = new DPBudgetAccountant(10.0, 1e-5);
+ oneArg.compose(0.5, 1e-5, 1.0);
+ twoArg.compose(0.5, 1e-5, 1.0);
+ assertEquals("Convenience constructor must default to delta=1e-5",
+ twoArg.totalEpsilonSpent(), oneArg.totalEpsilonSpent(), EPS);
+ }
+
+ @Test(expected = DMLRuntimeException.class)
+ public void testGaussianBudgetExhaustionThrows() {
+ // Budget = 0.1. Each Gaussian release costs more than 0.005, so 20
+ // releases must exceed the budget well before the loop ends.
+ DPBudgetAccountant acc = new DPBudgetAccountant(0.1, 1e-5);
+ for (int i = 0; i < 20; i++) {
+ acc.compose(0.3, 1e-5, 1.0);
+ }
+ }
+
+ @Test
+ public void testMixedCompositionExceedsEitherAlone() {
+ // Compose one Laplace and one Gaussian release. The total cost must
+ // exceed what either mechanism contributes alone, exercising the
+ // _pureEpsilonSum + gaussianEps addition path in totalEpsilonSpent().
+ DPBudgetAccountant mixed = new DPBudgetAccountant(100.0, 1e-5);
+ DPBudgetAccountant lapOnly = new DPBudgetAccountant(100.0, 1e-5);
+ DPBudgetAccountant gauOnly = new DPBudgetAccountant(100.0, 1e-5);
+
+ mixed.compose(0.5, 0.0, 1.0); // Laplace
+ mixed.compose(0.5, 1e-5, 1.0); // Gaussian
+
+ lapOnly.compose(0.5, 0.0, 1.0);
+ gauOnly.compose(0.5, 1e-5, 1.0);
+
+ assertTrue("Mixed cost must exceed Laplace-only cost",
+ mixed.totalEpsilonSpent() > lapOnly.totalEpsilonSpent());
+ assertTrue("Mixed cost must exceed Gaussian-only cost",
+ mixed.totalEpsilonSpent() > gauOnly.totalEpsilonSpent());
+ }
+
+ // --- Release count across multiple mixed releases --------------
+
+ @Test
+ public void testReleaseCountTracksAllReleases() {
+ DPBudgetAccountant acc = new DPBudgetAccountant(100.0, 1e-5);
+ assertEquals(0, acc.releaseCount());
+ acc.compose(0.1, 0.0, 1.0); // Laplace
+ assertEquals(1, acc.releaseCount());
+ acc.compose(0.1, 1e-5, 1.0); // Gaussian
+ assertEquals(2, acc.releaseCount());
+ acc.compose(0.1, 0.0, 1.0); // Laplace
+ acc.compose(0.1, 0.0, 1.0); // Laplace
+ acc.compose(0.1, 1e-5, 1.0); // Gaussian
+ assertEquals(5, acc.releaseCount());
+ }
+
+ // --- Edge-case inputs for rdpGaussian / gaussianSigma ----------
+
+ @Test
+ public void testGaussianSensitivityCancelsInRDP() {
+ // For the Gaussian mechanism: σ = Δf·sqrt(2·ln(1.25/δ))/ε, so
+ // D_α = α·Δf²/(2σ²) = α·ε²/(4·ln(1.25/δ)).
+ // Sensitivity cancels. Two accountants with the same (ε,δ) but
+ // different sensitivity must report identical totalEpsilonSpent().
+ DPBudgetAccountant acc1 = new DPBudgetAccountant(100.0, 1e-5);
+ DPBudgetAccountant acc2 = new DPBudgetAccountant(100.0, 1e-5);
+ acc1.compose(0.5, 1e-5, 1.0); // sensitivity = 1
+ acc2.compose(0.5, 1e-5, 100.0); // sensitivity = 100, same (ε,δ)
+ assertEquals("Gaussian RDP cost must be independent of sensitivity when (ε,δ) are fixed",
+ acc1.totalEpsilonSpent(), acc2.totalEpsilonSpent(), EPS);
+ }
+
+ @Test
+ public void testGaussianLargerEpsilonCostsMoreBudget() {
+ // D_α ∝ ε², so a release declared at a higher ε (less noise, more
+ // privacy loss) must cost more budget than one at a lower ε.
+ DPBudgetAccountant lowEps = new DPBudgetAccountant(100.0, 1e-5);
+ DPBudgetAccountant highEps = new DPBudgetAccountant(100.0, 1e-5);
+ lowEps.compose(0.1, 1e-5, 1.0);
+ highEps.compose(0.5, 1e-5, 1.0);
+ assertTrue("Larger epsilon per Gaussian release must cost more budget",
+ highEps.totalEpsilonSpent() > lowEps.totalEpsilonSpent());
+ }
+
+ // =======================================================================
+ // 2. Noise distribution tests (statistical sanity checks)
+ // =======================================================================
+ // These tests generate many samples and verify that the empirical mean
+ // is near zero and the empirical variance matches the theoretical value
+ // within a reasonable tolerance.
+ //
+
+ @Test
+ public void testLaplaceNoiseMeanNearZero() {
+ // For 10000 samples the empirical mean should be within 3σ/√n of 0.
+ int n = 10_000;
+ double scale = 2.0;
+ double[] samples = sampleLaplace(n, scale);
+ double mean = mean(samples);
+ double theoreticalStdErr = scale * Math.sqrt(2.0) / Math.sqrt(n);
+ assertTrue("Laplace mean should be near 0",
+ Math.abs(mean) < 5 * theoreticalStdErr);
+ }
+
+ @Test
+ public void testLaplaceNoiseVarianceCorrect() {
+ // Var[Laplace(0, b)] = 2b². Allow 10% relative error for n=10000.
+ int n = 10_000;
+ double scale = 1.5;
+ double[] samples = sampleLaplace(n, scale);
+ double variance = variance(samples);
+ double expected = 2.0 * scale * scale;
+ assertEquals("Laplace variance", expected, variance, 0.1 * expected);
+ }
+
+ @Test
+ public void testGaussianNoiseMeanNearZero() {
+ int n = 10_000;
+ double sigma = 3.0;
+ double[] samples = sampleGaussian(n, sigma);
+ double mean = mean(samples);
+ double theoreticalStdErr = sigma / Math.sqrt(n);
+ assertTrue("Gaussian mean should be near 0",
+ Math.abs(mean) < 5 * theoreticalStdErr);
+ }
+
+ @Test
+ public void testGaussianNoiseVarianceCorrect() {
+ int n = 10_000;
+ double sigma = 2.0;
+ double[] samples = sampleGaussian(n, sigma);
+ double variance = variance(samples);
+ double expected = sigma * sigma;
+ assertEquals("Gaussian variance", expected, variance, 0.1 * expected);
+ }
+
+ // -----------------------------------------------------------------------
+ // Helpers for noise distribution tests
+ // -----------------------------------------------------------------------
+
+ /** Sample n Laplace(0, scale) values using the inverse-CDF method. */
+ private static double[] sampleLaplace(int n, double scale) {
+ java.util.concurrent.ThreadLocalRandom rng =
+ java.util.concurrent.ThreadLocalRandom.current();
+ double[] out = new double[n];
+ for (int i = 0; i < n; i++) {
+ double u = rng.nextDouble();
+ double v = u - 0.5;
+ if (v == 0.0) v = 1e-15;
+ out[i] = -scale * Math.signum(v) * Math.log(1.0 - 2.0 * Math.abs(v));
+ }
+ return out;
+ }
+
+ /** Sample n N(0, sigma²) values. */
+ private static double[] sampleGaussian(int n, double sigma) {
+ java.util.concurrent.ThreadLocalRandom rng =
+ java.util.concurrent.ThreadLocalRandom.current();
+ double[] out = new double[n];
+ for (int i = 0; i < n; i++) {
+ out[i] = sigma * rng.nextGaussian();
+ }
+ return out;
+ }
+
+ private static double mean(double[] xs) {
+ double s = 0;
+ for (double x : xs) s += x;
+ return s / xs.length;
+ }
+
+ private static double variance(double[] xs) {
+ double m = mean(xs);
+ double s = 0;
+ for (double x : xs) s += (x - m) * (x - m);
+ return s / (xs.length - 1);
+ }
+}
diff --git a/src/test/java/org/apache/sysds/test/functions/privacy/dp/DPBuiltinDMLTest.java b/src/test/java/org/apache/sysds/test/functions/privacy/dp/DPBuiltinDMLTest.java
new file mode 100644
index 00000000000..048c9501e21
--- /dev/null
+++ b/src/test/java/org/apache/sysds/test/functions/privacy/dp/DPBuiltinDMLTest.java
@@ -0,0 +1,130 @@
+// ==========================================================================
+// DML integration test
+// ==========================================================================
+//
+// Full integration tests extend AutomatedTestBase and drive the DML runner.
+// Each test:
+// (a) Writes a DML script to a temp file.
+// (b) Provides input matrices via TestUtils.
+// (c) Calls runTest() and reads the output MatrixBlock.
+// (d) Verifies that the noisy result differs from the clean result by a
+// statistically plausible amount (not zero, not astronomically large).
+//
+
+
+package org.apache.sysds.test.functions.privacy.dp;
+
+import static org.junit.Assert.assertTrue;
+
+import java.io.File;
+import java.io.IOException;
+import java.nio.file.Files;
+import java.util.HashMap;
+
+import org.apache.sysds.runtime.matrix.data.MatrixValue.CellIndex;
+import org.apache.sysds.test.AutomatedTestBase;
+import org.apache.sysds.test.TestConfiguration;
+import org.apache.sysds.test.TestUtils;
+import org.junit.Test;
+
+public class DPBuiltinDMLTest extends AutomatedTestBase {
+
+ private static final String TEST_DIR = "functions/privacy/dp/";
+ private static final String TEST_CLASS = TEST_DIR + DPBuiltinDMLTest.class.getSimpleName() + "/";
+ private static final int ROWS = 100;
+ private static final int COLS = 10;
+ private static final String DML_LAPLACE =
+ "X = read($1);\n"
+ + "result = dp_laplace(colMeans(X), sensitivity=1.0, epsilon=$2);\n"
+ + "write(result, $3, format=\"text\");\n";
+ private static final String DML_GAUSSIAN =
+ "X = read($1);\n"
+ + "result = dp_gaussian(colMeans(X), sensitivity=1.0, epsilon=$2, delta=1e-5);\n"
+ + "write(result, $3, format=\"text\");\n";
+
+ @Override
+ public void setUp() {
+ addTestConfiguration("DPLaplace", new TestConfiguration(TEST_CLASS, "DPLaplace"));
+ addTestConfiguration("DPGaussian", new TestConfiguration(TEST_CLASS, "DPGaussian"));
+ }
+
+ @Test
+ public void testLaplaceOutputDiffersFromCleanMean() {
+ runDPTest("DPLaplace", DML_LAPLACE, "0.5");
+ }
+
+ @Test
+ public void testGaussianOutputDiffersFromCleanMean() {
+ runDPTest("DPGaussian", DML_GAUSSIAN, "0.5");
+ }
+
+ @Test
+ public void testHighEpsilonIsCloserToTruth() {
+ double[][] data = TestUtils.generateTestMatrix(ROWS, COLS, 0, 1, 1.0, 42);
+ // Higher ε → less noise → result closer to the true mean.
+ // NOTE: the DPBudgetAccountant caps total spend at the default budget
+ // (ε = 1.0) regardless of the per-release ε requested, so ε values
+ // here must stay well under that cap or the release is rejected.
+ double noisyLow = runAndGetMaxAbsColMeansDiffFromClean(data, "DPGaussian", DML_GAUSSIAN, "0.1");
+ double noisyHigh = runAndGetMaxAbsColMeansDiffFromClean(data, "DPGaussian", DML_GAUSSIAN, "0.5");
+ assertTrue("ε=0.5 should give less noise than ε=0.1", noisyHigh < noisyLow);
+ }
+
+ private void runDPTest(String testName, String dml, String epsilonStr) {
+ double[][] data = TestUtils.generateTestMatrix(ROWS, COLS, 0, 1, 1.0, 42);
+ HashMap
+ * dp_gaussian°target=mVar1·MATRIX·FP64°sensitivity=1.0·SCALAR·FP64·true
+ * °epsilon=0.5·SCALAR·FP64·true°delta=1e-5·SCALAR·FP64·true
+ * °_mVar2·MATRIX·FP64
+ *
+ *
+ * The first token is always the opcode; the last token is always the
+ * output operand; the tokens in between are key=value pairs. This matches
+ * the convention used by ParameterizedBuiltinCPInstruction exactly.
+ */
+ public static DPBuiltinCPInstruction parseInstruction(String str) {
+ String[] parts = InstructionUtils.getInstructionPartsWithValueType(str);
+ InstructionUtils.checkNumFields(parts, 4, 5); // laplace=4, gaussian=5
+ String opcode = parts[0];
+
+ // Output operand is always the last token.
+ CPOperand output = new CPOperand(parts[parts.length - 1]);
+
+ // The "target" parameter holds the variable name of the input matrix.
+ // ParameterizedBuiltinCPInstruction.constructParameterMap strips the
+ // type suffixes and returns bare key=value pairs.
+ LinkedHashMap
+ *
+ */
+ @Override
+ public void processInstruction(ExecutionContext ec) {
+
+ // ── 1. Read aggregate input ─────────────────────────────────────────
+ // getMatrixInput pins the block in memory and increments the
+ // reference count; we must call releaseMatrixInput afterwards.
+ MatrixBlock inBlock = ec.getMatrixInput(_params.get("target"));
+
+ // ── 2. Parse DP parameters ──────────────────────────────────────────
+ double epsilon = parsePositiveDouble("epsilon");
+ double delta = instOpcode.equals(OPCODE_GAUSSIAN)
+ ? parsePositiveDouble("delta") : 0.0;
+
+ // ── 3. Determine sensitivity (Phase-1: caller-supplied constant) ────
+ double sensitivity = sensitivityOf(inBlock);
+
+ // ── 4. Generate and add noise ────────────────────────────────────────
+ MatrixBlock noiseBlock = generateNoise(inBlock, sensitivity, epsilon, delta);
+
+ // Element-wise addition via the standard binary-operator path.
+ // binaryOperations allocates the output block internally.
+ BinaryOperator plusOp = new BinaryOperator(Plus.getPlusFnObject());
+ MatrixBlock outBlock = new MatrixBlock();
+ inBlock.binaryOperations(plusOp, noiseBlock, outBlock);
+
+ // ── 5. Record release and enforce budget ────────────────────────────
+ // getDPBudgetAccountant() returns a lazy-initialised DPBudgetAccountant that is
+ // owned by this ExecutionContext (added in a companion EC patch).
+ DPBudgetAccountant accountant = ec.getDPBudgetAccountant();
+ accountant.compose(epsilon, delta, sensitivity); // throws on exhaustion
+
+ // ── 6. Write output and release input pin ───────────────────────────
+ ec.releaseMatrixInput(_params.get("target"));
+ ec.setMatrixOutput(output.getName(), outBlock);
+ }
+
+ // -----------------------------------------------------------------------
+ // Sensitivity seam (Phase-1 stub; Phase-2 replaces this body only)
+ // -----------------------------------------------------------------------
+
+ /**
+ * Returns the sensitivity of {@code aggregate} to a single-record change,
+ * in the norm required by the mechanism actually invoked: L1 for
+ * {@code dp_laplace}, L2 for {@code dp_gaussian} (see the class
+ * Javadoc). The two only coincide when {@code aggregate} is scalar.
+ *
+ * Purpose
+ * Tracks composition of DP releases across the lifetime of a DML script
+ * execution. Each call to {@link #compose} records one release and checks
+ * whether the cumulative privacy cost has exceeded the user-specified budget.
+ *
+ * Composition strategy
+ * The mechanism type (Laplace vs Gaussian) is inferred from the {@code delta}
+ * argument passed to {@link #compose}:
+ *
+ *
+ *
+ *
+ *
+ * ε_total = ε_Laplace_sum + ε_Gaussian_RDP
+ *
+ * This follows from basic composition of a pure-DP mechanism with an
+ * approximate-DP mechanism, which is additive in ε.
+ *
+ * Rényi orders tracked (Gaussian path)
+ * α ∈ {2, 4, 8, 16, 32, 64, 128, 256, 512, 1024}. At query time the minimum
+ * converted ε across all orders is taken as the tightest available bound.
+ *
+ * Gaussian RDP divergence
+ * For the Gaussian mechanism with noise scale σ and L2 sensitivity Δf:
+ *
+ * D_α = α · Δf² / (2σ²)
+ *
+ * σ is back-derived from the caller's (ε, δ) via the standard calibration
+ * formula (see {@link #gaussianSigma}). Note that sensitivity cancels in the
+ * final expression, so the RDP cost depends only on the (ε, δ) parameters.
+ *
+ * RDP → (ε, δ) conversion (Mironov 2017, Proposition 3)
+ *
+ * ε(α) = R[α] + log(1 − 1/α) − log(δ·(α−1)) / α
+ *
+ *
+ * Lifecycle
+ * One instance is created per {@code ExecutionContext} (lazy init). It is
+ * garbage-collected with the context when the script finishes; no state
+ * leaks between script executions or between concurrent scripts.
+ *
+ * Thread safety
+ * Not thread-safe. A single DML script executes instructions sequentially
+ * on one thread, so no synchronisation is needed.
+ *
+ * @see DPBuiltinCPInstruction
+ */
+public class DPBudgetAccountant {
+
+ // -----------------------------------------------------------------------
+ // Rényi orders used for Gaussian composition
+ // -----------------------------------------------------------------------
+
+ /**
+ * Discrete set of Rényi orders α. All must be > 1.
+ * Finer grids give tighter bounds; this set covers the range relevant
+ * for typical ML workloads.
+ */
+ private static final double[] ORDERS = {
+ 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024
+ };
+
+ // -----------------------------------------------------------------------
+ // State
+ // -----------------------------------------------------------------------
+
+ /** Accumulated Rényi divergence at each order (Gaussian releases only). */
+ private final double[] _rdpSum = new double[ORDERS.length];
+
+ /**
+ * Running sum of pure ε from Laplace releases.
+ *
+ *
+ *
+ *
+ * @param epsilon per-release ε parameter (must be > 0)
+ * @param delta per-release δ parameter (0 for Laplace, >0 for Gaussian)
+ * @param sensitivity sensitivity Δf of the released quantity (must be > 0).
+ * The norm depends on the mechanism selected by
+ * {@code delta}: callers must supply the
+ * L1 sensitivity ‖f(D) − f(D′)‖₁ when
+ * {@code delta == 0} (Laplace), and the L2
+ * sensitivity ‖f(D) − f(D′)‖₂ when {@code delta > 0}
+ * (Gaussian). The two coincide for scalar-valued
+ * releases but diverge for vector-valued ones, so
+ * passing the wrong norm silently under- or
+ * over-calibrates the noise.
+ * @throws DMLRuntimeException if the cumulative ε after this release
+ * would exceed the budget
+ */
+ public void compose(double epsilon, double delta, double sensitivity) {
+ _releaseCount++;
+
+ if (delta == 0.0) {
+ // Laplace: pure ε-DP, basic composition — cost is exactly epsilon.
+ _pureEpsilonSum += epsilon;
+ } else {
+ // Gaussian: accumulate Rényi divergence at each order, then convert.
+ for (int i = 0; i < ORDERS.length; i++) {
+ double sigma = gaussianSigma(sensitivity, epsilon, delta);
+ _rdpSum[i] += rdpGaussian(ORDERS[i], sensitivity, sigma);
+ }
+ }
+
+ double spentEpsilon = totalEpsilonSpent();
+ if (spentEpsilon > _epsilonBudget) {
+ throw new DMLRuntimeException(String.format(
+ "Privacy budget exhausted after %d release(s): "
+ + "spent ε ≈ %.6f exceeds budget ε = %.6f (δ = %.2e). "
+ + "Reduce the number of releases or widen the budget.",
+ _releaseCount, spentEpsilon, _epsilonBudget, _delta));
+ }
+ }
+
+ // -----------------------------------------------------------------------
+ // Inspection
+ // -----------------------------------------------------------------------
+
+ /**
+ * Returns the current total privacy cost as an ε value.
+ *
+ *
+ * D_α = α · Δf² / (2σ²)
+ *
+ */
+ private static double rdpGaussian(double alpha, double sensitivity, double sigma) {
+ return alpha * (sensitivity * sensitivity) / (2.0 * sigma * sigma);
+ }
+
+ /**
+ * Gaussian noise scale σ calibrated to (ε, δ)-DP:
+ *
+ * σ = Δf · sqrt(2 · log(1.25 / δ)) / ε
+ *
+ * Must match the formula used in {@link DPBuiltinCPInstruction} so that
+ * the RDP cost recorded here is consistent with the noise actually injected.
+ */
+ private static double gaussianSigma(double sensitivity, double epsilon, double delta) {
+ return sensitivity * Math.sqrt(2.0 * Math.log(1.25 / delta)) / epsilon;
+ }
+}
diff --git a/src/test/java/org/apache/sysds/test/component/cp/DPBuiltinCPInstructionTest.java b/src/test/java/org/apache/sysds/test/component/cp/DPBuiltinCPInstructionTest.java
new file mode 100755
index 00000000000..750451b2991
--- /dev/null
+++ b/src/test/java/org/apache/sysds/test/component/cp/DPBuiltinCPInstructionTest.java
@@ -0,0 +1,362 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied. See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+
+package org.apache.sysds.test.component.cp;
+
+import org.apache.sysds.runtime.DMLRuntimeException;
+import org.apache.sysds.runtime.privacy.dp.DPBudgetAccountant;
+import org.junit.Test;
+
+import static org.junit.Assert.*;
+
+/**
+ * Tests for {@code DPBuiltinCPInstruction} and {@code DPBudgetAccountant}.
+ *
+ *
+ *
+ *
+ *