This paper develops a novel differentially private framework to solve convex optimization problems with sensitive optimization data and complex physical or operational constraints. Unlike standard noise-additive algorithms, that act primarily on the problem data, objective or solution, and disregard the problem constraints, this framework requires the optimization variables to be a function of the noise and exploits a chance-constrained problem reformulation with formal feasibility guarantees. The noise is calibrated to provide differential privacy for identity and linear queries on the optimization solution. For many applications, including resource allocation problems, the proposed framework provides a trade-off between the expected optimality loss and the variance of optimization results.

§ INTRODUCTION

Differential privacy <cit.> is a rigorous definition of privacy that quantifies and bounds the risk of disclosing sensitive attributes of datasets used in computations. Differentially private algorithms ensure privacy by introducing a calibrated noise to the inputs, outputs, or objectives of computations. It has been successfully applied to a variety of contexts, including histogram queries <cit.>, census surveys <cit.>, linear regression <cit.> and deep learning <cit.> to name but a few examples. However, their applications to constrained optimization problems remains limited, because it is generally hard to certify the feasibility of differentially private optimization solution.

This paper considers a parametric constrained optimization problems of the form

\begin{align}\label{base_problem}
\minimize{z}\quad& c(z)\quad
\st\quad \mathcal{Z}\triangleq\{z\;|\;Az\leqslant b,\; Gz=d\},
\end{align}

with variables $z \in \mathbb{R}_+^{n}$, convex cost function $c:\mathbb{R}^{n}\mapsto\mathbb{R}$, and convex, compact and non-empty feasible space $\mathcal{Z}$ with parameters $A\in\mathbb{R}^{m\times n}$, $b\in\mathbb{R}^{m}$, $G\in\mathbb{R}^{\ell\times n}$, and $d\in\mathbb{R}^{\ell}$, with $m > 0$ and $\ell > 0$. The paper assumes elements $c, A, b$ and $G$ as public, non-sensitive information about the system design, whereas vector $d=\{d_{i}\}_{i=1}^{\ell}$ contains private, sensitive data of every user $i$, e.g. the customer loads in an electrical power system. In such applications, the feasible space $\mathcal{Z}$ encodes hard operational constraints or physical laws that must be satisfied.

Releasing queries over the solutions of problem (<ref>) may leak information about the sensitive data $d$. For example, in energy network operations, releasing the nodal energy supplies using identity queries, or aggregated supply quantities using sum queries, exposes the allocation of energy demand $d$ across the network <cit.>. Therefore, the goal of this work is to compute such solution $z$ that makes queries over $z$ differentially private, while also being feasible for problem constraints.

While there exist various differential privacy algorithms to solve convex optimization problems, their application to constrained problems is limited because they do not generally guarantee that the privacy-preserving result necessarily satisfies the feasibility conditions. Algorithms based on input perturbation of the sensitive data $d$ modify the feasible space $\mathcal{Z}$ <cit.>,
thus returning an approximate solution to (<ref>) that may not satisfy the original constraints. The output perturbation mechanisms <cit.>, that add noise to the optimization results, generally cannot be certified feasible for $\mathcal{Z}$; see, for example, the impossibility results of <cit.>. The feasibility and near-optimality of the privacy-preserving results can be restored by leveraging the post-processing immunity of differential privacy. This, however, requires solving bilevel optimization problems <cit.>.

This paper addresses these limitations and develops a new framework that provides both privacy and feasibility guarantees for constrained optimization problems. Instead of applying the noise to either the parameters or the results of the optimization, the framework solves a stochastic chance-constrained optimization problem whose solution is used to sample a solution to (<ref>), which guarantees privacy and ensures feasibility with high probability. The approach requires a linear dependency between the optimization variables and the noise <cit.> and reveals a novel connection between differential privacy and stochastic chance-constrained optimization.

hide y axis,
hide x axis,
xtick = , ytick = , ticks = none, no markers, samples = 301,
width = 6cm, height = 8cm,

[name path = x_axis, domain=3:14, color = black, line width = 0.015cm,->,>=stealth] 0*x node[right,above] $z$;
[name path = y_axis_cost, restrict y to domain=0:90,draw=black, line width = 0.015cm,->,>=stealth] (3,x) node[midway,above,rotate=90] $c\braceit{z}$;
[name path = y_axis_dens, restrict y to domain=-90:0,draw=black, line width = 0.015cm,<-,>=stealth] (3,x)node[midway,above,rotate=90] density of  $z(\xi)$;

[name path = cost_feas, domain=5.25:10, smooth, color = black, line width = 0.03cm] 10+5*(x-6)^2;
[name path = cost_infeas, domain=3:5.25, smooth, color = black, line width = 0.015cm, densely dashed] 10+5*(x-6)^2;
[black] coordinates (6,10) node[scale=0.6];
[black,line width = 0.015cm] coordinates (6,10) (5,30) node[above] $\optimal{z}$;
[black] coordinates (8,30) node[scale=0.5];
[black,line width = 0.015cm] coordinates (8,30) (7,50) node[above] $\tilde{z}$;

[name path = feas_region, domain=5.25:10, color = maincolor, line width = 0.05cm] 0*x node[right,above,xshift=-1em,yshift=-0.2em] $\mathcal{Z}$;

[name path = projection_x_star, restrict y to domain=-80:10,draw=gray,line width = 0.015cm, densely dotted] (6,x);
[name path = projection_x,      restrict y to domain=-80:30,draw=gray,line width = 0.015cm, densely dotted] (8,x);
[name path = projection_x_star, restrict y to domain=-5:13,draw=gray,line width = 0.015cm, densely dotted] (5.25,x);
[name path = projection_x_star, restrict y to domain=-11:90,draw=gray,line width = 0.015cm, densely dotted] (10,x);

[name path = opt_loss_up, domain=8:11.5, draw=gray,line width = 0.015cm, densely dotted] 30;
[name path = opt_loss_dw, domain=6:11.5, draw=gray,line width = 0.015cm, densely dotted] 10;
[name path = y_axis_dens, restrict y to domain=10:30,draw=gray, line width = 0.015cm,<->,>=stealth] (11,x)node[midway,below,rotate=90] $\loss$;

[name path = density_feas, black, domain=5.25:10, line width = 0.03cm] 160*laplace(8,1);
[name path = density_infeas_one, black, domain=3:5.25, line width = 0.03cm] 160*laplace(8,1);
[name path = density_infeas_two, black, domain=10:13, line width = 0.03cm] 160*laplace(8,1);

[name path = density_neighbouring, black, domain=3:13, line width = 0.015cm, densely dashed] 160*laplace(6,1);

[maincolor!20] fill between[of=density_feas and feas_region,reverse=true];
[name path = x_axis_, domain=3:5.25, draw=none,] 0*x;
[name path = x_axis__, domain=10:14, draw=none,] 0*x;
[red!30,domain=10:11.5] fill between[of=density_infeas_one and x_axis_,reverse=true];
[red!30,domain=10:11.5] fill between[of=density_infeas_two and x_axis__,reverse=true];

[maincolor, line width = 0.015cm] coordinates (11,-35) node[right] $\mathbb{P}\bracketit{z\in\mathcal{Z}}$;
[red,line width = 0.015cm] coordinates (11,-20) node[right] $\mathbb{P}\bracketit{z\notin\mathcal{Z}}$;

The functioning of the framework is illustrated in the adjacent figure, showing the projections of solutions $z$ onto cost function $c$ and feasible space $\mathcal{Z}$.
Consider the optimal solution $\optimal{z}\;$ returned by problem (<ref>). The output perturbation of $\optimal{z}\;$ results in solutions whose density $\optimal{z}(\xi)$ (dashed line on the bottom of the figure) is prone to lie outside the feasible space $\mathcal{Z}$.
By restricting the optimization variables to be a linear function of the noise, e.g., $z(\xi)=\tilde{z}+f(\xi)$ where $\tilde{z}$ is the expected value of the solution with respect to the random variable distribution and $f(\xi)$ is the linear functional recourse, the stochastic problem optimizes $\tilde{z}$ and $f(\xi)$ providing a new probability density of solution $z$ (solid line). This new density renders any realization of the noise $\xi$ feasible within a prescribed probability $\mathbb{P}\bracketit{\tilde{z}+f(\xi)\in\mathcal{Z}}$, specified by the curator of problem (<ref>). If the functional recourse $f(\xi)$ is made independent from the sensitive data $d$, the framework enjoys both privacy and feasibility guarantees. However, it introduces a trade-off between privacy and the optimality loss with expected value  $\loss=\mathbb{E}[c(z(\xi)) - c(\optimal{z}\;)]$.

The chance-constrained optimization always ensures the satisfaction of the system $G z = d$, which may represent flow conservation constraints and other physical laws that cannot be violated.
This setting, however, restricts the queries to be made only on strict subsets of solution $z$, because the perturbation needs to be redistributed among the variables to guarantee the feasibility of the system. By separating the solution into released and non-released variables, the latter can be optimized to provide trade-offs between the optimality loss and its variance, as well as the trade-offs between the optimality loss and the overall solution variance. The contributions of this work can be summarized as follows:

(1) It develops a novel differentially private framework for the release of identity and sum queries over the solutions of constrained convex programs, using stochastic chance-constrained optimization (Section <ref>).
(2) The released solutions are guaranteed feasible with a high probability. The feasibility guarantees are studied for individual and joint constraint satisfaction, providing higher or lower optimality losses, respectively (Section <ref>).
(3) The framework establishes a trade-off between the expected and the worst-case errors by controlling the variance of the optimization results and the optimality loss (Section <ref>).
(4) On the benchmark energy optimization datasets, the framework is shown to outperform the standard output perturbation algorithm (Section <ref>).

Notation Upper and lower case symbols are used to denote, respectively, matrices and vectors. The indexed notation $A_i$ is used to denote the $i^{\text{th}}$ row vector of matrix $A$. The operator $\text{diag}(a)$ returns the diagonal matrix with entries of vector $a$, and $\text{diagv}(A)$ returns the vector of diagonal elements of matrix $A$. The ceil $\lceil r \rceil$ maps real number $r$ into the least succeeding integer. Notation $\circ$ denotes a Schur product.
$\mathbb{0}$ and $\mathbb{1}$ respectively denote vectors of zeros and ones of proper dimensions.

§ PRELIMINARIES

Across the paper, it is assumed that the optimal solution to problem (<ref>) exists and is unique and that the data $d_i$, contributed by each individual $i$, in $d$ are not correlated. Thus, the problem can be seen as an algorithm
\cM: \mathbb{R}^{\ell}\mapsto\mathbb{R}^{n}
with a unique mapping of datasets $d$ to optimization results. To enable private queries over the optimization results, this work considers a differentially private counterpart $\tilde{\mathcal{M}}$ of $\mathcal{M}$. While the traditional differential privacy definition aims at protecting the participation of an individual data <cit.>,
this work focuses on obfuscating the magnitude $d_i$ associated with participant $i$ of the input vector $d$. To capture this privacy notion, the paper focuses on the indistinguishability framework proposed by <cit.>, which protects the sensitive data of each individual up to some measurable quantity $\alpha > 0$ and defines two neighboring datasets $d, d'$ (written $\sim_\alpha$) as

\begin{align*}
d\sim_{\alpha}d' \Longleftrightarrow \exists i\;\text{s.t.}\;
| d_{i}-d_{i}' | \leqslant \alpha \;\wedge\;d_{j}=d_{j}',\forall j \neq i,
\end{align*}

where $d$ and $d'$ are input vectors to problem (<ref>) and $\alpha$ is a positive real value. Following previous work <cit.>, this relation requires the assumption that the neighboring datasets are feasible for problem (<ref>), which is not restrictive, as only feasible solutions are of interest to release.

Differential privacy requires that the maximum divergence of the algorithm output distributions on neighboring inputs to be bounded by privacy parameters $\varepsilon$ and $\delta$, such that

\begin{align*}
\mathbb{P}\bracketit{\tilde{\mathcal{M}}(d) = O} \leqslant
\mathbb{P}\bracketit{\tilde{\mathcal{M}}(d')= O}\exp(\varepsilon) + \delta
\end{align*}

for a random algorithm $\tilde{\mathcal{M}}$ and any output $O$, where $\mathbb{P}$ denotes the probability over runs of $\tilde{\mathcal{M}}$.
If $\delta =0$, $\tilde{\mathcal{M}}$ is said to be $\varepsilon$-differentially private.

The global sensitivity methods are known to provide differential privacy by augmenting the output of computations with the noise  calibrated to the $\ell_1$- or $\ell_2$-sensitivity.
The $\ell_{p}-$sensitivity

$$\Delta_\alpha \triangleq \max_{d \sim_\alpha d'}\norm{\mathcal{M}(d)-\mathcal{M}(d')}_{p}, \quad p = 1,2,$$

is used to bound the change in the algorithm output induced by any two $\alpha$-indistinguishable inputs.
In many applications of interest,
$G\in\{0,1\}^{\ell\times n}$ is a binary matrix and the domain of datasets $d$ is normalized in $[0,1]$. Thus, $\Delta_{\alpha}$ is directly upper-bounded by $\alpha$.

Let $\text{Lap}(\lambda)^n$ denote the i.i.d. Laplace distribution over $n$ dimensions with $0$ mean and scale $\lambda$. The following ubiquitous result provides an $\varepsilon$-differentially private algorithm <cit.>.

Let $\mathcal{M}$ be an algorithm with $\ell_1$ sensitivity $\Delta_\alpha$ that maps datasets $d$ to $\mathbb{R}^{n}$. The Laplace mechanism $\mathcal{M}(d) + \xi$, with $\xi\sim\text{Lap}(\nicefrac{\Delta_\alpha}{\varepsilon})^{n}$, attains $\varepsilon$-differential privacy.

§ INTERNALIZING GLOBAL SENSITIVITY METHODS INTO CONSTRAINED OPTIMIZATION

The direct application of Theorem <ref> to the optimal optimization solution may produce a result that violates the problem constraints. This section introduces a suitable transformation of problem (<ref>) into a stochastic chance-constrained problem that internalizes the global sensitivity methods to establish both privacy and constraint feasibility guarantees. This section first discusses the Laplace mechanism and then extends the results to the Gaussian mechanism.

Consider a random perturbation $\xi\in\mathbb{R}^{p}$ calibrated to Laplace distribution $\mathbb{P}_{\xi} = \text{Lap}(\nicefrac{\Delta_{\alpha}}{\varepsilon})^p$ for some arbitrary dimension $p$, and assume that the solution $z$ depends on the realization of $\xi$ as

\begin{align}\label{rand_sol}
z(\xi)= \tilde{z} + f(\xi) = \tilde{z} + Z\xi,
\end{align}

where $\tilde{z}\in\mathbb{R}^{n}$ is the expected value of the solution with respect to distribution $\mathbb{P}_{\xi}$ and $Z\xi$ is the linear functional recourse with recourse decision $Z\in\mathbb{R}^{n\times p}$, which is used to adjust the expected solution to any realization of $\xi$. Therefore, any query made over $z(\xi)$ will constitute the expected and random components. To provide privacy guarantees, the random component is required to be  independent from data $d$. This can be achieved by enforcing additional, query-specific, constraints $\mathcal{Q}$ on the recourse $Z$. While the framework can accommodate the general class of linear queries over solution $z(\xi)$, for ease of presentation, this work focuses on identity and sum queries.

This query releases a specified subset of solution $z(\xi)$.
Consider a random perturbation $\xi\in\mathbb{R}^{n}$ and a diagonal matrix $I\in\mathbb{R}^{n\times n}$, such that

\begin{align*}
\xi_{i} = \left\{
\begin{array}{l}
\text{Lap}(\nicefrac{\Delta_{\alpha}}{\varepsilon}), \\
\end{array} \right.
I_{ii} = \left\{
\begin{array}{rl}
1, & \text{if $i^{\text{th}}$ element of vector $z(\xi)$ is subject to release},\\
0, & \text{otherwise.}
\end{array} \right.
\end{align*}

The identity query release is thus $Iz(\xi) = I\tilde{z} + IZ\xi$ for $Z\in\mathbb{R}^{n\times n}$.
The random component $IZ\xi$ is made independent from the dataset $d$ if the recourse decision $Z$ is constrained as follows

Z\in\mathcal{Q}(I)\triangleq\{Z\;|Z\circ I = I,\;Z\circ(I-\text{diagv}(I)\mathbb{1}_{n}^{\top}) = \mathbb{0}_{n\times n}\},

which yields $Iz(\xi)=\left\{
\begin{array}{ll}
\tilde{z}_{i} + \xi_{i}, & \text{if $I_{ii}=1$,}\\
0, & \text{otherwise.}
\end{array} \right.$

This query releases $p$ sum statistics over non-intersecting subsets of $z(\xi)$.  Consider a random perturbation $\xi\in\mathbb{R}^{p}$ and a matrix $S\in\mathbb{R}^{p\times n}$, such that

\begin{align*}
S_{ij} = \left\{
\begin{array}{rl}
1, & \text{if element $z_{j}(\xi)$ participates in sum statistic $i$},\\
0, & \text{otherwise.}
\end{array} \right.
\end{align*}

The sum query releases $Sz(\xi)=S\tilde{z} + SZ\xi$ for $Z\in\mathbb{R}^{n \times p}$.
The random component $SZ\xi$ is made independent from the dataset $d$ if the recourse decision $Z$ is constrained as follows

Z\in\mathcal{Q}(S)\triangleq\{Z\;|SZ = \text{diag}(\mathbb{1}_{p})\},

which yields $S z(\xi)= S \tilde{z} + \xi$.

To produce random solutions to (<ref>), function (<ref>) is optimized using the following stochastic program

\begin{align}
\minimize{\tilde{z},Z\in\mathcal{Q}}\quad &\mathbb{E}^{\mathbb{P}_\xi}\bracketit{c(\tilde{z} + Z\xi)}\label{cc_problem_v2_obj}\\
\st\quad &\mathbb{P}_{\xi}\bracketit{A(\tilde{z} + Z\xi)\leqslant b}\geqslant 1-\eta,\label{cc_problem_v2_cc}\\
&G(\tilde{z} + Z\xi)=d \quad\mathbb{P}_{\xi}\text{-a.s.},\label{cc_problem_v2_asc}
\end{align}

which optimizes $\tilde{z}$ and $Z$ by anticipating all realizations of the random variable $\xi$. This problem minimizes the expected value of the convex cost function (<ref>) with respect to the random variable $\xi$.

The problem constraints are given by a set of probabilistic constraints. The joint chance constraint (<ref>) requires the satisfaction of the inequality constraints with a prescribed probability $1-\eta$, specified by the curator of problem (<ref>). The almost sure constraint (<ref>) requires the equality constraints to hold with probability 1. Note that, if problem (<ref>) is infeasible, it follows that the privacy parameters $\alpha$ and $\varepsilon$ are too strong for the feasibility requirement $\eta$.

As the recourse of problem (<ref>) amounts to finitely-dimensional linear functions, objective function (<ref>) and chance constraint (<ref>) admit computationally tractable reformulations <cit.> (additional details will be given in Section <ref>).
The almost sure constraint (<ref>) includes a random variable and, therefore, satisfying it is computationally intractable. However, it can be equivalently reformulated using the following set of equations:

\begin{align}\label{balance}
G\tilde{z} = d,\quad GZ = \mathbb{0}.
\end{align}

If variables $\tilde{z}$ and $Z$ are subject to (<ref>), their optimal solution satisfies the equality constraint (<ref>) for any realization of $\xi$. The structural properties of $G$ restrict the set of potential queries and a query is said to be implementable if there exists $Z\in\mathcal{Q}$ such that $GZ=\mathbb{0}$ holds.

[Flow conservation constraint]Assume $G\in\mathbb{R}^{\ell\times n}$ represents the incidence matrix of a fully connected graph $\mathcal{G}$ (its rank is $n-1$) and that $G z = d$ represents a flow conservation constraint. Consider an identity query $Iz(\xi) = I\tilde{z} + IZ\xi$ and $Z\in \mathcal{Q}(I)$ as in Definition <ref>. The identity query is implementable if $\Tr[I] < n$, i.e., not all elements of $Z$ are constrained by $\mathcal{Q}(I)$ and $GZ = \mathbb{0}$ holds.

The constraint $GZ = \mathbb{0}$ plays a critical role: it balances the perturbation between the released and non-disclosed variables. As a result, the query should leave enough degree of freedom to satisfy the equality constraint.

This limitation is solely induced by the need to preserve the satisfaction of the equality constraint and it is not seen as a limiting factor for many applications

(see Section <ref>).

Private identity query (PIQ) algorithm The procedure is summarized in Algorithm <ref>, which takes as inputs the dataset $d$, the $\ell_{1}$-sensitivity of the identity query, the privacy $\varepsilon$ and feasibility $\eta$ requirements, the known covariance $\Sigma$, and the query specification $I$. Upon receiving the optimal chance-constrained solution (line 2), the algorithm draws a sample from the Laplace distribution (line 3) and computes a $(1-\eta)$-feasible solution for problem (<ref>) (line 4).
The algorithm returns an $\varepsilon$-differentially private identity query which satisfies problem (<ref>) constraints with probability $(1-\eta)$.

Algorithm <ref> is $\varepsilon$-differentially private, i.e.,

\begin{align*}
\mathbb{P}\bracketit{I\braceit{\optimal{z}(d) + \optimal{Z}(d)\xi} = O}\leqslant\mathbb{P}\bracketit{I\braceit{\optimal{z}(d') + \optimal{Z}(d')\xi} = O}\text{exp}(\varepsilon)
\end{align*}

for any two $\alpha-$neighboring datasets $d$ and $d'$ and output solutions $O$.

Private sum query (PSQ) algorithm The procedure is summarized in Algorithm <ref>, which differs from Algorithm <ref> by the query specification $S$ and the noise dimension.

Algorithm <ref> is $\varepsilon$-differentially private, i.e.,

\begin{align}
\mathbb{P}\bracketit{S(\optimal{z}(d)+\optimal{Z}(d))\xi = O}\leqslant\mathbb{P}\bracketit{S(\optimal{z}(d')+\optimal{Z}(d'))\xi = O}\text{exp}(\varepsilon)
\end{align}

for any two $\alpha-$neighboring datasets $d$ and $d'$ and output solutions $O$.

Input: $d,\Delta_\alpha, \varepsilon, \eta, I$

$(\optimal{\tilde{z}}\;, \optimal{Z}\;) \gets $
Solve (<ref>) for $Z\in\mathcal{Q}(I)$

$\hat{\xi} \gets$ Sample from $\text{Lap}(\nicefrac{\Delta_\alpha}{\varepsilon})^{n}$

Compute solution to (<ref>) as $\hat{z}=\optimal{\tilde{z}}\;+\optimal{Z}\hat{\xi}$

Release: $I\hat{z}$

Private identity query (PIQ)

Input: $d,\Delta_\alpha, \varepsilon, \eta,  S$

$(\optimal{\tilde{z}}\;, \optimal{Z}\;) \gets $
Solve (<ref>) for $Z\in\mathcal{Q}(S)$

$\hat{\xi} \gets$ Sample from $\text{Lap}(\nicefrac{\Delta_\alpha}{\varepsilon})^{p}$

Compute solution to (<ref>) as $\hat{z}=\optimal{\tilde{z}}\;+\optimal{Z}\hat{\xi}$

Release: $S\hat{z}$

Private sum query (PSQ)

In addition to releasing a privacy-preserving answer with a probabilistic feasibility certificate, the proposed framework also allows to verify the feasibility of the sampled solution $\hat{z}$ without incurring an additional privacy loss. Since the equality constraint holds due to (<ref>), it is sufficient to verify the feasibility of constraint $A\hat{z}\leqslant b$ without accessing the original data $d$. The operation is private by post-processing immunity of differential privacy <cit.>.

Furthermore, since formulation (<ref>) is independent from the distribution of the noise, the framework can accommodate other global sensitivity methods. In particular, the following result holds.
Let $\delta,\varepsilon\in(0,1)$ and let $\Delta_{\alpha}^{2}$ be the $\ell_{2}-$sensitivity.
Algorithms <ref> and <ref> that calibrate $\xi\in\mathcal{N}(0,\sigma^2)$ to the Gaussian distribution with $\sigma\geqslant \Delta_{\alpha}^{2} \sqrt{2\ln(1.25/\delta)}/\varepsilon$ satisfy $(\varepsilon,\delta)$- differential privacy.

The relation between the feasibility requirement $\eta$ and the privacy parameters $\varepsilon$ and $\delta$ is implicit in the formulation of the chance-constrained problem: the variance of the noise affects the ability to satisfy the problem constraints within the feasibility requirement and vice-versa. To render this relation explicit, Appendix <ref> discusses a version of Algorithms <ref> and <ref> that iterates lines 3 to 5 an optimal number $T$ of times to guarantee the release of a feasible solution with probability $1-\mu$, for some $0 < \mu < 1$.

Given feasibility requirement $\eta$, privacy parameter $\varepsilon/T$, and value $0 < \mu < 1$, the iterative variants of Algorithms <ref> and <ref> return an $\varepsilon$-differentially private solution that is feasible with probability at least $1-\mu$ within $T = \lceil \frac{\log(\mu)}{\log(\eta)} \rceil$ steps.

§ REFORMULATIONS AND FEASIBILITY GUARANTEES

The optimization problem (<ref>) is intractable because it constitutes the optimization of a random variable. However, due to the convexity assumption on (<ref>), linear functional recourse and known distribution of $\xi$, problem (<ref>) admits tractable reformulations. There are several avenues to reformulate the joint chance constraint (<ref>) with different degrees of conservatism in terms of expected optimality loss <cit.>. This work provides a conservative joint constraint satisfaction guarantee, using a sample approximation, and a less conservative individual constraint satisfaction guarantee, using an analytic reformulation. The objective function is reformulated as follows.

Objective function reformulation Consider a quadratic cost function with first- and second-order coefficients $c_{1}\in\mathbb{R}^{n}$ and $c_{2}\in\mathbb{R}^{n}$, and a diagonal covariance matrix $\Sigma=\mathbb{E}[\xi\xi^{\top}]$ with diagonal elements being equal to $\lambda=\Delta_\alpha/\varepsilon$. Then, the objective function (<ref>) reformulates as

\begin{align*}
&\mathbb{E}^{\mathbb{P}_\xi}\bracketit{c_{1}^{\top}(z+Z\xi) + (z+Z\xi)^{\top}\text{diag}(c_{2})(z+Z\xi)}
= c_{1}^{\top}z + z^{\top}\text{diag}(c_{2})z + \Tr\bracketit{Z^{\top}\text{diag}(c_{2})Z\Sigma},
\end{align*}

which follows from the zero-mean distribution $\mathbb{P}_{\xi}$ and the fact that $\mathbb{E}[\xi\xi^{\top}]=\Sigma.$ Notice that for the affine cost functions, the analytic reformulation of (<ref>) reduces to $c_{1}^{\top}z$.

Sample approximation This approximation substitutes the chance constraint (<ref>) with a finite number of deterministic constraints, each enforced on a specific realization of random perturbation <cit.>. This work invokes the sample approximation method from <cit.>, which enforces (<ref>) on the vertices of the rectangular sample set extracted from distribution $\mathbb{P}_{\xi}$, i.e., for $\xi\in\mathbb{R}^{p}$

\begin{align}\label{sample_ref}
\mathbb{P}_{\xi}\bracketit{
\braceit{A(z+Z\xi)\leqslant b}}\geqslant 1-\eta
\quad\equiv\quad
Az \leqslant b - AZ\hat{\xi}^{v},\quad\forall v=1,\dots,2^{p},
\end{align}

where $\hat{\xi}^{v}\in\mathbb{R}^{2^p}$ is the $v^{\text{th}}$ vertex of the extracted sample set. <cit.> show that the joint constraint satisfaction is attained if the number of samples $N$ from $\mathbb{P}_{\xi}$ is properly chosen.

The equivalence (<ref>) holds with confidence $(1-\beta)$ if the rectangular set is built upon $S$ samples extracted from $\mathbb{P}_{\xi}$, with $N$ at least as much as

N \geqslant \left\lceil
\frac{1}{\eta}\frac{e}{e-1}\braceit{2p-1+\text{ln}\frac{1}{\beta}}
\right\rceil.

Finally, notice that this approximation requires an additional input $\beta$ to Algorithms <ref> and <ref> to accommodate the confidence level of the data curator.

Analytic reformulation

The joint chance constraint (<ref>) can be rewritten as a union of individual chance constraints. For some vector $\overline{\eta}\in\mathbb{R}_{+}^{m}$ of individual constraint violation probabilities, the individual chance constraints can be reformulated exactly using second-order cone constraints <cit.>:

\begin{align}\label{analytic_ref}
%     \textstyle
A_{i}z \leqslant b_{i} - f(1-\overline{\eta}_{i})\norm{A_{i}Z\Sigma^{\nicefrac{1}{2}}}_{2}, \quad\forall i=1,\dots,m,
\end{align}

where $f(1-\overline{\eta}_{i})$ is a distribution-dependent safety parameter and $\Sigma^{\nicefrac{1}{2}}$ is the lower triangular matrix resulting from the Cholesky factorization of $\Sigma$. For any symmetric and unimodal distribution of $\xi$, $f(1-\overline{\eta}_{i})$ amounts to $(\nicefrac{2}{9\overline{\eta}_{i}})^{\nicefrac{1}{2}}$ if $0\leqslant\overline{\eta}_{i}\leqslant\nicefrac{1}{6}$ (see the result from <cit.>). For the Gaussian distribution of $\xi$, $f(1-\overline{\eta}_{i})$ amounts to the inverse CDF of the standard Gaussian distribution at $(1-\overline{\eta}_{i})-$quantile <cit.>. Observe that, for the fixed parameter $\overline{\eta}_{i}$, the last term in the right-hand side of (<ref>) is a safety margin, which reduces the feasible space of the original problem (<ref>) to guarantee individual constraint feasibility for $(1-\overline{\eta}_{i})$ realizations of $\xi$.

If $\mathbb{1}^{\top}\overline{\eta}\leqslant\eta$ holds, the individual chance constraints guarantee joint constraint satisfaction probability $\eta$. Yet, finding the optimal value $\overline{\eta}$ is an NP- hard problem <cit.>. Section <ref>, shows that, for small problem instances, the choice $\overline{\eta}=\eta$ results in the desired joint constraint satisfaction while providing a significantly less conservative solution than the sample approximation.

§ VARIANCE-AWARE DIFFERENTIALLY PRIVATE ALGORITHMS

For many systems governed by the solution of problem (<ref>), e.g. energy networks, it is important to control the impact of the differentially private solutions on the optimality loss (e.g., extra supply cost) and the variance of the state variables (e.g., supply and flow allocations). This section extends Algorithms <ref> and <ref> to provide a minimal variance solution without affecting the privacy guarantees.

Minimal variance of optimality loss
The chance-constrained problems of Algorithms <ref> and <ref> optimize against the expected value of the cost function. Its solution provides the estimate of the expected optimality loss relative to the solution of problem (<ref>). The worst-case outcome of the optimality loss, however, may significantly exceed the expected value. A trade-off between the expected and worst-case outcomes can be attained by controlling the variance of the optimality loss.

As the value of the cost function of problem (<ref>) is deterministic, it is sufficient to control the variance of (<ref>) to attain the desired result. For a linear cost function, the variance admits a convex expression $\text{Var}\bracketit{c_{1}^{\top}(z + Z\xi)}=\Tr\bracketit{Z^{\top}\text{diag}(c_{1})\text{diag}(c_{1})Z\Sigma}$ in recourse variable $Z$. Therefore, it can be minimized by optimizing,  instead, the following objective function

\begin{align}\label{var_problem_loss}
\minimize{z,Z\in\mathcal{Q}}\quad& (1-\varphi)\mathbb{E}^{\mathbb{P}_\xi}\bracketit{c_{1}^{\top}(z + Z\xi)} + \varphi\norm{\Sigma^{\nicefrac{1}{2}}Z^{\top}c_{1}}_{2},
\end{align}

which optimizes the trade-off between the expected value and the standard deviation of the cost function for some factor $\varphi\in[0,1]$.
Thus, varying the factor $\varphi$ establishes a Pareto frontier between the optimality loss and its variance. Since the recourse decision $Z$ is subject to query-specific constraints $\mathcal{Q}$, the results of Theorems <ref> and <ref> hold. Finally, the variance of the non-affine cost functions does not permit convex formulations and is not considered in this paper.

Minimal variance of optimization variables  The variance of the  optimization solution $z\braceit{\xi}=z+Z\xi$ admits a convex expression $\text{Var}\bracketit{z\braceit{\xi}} = \Tr\bracketit{Z^{\top}\Sigma Z}$ in $Z$. Therefore, it can be controlled by optimizing the recourse decision $Z$ using the following objective function

\begin{align}\label{var_problem_variables}
\minimize{z,Z\in\mathcal{Q}}\quad& (1-\varphi)\mathbb{E}^{\mathbb{P}_\xi}\bracketit{c(z + Z\xi)} + \varphi\norm{Z\Sigma^{\nicefrac{1}{2}}\mathbb{1}}_{2},
\end{align}

which finds the optimal trade-off between the expected cost and the standard deviation of the optimization variables for some factor $\varphi\in[0,1]$. Since the optimal recourse is still guided by the query-specific constraints, the privacy guarantees provided by Theorems <ref> and <ref> are preserved.

§ EXPERIMENTS

Problem description
The proposed framework is applied to the energy resource allocation problem using a set of benchmark networks from <cit.>. The problem goal is to compute the cost-optimal supply allocations across the network to satisfy nodal demands while respecting the supply and network limits. The problem is described by an undirected graphs $\mathcal{G}(N,E)$ with a set of nodes $N$ and a set of edges $E$, connecting those nodes. The graph typology is represented by the weighted Laplacian matrix $B$ formed from non-negative edge weights $\beta\in\mathbb{R}_{+}^{\setcard{E}}$. The nodal supply  $p\in\mathbb{R}_{+}^{\setcard{N}}$ is allocated in the network to meet nodal demand $d\in\mathbb{R}_{+}^{\setcard{N}}$. The flow along the edges is modeled considering a vector of nodal potentials $\theta\in\mathbb{R}^{\setcard{N}}$, their difference is proportional to the network flows, i.e., the flow in edge $\ell$ amounts to $f_{\ell}(\theta)=\beta_{\ell}(\theta_{\se{\ell}}-\theta_{\re{\ell}}),\forall \ell\in E$, with operators $\se{\ell}$ and $\re{\ell}$ returning the sending and receiving nodes of edge $\ell$, respectively. Finally, the nodal supply incurs costs computed by function $c:\mathbb{R}^{\setcard{N}}\mapsto\mathbb{R}$. This allocation problem gives rise to the following optimization

\begin{align}
\minimize{p,\theta} \quad & c(p)\\
\st \quad& B\theta = p - d \label{problem_det_bal}\\
&\underline{p} \leqslant p \leqslant \overline{p} \label{problem_det_supply}\\
&\underline{f} \leqslant f(\theta) \leqslant \overline{f}. \label{problem_det_flow}
\end{align}

The objective function minimizes the total supply cost, while the equality constraint balances nodal demand, supply, and net flow injection. The inequality constraints respect the minimum and maximum nodal supply and the network flow limits $\underline{p},\overline{p}\in\mathbb{R}_{+}^{\setcard{N}}$, $\underline{f}\in\mathbb{R}_{-}^{\setcard{E}},\overline{f}\in\mathbb{R}_{+}^{\setcard{E}}$. A rearrangement of the terms in (<ref>)-(<ref>) makes the problem representable in the form expressed by problem (<ref>), thus its chance-constrained counterpart is achieved as detailed in Sections <ref> and <ref>.

The experiments concern the identity and sum queries made over the subset of nodal supplies $p$ and make use of the $\ell_{1}-$sensitivity $\Delta_{\alpha}$ of vector $p$ on the two $\alpha-$indistinguishable datasets $d$ and $d'$. Consider the optimal supply allocations $\optimal{p}\;$ and $\optimal{p}'$ obtained, respectively, on datasets $d$ and $d'$.

$\Delta_{\alpha}=\norm{\optimal{p}-\optimal{p}'}_{1} \leqslant \alpha.$

The equality constraint (<ref>) requires the balance between the total supply and total demand. By construction, $\sum_{i\in N}[B\optimal{\theta}\;]_{i} = 0,$ thus from (<ref>) it follows

\begin{align*}
&\sum_{i\in N} \optimal{p}_{i} - \sum_{i\in N} d_{i} = 0 \\
&\sum_{i\in N} \optimal{p}_{i}^{\prime} - \sum_{i\in N} d_{i}^{\prime} = 0,
\end{align*}

\sum_{i\in N}\optimal{p}_{i} - \sum_{i\in N}\optimal{p}_{i}^{\prime} = \sum_{i\in N} d_{i} - \sum_{i\in N} d_{i}^{\prime} \leqslant \alpha,

because the datasets $d$ and $d'$ differ by at most $\alpha$ in one entry, i.e., $\norm{d-d'}_{1}\leqslant\alpha$.
Therefore, $\Delta_{\alpha}=\norm{\optimal{p}-\optimal{p}'}_{1} \leqslant \alpha$.

Experimental setup
The experiments are organized as follows. For every network instance, the variable limits are fixed, while cost coefficients and nodal demands are i.i.d. drawn from the following uniform distributions  $c_{1}\sim U(1,3)$, $c_{2}\sim U(\nicefrac{1}{10},\nicefrac{3}{10})$, and $d\sim U(\nicefrac{1}{2},1).$ The results are thus reported for 100 independent simulation runs. The identity and sum queries are made over an arbitrary set of 30% of nodal supplies, which is sampled at every simulation run. The privacy loss parameter $\varepsilon$ is set to $1$ and the indistinguishability parameter $\alpha$ is set to $0.1$ for identity queries and $0.5$ for sum queries. As $\ell_{1}$-sensitivity $\Delta_{\alpha}$ is bounded by $\alpha$, random perturbations thus obey the Laplace distribution $\text{Lap}(\alpha)$. The feasibility requirements for the joint and individual constraint satisfaction are set uniformly at $\eta=\overline{\eta}=2.5\%$, and the out-of-sample empirical constraint violation probability is obtained for 1000 samples at every simulation run.

Privacy-preserving algorithms The algorithm abbreviations SQ   are appended by s   to indicate whether the chance constraints are reformulated, respectively, analytically or by samples (Section <ref>). They are compared with the output perturbation P  solution is said to be feasible if problem (<ref>) returns a feasible solution for the fixed solution of the _FG.

Algorithm comparison The algorithms are compared in terms of their ability to release private queries while satisfying the feasibility requirement. Table <ref> summarizes the results for identity query answers obtained on several networks differing by the number of variables ($n$) and constraints ($|\mathcal{Z}|$). The results indicate that the P  and SQ-a  algorithm releasing $9$ sum statistics for the 8_ieee case and for the various assignments of the trade-off parameter $\varphi\in[0,1]$. The left plot shows that independently of cost sparsity, the algorithm can produce a differentially private output at zero variance of the optimality loss. The right plot demonstrates the drastic reduction of the overall solution variance $\text{Var}\bracketit{z(\xi)}=\Tr\bracketit{Z^{\top}\Sigma Z}$, almost to the $\Tr[\Sigma]$ total variance of the $9$ random perturbations. Both results, however, require larger conservatism of the solution in terms of optimality loss.

Identity query summary for 100 network data samples

max width=0.9

1c3*Case ID     3*$n\times|\mathcal{Z}|$     6cEmpirical constraint violation $\mathbb{P}[\hat{z}\notin\mathcal{Z}]$ [%]     4cOptimality loss $\Delta c$ [%]

(r)3-8 (r)9-12
1c          2c OP      2cPIQ-a     2cPIQ-s     2cPIQ-a     2cPIQ-s

1c          mean     std     mean     std     mean     std     mean     std     mean     std

.07     0.64     0.34     0.27     0.31     3.42     3.55     5.72     7.64

2.24     26.85     1.48     0.78     0.27     0.25     1.55     1.33     3.45     2.88

2     0.49     0.35     2.17     0.81     4.7     1.68

9     1.91     7.17     1.50     1.28     1.06     2.4     0.78     5.51     5.06

8_ieee     236$\times$728     99.99     0.02     14.35     2.06     1.51     0.47     2.46     0.56     4.89     1.20

Sum query summary for 100 network data samples

max width=0.72

[t]3mm3*[origin=c]90queries #     6cEmpirical constraint violation  $\mathbb{P}[\hat{z}\notin\mathcal{Z}]$ [%]     4cOptimality loss $\Delta c$ [%]

2c OP      2cPLQ-a     2cPLQ-s     2cPLQ-a     2cPLQ-s

mean     std     mean     std     mean     std     mean     std     mean     std

1     0.00     0.00     0.00     0.00     0.00     0.00     0.00     0.00     0.00     0.00

3     0.62     0.03     0.00     0.00     0.00     0.00     0.07     0.04     0.12     0.11

6     16.48     20.2     1.39     0.84     0.2     0.24     0.58     0.32     1.81     0.82

9     59.83     23.6     4.94     1.11     1.00     0.92     2.36     1.01     13.35     7.31

xlabel near ticks,
ylabel near ticks,
ytick=0, 2.5, 5, 7.5,
y tick label style=/pgf/number format/.cd,fixed,fixed zerofill,precision=1,/tikz/.cd,
label style=font=,
tick label style=font=,
legend pos = south west,
legend style=font=,
legend cell align=left,fill=none,xshift=4.55cm,yshift=0cm,
width=5.6cm, height=3.5cm]
[black, line width=0.15mm, mark=pentagon*,every mark/.append style=solid, fill=gray, thick,mark size=1.5pt,
mark repeat*=20, smooth]
table [x index = 1, y index = 2] data_plots/summary_loss.dat;
[black, line width=0.15mm, mark=*,every mark/.append style=solid, fill=gray, thick,mark size=1.5pt,
mark repeat*=20, smooth]
table [x index = 3, y index = 4] data_plots/summary_loss.dat;
[black, line width=0.15mm, mark=otimes*,every mark/.append style=solid, fill=gray, thick,mark size=1.5pt,
mark repeat*=20, smooth]
table [x index = 5, y index = 6] data_plots/summary_loss.dat;
[black, line width=0.15mm, mark=triangle*,every mark/.append style=solid, fill=gray, thick,mark size=1.5pt,
mark repeat*=20, smooth]
table [x index = 7, y index = 8] data_plots/summary_loss.dat;
[black, line width=0.15mm, mark=square*,every mark/.append style=solid, fill=gray, thick,mark size=1.5pt,
mark repeat*=20, smooth]
table [x index = 9, y index = 10] data_plots/summary_loss.dat;
[->,>=stealth, gray] coordinates
(25,7.5) (29,1.5) node [midway, above, sloped] $\varphi$ increases;

xshift = 8cm,
xlabel near ticks,
ylabel near ticks,
x tick label style=/pgf/number format/1000 sep= ,
label style=font=,
tick label style=font=,
legend pos = south west,
legend style=draw=none, font=,
legend cell align=left,fill=none,
width=5.6cm, height=3.5cm]
[black, line width=0.15mm, mark=pentagon*,every mark/.append style=solid, fill=gray, thick,mark size=1.5pt,
mark repeat*=20, smooth]
table [x index = 1, y index = 2] data_plots/summary_var_sol.dat;
[black, line width=0.15mm, mark=*,every mark/.append style=solid, fill=gray, thick,mark size=1.5pt,
mark repeat*=20, smooth]
table [x index = 3, y index = 4] data_plots/summary_var_sol.dat;
[black, line width=0.15mm, mark=otimes*,every mark/.append style=solid, fill=gray, thick,mark size=1.5pt,
mark repeat*=20, smooth]
table [x index = 5, y index = 6] data_plots/summary_var_sol.dat;
[black, line width=0.15mm, mark=triangle*,every mark/.append style=solid, fill=gray, thick,mark size=1.5pt,
mark repeat*=20, smooth]
table [x index = 7, y index = 8] data_plots/summary_var_sol.dat;
[black, line width=0.15mm, mark=square*,every mark/.append style=solid, fill=gray, thick,mark size=1.5pt,
mark repeat*=20, smooth]
table [x index = 9, y index = 10] data_plots/summary_var_sol.dat;
[->,>=stealth, gray] coordinates
(23,120) (26,20) node [midway, above, sloped] $\varphi$ increases;
[red,dashed] coordinates
(16,5) (35,5) node [midway, above, sloped, xshift=1.0cm] $\Tr[\Sigma]=4.5$;

Trade-offs: expected value of optimality loss vs. variance (left) and vs. variance of optimization solution (right). The results are given for $c_{1} \!\sim\! U[1,\overline{c}]$ and averaged over 100 runs.

§ RELATED WORK

There is a large body of work on differentially private algorithms for convex optimization in the context of empirical risk minimization (ERM) problems. Output perturbation algorithms <cit.> focus on adding the noise to the optimization results. Objective perturbation algorithms <cit.> perturb the optimization objective and perform well for smooth loss functions. Exponential sampling algorithms <cit.> rely on an evaluation function to select a candidate output, while achieving $(\varepsilon, 0)$-differential privacy that may be difficult to implement due to the exponential nature of evaluation function. Finally, noisy stochastic gradient descent (SGD) algorithms <cit.> provide a privacy-preserving version of SGD that can be combined with accountant methods to provide tight bounds. All these algorithms, however, are meant for a particular class of unconstrained or regularized convex optimization problems and do not focus on reporting solutions that must satisfy problem constraints.

The contributions on differentially private constrained convex optimization for generic decision-making problems are much more sparse. <cit.> studied differential privacy in combinatorial optimization problems and derived information-theoretic bounds on the task utility. <cit.> proposed to solve linear programs privately using a differentially private variant of the multiplicative weights mechanism. <cit.> focused on a particular class of convex optimization problems whose objective function is piecewise affine, with the possibility of including linear inequality constraints.
<cit.> proposed a private data-release mechanism relying on projections to restore the feasibility of the violated constraints due to input perturbation.

Finally, <cit.> developed a differentially private algorithm for a class of linear programs that solely include the inequality constraints whose right-hand side contains sensitive data. The work relies on the input perturbation of the inequality right-hand sides to achive $(\varepsilon, \delta)-$differential privacy.

There are also differential privacy proposals for the distributed convex optimization. A privacy-preserving version of the alternating direction method of multipliers <cit.> has been studied in the context of the unconstrained ERM problem  <cit.> and constrained energy resource allocation problem <cit.>. <cit.> proposed a private distributed projected gradient descent algorithm for constrained convex optimization problems. This collection of work, however, minimizes the privacy leakage by acting on the information exchanged by agents during the coordination process and does not provide privacy guarantees for the release of optimization solution.

§ CONCLUSION

The paper proposed a novel framework to release privacy-preserving solutions of constrained convex optimization problems that contain complex feasibility constraints. The framework relies on a combination of differential privacy and stochastic optimization theory and provides the foundations for two algorithms answering privacy-preserving identity and sum queries over the optimization solutions. The feasibility guarantees were studied for both individual and joint constraint satisfaction and the paper examined the trade-off between the expected and the worst-case errors by controlling the variance of the solutions and the optimality loss. Finally, the proposed framework was shown to outperform standard output perturbation algorithms on several energy benchmark networks.

§ PROOF OF THEOREM <REF>

Without loss of generality, consider that the identity query $Iz(\xi)=I\tilde{z}+IZ\xi$ requires releasing first $k$ items of $z(\xi)$, such that the diagonal matrix $I$ can be described as

\begin{align*}
I =
\begin{bmatrix}
\text{diag}(\mathbb{1})_{k\times k} & \\
& \text{diag}(\mathbb{0})_{n-k\times n-k}
\end{bmatrix},
\end{align*}

the perturbation vector $\xi$ as

\begin{align*}
\xi =
\begin{bmatrix}
\xi_{1}, \dots, \xi_{k}, \mathbb{0}_{n-k}^{\top}
\end{bmatrix}^{\top},
\end{align*}

and an arbitrary identity outcome $O$ as

\begin{align*}
O =
\begin{bmatrix}
O_{1}, \dots, O_{k}, \mathbb{0}_{n-k \times n-k}^{\top}
\end{bmatrix}^{\top}.
\end{align*}

Denote the optimal solution of the chance-constrained problem (<ref>)  by $\optimal{\tilde{z}}$ and $\optimal{Z}$ .
It needs to be shown that the ratio of probabilities
that the algorithm returns the same outcome $O$ on two $\alpha-$indistinguishable input datasets $d$ and  $d'$ is bounded by a constant $\text{exp}\braceit{\varepsilon}$:

\begin{align*}
\mathbb{P}\bracketit{I\braceit{\optimal{\tilde{z}}(d) + \optimal{Z}(d)\xi} = O} /\;
\mathbb{P}\bracketit{I\braceit{\optimal{\tilde{z}}(d') + \optimal{Z}(d')\xi} = O}
\leqslant \text{exp}\braceit{\varepsilon}.
\end{align*}

It follows that:

\begin{align*}
&\mathbb{P}\bracketit{I\optimal{\tilde{z}}(d) + I\optimal{Z}(d)\xi = O} /\;
\mathbb{P}\bracketit{I\optimal{\tilde{z}}(d') + I\optimal{Z}(d')\xi = O}
\\
\frac{\mathbb{P}\bracketit{\begin{bmatrix}\optimal{\tilde{z}}_{1}(d) \\ \vdots \\ \optimal{\tilde{z}}_{k}(d) \\ \mathbb{0}_{n-k}\end{bmatrix} + \begin{bmatrix}
\xi_{1} \\ \vdots \\ \xi_{k} \\ \mathbb{0}_{n-k}
\end{bmatrix} =
\begin{bmatrix}
O_{1} \\ \vdots \\ O_{k} \\ \mathbb{0}_{n-k}
\end{bmatrix}}}
{\mathbb{P}\bracketit{\begin{bmatrix}\optimal{\tilde{z}}_{1}(d') \\ \vdots \\ \optimal{\tilde{z}}_{k}(d') \\ \mathbb{0}_{n-k}\end{bmatrix} + \begin{bmatrix}
\xi_{1} \\ \vdots \\ \xi_{k} \\ \mathbb{0}_{n-k}
\end{bmatrix} =
\begin{bmatrix}
O_{1} \\ \vdots \\ O_{k} \\ \mathbb{0}_{n-k}
\end{bmatrix}
\overset{(\text{ii})}{=}
\frac{
\mathbb{P}\bracketit{
\begin{bmatrix}
\xi_{1} \\ \vdots \\ \xi_{k}
\end{bmatrix}
\begin{bmatrix}
O_{1} \\ \vdots \\ O_{k}
\end{bmatrix}
\begin{bmatrix}\optimal{\tilde{z}}_{1}(d) \\ \vdots \\ \optimal{\tilde{z}}_{k}(d) \end{bmatrix}
\mathbb{P}\bracketit{
\begin{bmatrix}
\xi_{1} \\ \vdots \\ \xi_{k}
\end{bmatrix}
\begin{bmatrix}
O_{1} \\ \vdots \\ O_{k}
\end{bmatrix}
\begin{bmatrix}\optimal{\tilde{z}}_{1}(d') \\ \vdots \\ \optimal{\tilde{z}}_{k}(d') \end{bmatrix}
\\
\frac{\prod_{i=1}^{k}\text{exp}\braceit{-\frac{\varepsilon\norm{O_{i}-\optimal{\tilde{z}}_{i}(d)}_{1}}{\Delta_{\alpha}}}}
\prod_{i=1}^{k}\text{exp}\braceit{\frac{\varepsilon\norm{O_{i}-\optimal{\tilde{z}}_{i}(d')}_{1} - \varepsilon\norm{O_{i}-\optimal{\tilde{z}}_{i}(d)}_{1}}{\Delta_{\alpha}}}
\\
\prod_{i=1}^{k}\text{exp}\braceit{\frac{\varepsilon\norm{\optimal{\tilde{z}}_{i}(d) - \optimal{\tilde{z}}_{i}(d')}_{1}}{\Delta_{\alpha}}}
\text{exp}\braceit{\frac{\varepsilon\norm{\optimal{\tilde{z}}(d) - \optimal{\tilde{z}}(d')}_{1}}{\Delta_{\alpha}}}
\overset{(\text{v})}{\leqslant}\text{exp}\braceit{\varepsilon},
\end{align*}

where (i) is obtained from the primal feasibility condition $Z\in\mathcal{Q}(I)$, which enforces independence between the query random component and the sensitive data (see Definition <ref>), (ii) comes from rearranging the terms and removing zero entries, (iii) is due to the definition of the probability density function of the Laplace distribution, (iv) follows the reverse inequality of norms, and (v) is from the definition of $\ell_{1}-$sensitivity on $\alpha-$indistinguishable input datasets.

§ PROOF OF THEOREM <REF>

Without loss of generality, consider that the sum query $Sz(\xi)=S(\tilde{z} + Z\xi)$ requires releasing $p$ amount of sum statistics over non-intersecting subsets of $z(\xi)$. We thus need to show that the ratio of probabilities that the algorithm returns the same outcome $O\in\mathbb{R}^{p}$, i.e.,

\begin{align*}
\mathbb{P}\bracketit{S(\optimal{\tilde{z}}(d)+\optimal{Z}(d))\xi = O}/\;\mathbb{P}\bracketit{S(\optimal{\tilde{z}}(d')+\optimal{Z}(d'))\xi = O}\leqslant\text{exp}(\varepsilon),
\end{align*}

is bounded by a constant $\text{exp}(\varepsilon)$, where $S\in\mathbb{R}^{p\times n}$ and $\xi\in\mathbb{R}^{p}$ as in Defintion <ref>. By denoting the optimal solution of the chance-constrained problem (<ref>)  by $\optimal{\tilde{z}}$ and $\optimal{Z}\;$, this ratio writes as

\begin{align*}
&\frac{\mathbb{P}\bracketit{S\optimal{\tilde{z}}(d)+S\optimal{Z}(d)\xi = O}}{\mathbb{P}\bracketit{S\optimal{\tilde{z}}(d')+S\optimal{Z}(d')\xi = O}}
\overset{(\text{i})}{=}
\frac{\mathbb{P}\bracketit{S\optimal{\tilde{z}}(d)+\xi = O}}{\mathbb{P}\bracketit{S\optimal{\tilde{z}}(d')+\xi = O}}
\frac{\mathbb{P}\bracketit{\xi = O - S\optimal{\tilde{z}}(d)}}{\mathbb{P}\bracketit{\xi = O - S\optimal{\tilde{z}}(d')}}
\\
\frac{\prod_{i=1}^{p}\text{exp}\braceit{-\frac{\varepsilon\norm{O_{i} - [S\optimal{\tilde{z}}\;(d)]_{i}}_{1}}{\Delta_{\alpha}}}}
{\prod_{i=1}^{p}\text{exp}\braceit{-\frac{\varepsilon\norm{O_{i} - [S\optimal{\tilde{z}}\;(d')]_{i}}_{1}}{\Delta_{\alpha}}}}
\prod_{i=1}^{p}\text{exp}\braceit{\frac{\varepsilon\norm{O_{i} - [S\optimal{\tilde{z}}\;(d')]_{i}}_{1}- \varepsilon\norm{O_{i} - [S\optimal{\tilde{z}}\;(d)]_{i}}_{1}}{\Delta_{\alpha}}}
\\
\prod_{i=1}^{p}\text{exp}\braceit{\frac{\varepsilon\norm{[S\optimal{\tilde{z}}\;(d)]_{i}-[S\optimal{\tilde{z}}\;(d')]_{i}}_{1}}{\Delta_{\alpha}}}
\text{exp}\braceit{\frac{\varepsilon\norm{S\optimal{\tilde{z}}\;(d)-S\optimal{\tilde{z}}\;(d')}_{1}}{\Delta_{\alpha}}}
\overset{(\text{iv})}{\leqslant}
\text{exp}\braceit{\varepsilon},
\end{align*}

where (i) follows from the primal feasibility condition $Z\in\mathcal{Q}(S)$, which requires the random component of the sum query to be independent from the data (see Definition <ref>), (ii) is due to the definition of the probability density function of the Laplace distribution, (iii) follows from the reserve inequality of norms, and (iv) is from the $\ell_{1}-$sensitivity of the sum query, which is identical to the $\ell_{1}-$sensitivity of the identity query.

§ PROOF OF THEOREM <REF>

Similarly to the proofs of Theorems 2 and 3, the random components of the identity and linear queries can be shown to be independent from a datasets $d$ and $d'$ using the query specific feasibility conditions $\mathcal{Q}$. The reminder of the proof can be obtained by following the same steps of the proof in <cit.>, using notation $f(d)=I\optimal{\tilde{z}}(d)$ for the identity query and $f(d)=S\optimal{\tilde{z}}(d)$ for the sum query, where $f(\cdot)$ is the function of interest in <cit.>.

§ PROOF OF THEOREM <REF>

Input: $d,\Delta_\alpha, \varepsilon, \eta, \mu, I$

$(\optimal{\tilde{z}}\;, \optimal{Z}\;) \gets $
Solve (<ref>) for $Z\in\mathcal{Q}(I)$
and $\xi \sim \text{Lap}(T\Delta_\alpha/\varepsilon)$

$i=1,\ldots,T = \lceil\frac{\log(\mu)}{\log(\eta)}\rceil$

$\hat{\xi} \gets$ Sample from $\text{Lap}

Compute solution to (<ref>) as $\hat{z}=\optimal{\tilde{z}}\;+\optimal{Z}\hat{\xi}$

$A\hat{z}\leqslant b\; \lor\; i=T$
Release: $I\hat{z}$

Private identity query (PIQ)

Input: $d,\Delta_\alpha, \varepsilon, \eta,  \mu, S$

$(\optimal{\tilde{z}}\;, \optimal{Z}\;) \gets $
Solve (<ref>) for $Z\in\mathcal{Q}(S)$
and $\xi \sim \text{Lap}(T\Delta_\alpha \varepsilon)$

$i=1,\ldots,T = \lceil\frac{\log(\mu)}{\log(\eta)}\rceil$

$\hat{\xi} \gets$ Sample from $\text{Lap}

Compute solution to (<ref>) as $\hat{z}=\optimal{\tilde{z}}\;+\optimal{Z}\hat{\xi}$

$A\hat{z}\leqslant b\; \lor\; i = T$
Release: $S\hat{z}$

Private sum query (PSQ)

The iterative versions of the Algorithms <ref> and <ref> are provided, respectively, in Algorithms <ref> and <ref>.

Consider the optimal solution $(\optimal{\tilde{z}}{}, \optimal{Z}\;)$ returned by the chance constraint problem (line 2) and recall that  the sampling process (lines 4–5) generates a $(1-\eta)$-feasible solution $\hat{z}$.

The new algorithms, illustrated in Algorighms <ref> and <ref>, alternate this step with a constraint satisfaction test (line 6) for a maximum number of number $T$ of times with the goal of generating a solution that satisfies the problem constraints with probability at least $1-\mu$. Recall that the constraint satisfaction test, performed in line 6, can be achieved at no extra privacy loss (see Section <ref> for details).

The repetition of such process can be seen as a sequence of independent Bernulli trials, each with probability $1-\eta$ of success (i.e., $\hat{z}$ satisfies the problem constraints) and probability $\eta$ of failure (i.e., $\hat{z}$ violates the problem constraints).
Let $\text{FAIL}$ be the discrete random variable describing the number of unsuccessful trials prior to the first success. Thus, $\text{FAIL}$ is described by a Geometric random variable with probability $(1-\eta)$.
Formally, the goal is described by the following problem:

\[
T\triangleq\arg\min_T \mathbb{P}(\text{FAIL} \geqslant T) \leqslant \mu,
\]

requiring that the first success is seen after $T$ trials with probability no larger than $\mu$. Using a Geometric distribution of order $T$, it follows that:

\begin{align*}
\mathbb{P}(\text{FAIL} \geqslant T) &= 1 - \mathbb{P}(\text{FAIL} < T) \\
&= 1 - (1-\eta) \sum_{i=1}^{T-1} \eta^i \\
&= 1 - (1-\eta) \frac{1-\eta^T}{1-\eta} = \eta^T
\end{align*}

Thus, the solution to the minimizer above is for
\(
T = \left\lceil \frac{\log(\mu)}{\log(\eta)} \right\rceil.
\)