[API Gap] Marginal scalar accounting yields optimistic estimates for jointly observed correlated Gaussian releases

## [API Gap] Marginal scalar accounting yields optimistic estimates for jointly observed correlated Gaussian releases

**Type:** Documentation / API modeling gap
**Severity:** Medium — marginal accounting silently optimistic when used on correlated multi-output releases

---

### Context

This issue is not a mathematical error in the implementation. The library accurately
computes privacy guarantees for scalar Gaussian mechanisms under the defined model.

The gap is at the API affordance level: the current API exposes only scalar Gaussian
accounting primitives, which users may naturally apply marginally to a correlated
multi-output deployment, obtaining optimistic ε-hockey-stick divergence estimates
relative to the true joint release observed by an adversary.

---

### Setup

Consider a mechanism releasing **two scalar outputs simultaneously**:

```
M(D) = (f_A(D) + Z_A, f_B(D) + Z_B)
```

where:
- `Z_A, Z_B ~ N(0, σ²)` **independently** (no correlated noise assumed)
- `f_A` has sensitivity `Δ`: changing one record shifts `f_A(D)` by at most `Δ`
- `f_B(D) = c · f_A(D)`: client B's output is coupled to client A's via shared model state (`c` = coupling coefficient, `0 ≤ c ≤ 1`)
- **Both outputs are released** and jointly observable by an adversary (e.g., the FL server)

**Adjacency:** `D → D'` removes one record, changing `f_A(D)` by `Δ` (and `f_B(D)` by `c·Δ` via coupling).

**Adversary model:** The adversary observes the full released vector `(output_A, output_B)` jointly, as would be natural for an FL server receiving all client updates before aggregation.

---

### What a User Naturally Does

```python
from dp_accounting.pld import privacy_loss_distribution as pld_lib

# User accounts for mechanism A only (standard API usage)
pld = pld_lib.from_gaussian_mechanism(standard_deviation=sigma, sensitivity=Delta)
delta_stated = pld.get_delta_for_epsilon(eps)
```

This correctly computes the ε-hockey-stick divergence for a **scalar Gaussian
mechanism in isolation**. When `c = 0` (no coupling), this is exact. When `c > 0`,
the correlated release of `output_B` provides additional discriminating information
that is not accounted for.

---

### Reproducible Demonstration (real library API, N=3M samples)

```python
import numpy as np
from numpy.random import multivariate_normal
from dp_accounting.pld import privacy_loss_distribution as pld_lib
import math

np.random.seed(42)
N_MC = 3_000_000
sigma, Delta, eps = 1.0, 1.0, 1.0

# Library (marginal accounting of scalar output A)
pld = pld_lib.from_gaussian_mechanism(standard_deviation=sigma, sensitivity=Delta)
delta_marginal = pld.get_delta_for_epsilon(eps)

def true_delta_joint(c, n=N_MC):
"""True ε-hockey-stick for joint release with coupling c."""
mu_D = np.array([Delta, c * Delta]) # D: f_A=Δ, f_B=c*Δ
mu_Dp = np.zeros(2) # D': f_A=0, f_B=0
Cov = np.eye(2) * sigma**2 # Independent noise
Cov_inv = Cov
X = multivariate_normal(mu_D, Cov, n)
L = 0.5 * (
np.einsum('ij,jk,ik->i', X-mu_Dp, Cov_inv, X-mu_Dp) -
np.einsum('ij,jk,ik->i', X-mu_D, Cov_inv, X-mu_D)
)
return float(np.mean(np.maximum(0, 1 - np.exp(eps - L))))

# c=0: no coupling → matches library exactly (baseline validation)
print(f"c=0.0: true={true_delta_joint(0.0):.4e}, marginal={delta_marginal:.4e}") # ~equal

# c>0: coupling → optimistic marginal estimate
print(f"c=0.3: true={true_delta_joint(0.3):.4e}, marginal={delta_marginal:.4e}") # 1.12x
print(f"c=0.5: true={true_delta_joint(0.5):.4e}, marginal={delta_marginal:.4e}") # 1.34x
print(f"c=0.7: true={true_delta_joint(0.7):.4e}, marginal={delta_marginal:.4e}") # 1.65x
print(f"c=1.0: true={true_delta_joint(1.0):.4e}, marginal={delta_marginal:.4e}") # 2.26x
```

**Verified output (N=3M samples, real library API):**

| Coupling `c` | True H_ε (joint) | Marginal H_ε (library) | Ratio |
|---|---|---|---|
| 0.0 | 1.268e-01 | 1.269e-01 | **1.00x** ← baseline match ✓ |
| 0.1 | 1.286e-01 | 1.269e-01 | 1.01x |
| 0.3 | 1.427e-01 | 1.269e-01 | **1.12x** |
| 0.5 | 1.700e-01 | 1.269e-01 | **1.34x** |
| 0.7 | 2.094e-01 | 1.269e-01 | **1.65x** |
| 1.0 | 2.862e-01 | 1.269e-01 | **2.26x** |

---

### Analytical Characterization

The joint PLRV for this mechanism is:

```
L(x_A, x_B) = Δ/σ² · (x_A + c·x_B) − Δ²(1+c²)/(2σ²)
```

The ε-hockey-stick satisfies:
- `H_ε(0) = delta_marginal` (exact, verified above)
- `dH_ε/dc|_{c=0} = 0` (first-order term cancels; `x_A, x_B` independent at `c=0`)
- `d²H_ε/dc²|_{c=0} > 0` (positive second derivative, verified numerically)

**Local Proposition.** For `c` in a neighbourhood of 0:

```
H_ε(c) = H_ε(0) + A · c² / 2 + O(c⁴), A > 0
```

In particular, `H_ε(c) > H_ε(0)` for all `c ≠ 0` near 0. The marginal estimate
is the strict minimum.

**Note:** Global monotonicity for all `c ∈ [0,1]` is supported numerically (monotone
in all tested cases) but a formal proof for all parameter values is not provided here.

---

### Connection to GCI Sign Theorem

For the special case of correlated noise `(Z_A, Z_B) ~ N(0, [[σ²,ρσ²],[ρσ²,σ²]])`,
the joint privacy event corresponds to a strip geometry in the output space.
The **Sign Theorem** (DOI: [10.5281/zenodo.20078486](https://zenodo.org/records/20078486)) provides an analytical criterion:

```
sign(dH_ε/dρ|_{ρ=0}) = sign(m₁ · m₂) < 0
```

where `m₁ = Δ/2 > 0`, `m₂ = −Δ/2 < 0` are the effective means in the shifted coordinate.
This gives a **first-order** monotone increase for the correlated-noise variant.

---

### Practical Relevance

This gap arises naturally in **synchronous federated learning** where:
- Multiple clients share a common model before each round
- Their gradient functions `f_A, f_B` are implicitly coupled via the shared state
- All client updates are released to the server and jointly observable
- A user may naturally call `from_gaussian_mechanism` per client and compose,
without accounting for the joint observability of correlated outputs

---

### What This Is Not

- This is **not** a mathematical error in the implementation
- This is **not** a claim that correlation alone breaks DP
- This is **not** applicable to the standard single-mechanism case
- The library correctly analyzes the mechanisms it is given

### Suggested Addition

```python
def from_gaussian_mechanism(standard_deviation, sensitivity=1, ...):
"""...existing docstring...

Note:
This function models a **single scalar** Gaussian release.
If the system releases multiple outputs jointly (e.g., correlated
client gradients in federated learning), and an adversary can observe
all outputs together, the true ε-hockey-stick divergence may exceed
this estimate. No current API primitive in this library models
multi-output correlated Gaussian releases.
See https://doi.org/10.5281/zenodo.20078486 for quantitative analysis.
"""
```

---

### References

- Analytical bound (sign theorem + GCI failure):
https://zenodo.org/records/20078486 — independently audited, 100/100, 21/21 prover PASS
- Critical correlation `ρ*(r)`:
https://zenodo.org/records/20078649

Start by locating the Python implementation of pld_lib.from_gaussian_mechanism and its docstring, then read the surrounding API documentation to confirm how scalar Gaussian releases are described. Add a warning that the function models a single scalar release and does not account for jointly observed correlated outputs; done means the limitation and referenced analysis are clearly documented without changing the accounting behavior.

Written by the indexing model from the issue text.