Information-Theoretic Privacy through Chaos Synchronization and Optimal Additive Noise

TL;DR

This paper introduces a privacy-preserving method using synchronized chaotic oscillators to generate optimal additive noise, minimizing information leakage in data queries over public channels.

Contribution

It proposes a novel approach combining chaos synchronization with convex optimization to enhance data privacy in communication systems.

Findings

Optimal noise distribution reduces mutual information effectively.

Chaotic oscillators can be synchronized to generate identical noise realizations.

Simulations demonstrate the method's effectiveness in privacy preservation.

Abstract

We study the problem of maximizing privacy of data sets by adding random vectors generated via synchronized chaotic oscillators. In particular, we consider the setup where information about data sets, queries, is sent through public (unsecured) communication channels to a remote station. To hide private features (specific entries) within the data set, we corrupt the response to queries by adding random vectors. We send the distorted query (the sum of the requested query and the random vector) through the public channel. The distribution of the additive random vector is designed to minimize the mutual information (our privacy metric) between private entries of the data set and the distorted query. We cast the synthesis of this distribution as a convex program in the probabilities of the additive random vector. Once we have the optimal distribution, we propose an algorithm to generate pseudo-random realizations from this distribution using trajectories of a chaotic oscillator. At the other end of the channel, we have a second chaotic oscillator, which we use to generate realizations from the same distribution. Note that if we obtain the same realizations on both sides of the channel, we can simply subtract the realization from the distorted query to recover the requested query. To generate equal realizations, we need the two chaotic oscillators to be synchronized, i.e., we need them to generate exactly the same trajectories on both sides of the channel synchronously in time. We force the two chaotic oscillators into exponential synchronization using a driving signal. Simulations are presented to illustrate our results.

Tables

Table 1. Table 1: Probability mass functions of X 𝑋 X and Y 𝑌 Y , and part of the one of ( X , Y ) 𝑋 𝑌 (X,Y) .

$X$	$\left[\begin{array}{c}0\\ 0\end{array}\right]$	$\left[\begin{array}{c}0\\ 1\end{array}\right]$	$\left[\begin{array}{c}0\\ 2\end{array}\right]$	$\left[\begin{array}{c}0\\ 3\end{array}\right]$	$\left[\begin{array}{c}0\\ 4\end{array}\right]$	$\left[\begin{array}{c}1\\ 0\end{array}\right]$	$\left[\begin{array}{c}1\\ 1\end{array}\right]$	$\left[\begin{array}{c}1\\ 2\end{array}\right]$	$\left[\begin{array}{c}1\\ 3\end{array}\right]$	$\left[\begin{array}{c}1\\ 4\end{array}\right]$
$p_X\left(x\right)$	0.5888	0.0200	0.0056	0.0560	0.0038	0.2616	0.0110	0.0042	0.0468	0.0022

Table 2. Table 2: Optimal distribution p V ∗ ​ ( v ) superscript subscript 𝑝 𝑉 𝑣 p_{V}^{*}(v) of the distorting additive random variable V 𝑉 V .

$V$	1	2	3	4	5	6	7	8	9
$p_V^{\ast }\left(v\right)$	0.1664	0.1522	0.1518	0.1355	0.1033	0.0832	0.0690	0.0591	0.0795

Equations

$$I[X;Y]:=\sum _{x\in \mathcal{X}}\sum _{y\in \mathcal{Y}}{p}_{X,Y}(x,y)\log \frac{p_{X,Y}(x,y)}{p_X(x)p_Y(y)}.$$

$$\dot{x}(t)=r(x(t),u(t)),$$

$$Q(x,u)=\frac{1}{2}\left(P\left(\frac{\partial r}{\partial x}(x,u)\right)+{\left(\frac{\partial r}{\partial x}(x,u)\right)}^{T}P\right),$$

$$\frac{d}{dt}\Big{(}\big{(}x_{1}(t)-x_{2}(t)\big{)}^{\top}P\big{(}x_{1}(t)-x_{2}(t)\big{)}\Big{)}\leq-\alpha\left|x_{1}(t)-x_{2}(t)\right|^{2},\hskip 2.84526ptt\in{\mathds{R}}_{\geq 0},$$

$$\left\{\begin{aligned} &p_{V}^{*}(v):=\operatorname*{arg\,min}_{p_{V}(v)}\ I[X;V+Y],\\ &\hskip 5.69054pt\text{\emph{s.t. }}V\rotatebox[origin={c}]{90.0}{$\models$}Y\text{ \emph{and} }p_{V}(v)\in\text{ \emph{Simplex}}.\\ \end{aligned}\right.$$

$$\{\begin{array}{rl}\dot{{\zeta }_1}(t)& =r({\zeta }_1(t),u(t)),\\ s_1(t)& =h({\zeta }_1(t)),\end{array}$$

$$\{\begin{array}{rl}\dot{\xi }(t)& =d(\xi (t)),\\ u(t)& =l(\xi (t)),\end{array}$$

$$\{\begin{array}{rl}\dot{{\zeta }_2}(t)& =r({\zeta }_2(t),u(t)),\\ s_2(t)& =h({\zeta }_2(t)),\end{array}$$

$$I[X;Z]$$

$$p_{Z|X}(z|x)$$

$$p_Z(z)$$

$$\text{Pr}[Z=z,Y=y]$$

$$\stackrel{\text{(a})}{=}\text{Pr}[V=z-y]\text{Pr}[Y=y]=p_V(z-y)p_Y(y),$$

$$p_{Z|Y}(z|y)$$

$$=\frac{p_V(z-y)p_Y(y)}{p_Y(y)}=p_V(z-y).$$

$$p_Z(z)$$

$$=\sum _{y\in \mathcal{Y}}\text{Pr}[V=z-y,Y=y]$$

$$\stackrel{(b)}{=}\sum _{y\in \mathcal{Y}}\text{Pr}[Y=y]\text{Pr}[V=z-y]=\sum _{y\in \mathcal{Y}}{p}_Y(y)p_V(z-y),$$

$$\{\begin{array}{rl}& p_V^{\ast }(v)=\underset{p_V(v)}{\mathrm{argmin}}\sum _{x\in \mathcal{X}}\sum _{y\in \mathcal{Y}}\sum _{z\in \mathcal{Z}}{p}_X(x)p_{Y|X}(y|x)p_V(z-y)\log \frac{{\sum }_{y\in \mathcal{Y}}{p}_{Y|X}(y|x)p_V(z-y)}{{\sum }_{y\in \mathcal{Y}}{p}_Y(y)p_V(z-y)},\\ & \text{ s.t. }{p}_V(v)\in \text{ Simplex}.\end{array}$$

$$\frac{1}{2}\left(P\left(\frac{\partial r}{\partial \zeta }(\zeta ,u)\right)+{\left(\frac{\partial r}{\partial \zeta }(\zeta ,u)\right)}^{T}P\right),$$

$$v_{k}=\psi(s_{k}):=\small\left\{\begin{array}[]{l}y_{1}$ \hskip 4.97922ptif $s_{k}\in c^{1},\\ \hskip 31.29802pt\vdots\\ y_{M}$ if $s_{k}\in c^{M}.\end{array}\right.$$

$$E[d(Y,\hat{Y})]$$

$$=\sum _{(y,\hat{y})\in {\mathcal{V}}_{\delta }}{p}_Y(y)p_{\hat{Y}|Y}(\hat{y}|y)|y-\hat{y}{|}^2\le \sum _{(y,\hat{y})\in {\mathcal{V}}_{\delta }}{p}_Y(y)|y-\hat{y}{|}^2=:{\bar{d}}_{\delta },$$

$$\{\begin{array}{rl}\dot{{\xi }_1}(t)& =10({\xi }_2(t)-{\xi }_1(t)),\\ \dot{{\xi }_2}(t)& =28{\xi }_1(t)-{\xi }_2(t)-{\xi }_1(t){\xi }_3(t),\\ \dot{{\xi }_3}(t)& =-\frac{8}{3}{\xi }_3(t)+{\xi }_1(t){\xi }_2(t),\\ u(t)& ={\xi }_1(t),\end{array}$$

$$\displaystyle C=\big{\{}$$

$$\displaystyle[2.3321,3.4341),[3.4341,4.5985),[4.5985,5.7743),[5.7743,\infty]\big{\}}.$$

Full text

**Information -Theoretic Privacy through Chaos Synchronization and

Optimal Additive Noise

**Carlos Murguia1,a, Iman Shames1,b, Farhad Farokhi1,2,c, and Dragan Nešić1,d

1Department of Electrical and Electronic Engineering, University of Melbourne, Australia

2The Commonwealth Scientific and Industrial Research Organisation (CSIRO), Data61, Australia

acarlos.murguia@unimelb.edu.au; biman.shames@unimelb.edu.au; cfarhad.farokhi@unimelb.edu.au;

ddnesic@unimelb.edu.au

1 Abstract

We study the problem of maximizing privacy of data sets by adding random vectors generated via synchronized chaotic oscillators. In particular, we consider the setup where information about data sets, queries, is sent through public (unsecured) communication channels to a remote station. To hide private features (specific entries) within the data set, we corrupt the response to queries by adding random vectors. We send the distorted query (the sum of the requested query and the random vector) through the public channel. The distribution of the additive random vector is designed to minimize the mutual information (our privacy metric) between private entries of the data set and the distorted query. We cast the synthesis of this distribution as a convex program in the probabilities of the additive random vector. Once we have the optimal distribution, we propose an algorithm to generate pseudorandom realizations from this distribution using trajectories of a chaotic oscillator. At the other end of the channel, we have a second chaotic oscillator, which we use to generate realizations from the same distribution. Note that if we obtain the same realizations on both sides of the channel, we can simply subtract the realization from the distorted query to recover the requested query. To generate equal realizations, we need the two chaotic oscillators to be synchronized, i.e., we need them to generate exactly the same trajectories on both sides of the channel synchronously in time. We force the two chaotic oscillators into exponential synchronization using a driving signal. Exponential synchronization implies that trajectories of the oscillators converge to each other exponentially fast for all admissible initial conditions and are perfectly synchronized in the limit only. Thus, in finite time, there is always a “small” difference between their trajectories. To implement our algorithm, we assume (as it is often done in related work) that systems have been operating for sufficiently long time so that this small difference is negligible and oscillators are practically synchronized. We quantify the worst-case distortion induced by assuming perfect synchronization, and show that this distortion vanishes exponentially fast. Simulations are presented to illustrate our results.

Keywords: Privacy; Data Sets, Queries, Mutual Information, Chaos.

2 Introduction

In a hyperconnected world, scientific and technological advances have led to an overwhelming amount of user data being collected and processed by hundreds of companies over public networks. Companies mine this data to provide targeted advertising and personalized services. However, these new technologies have also led to an alarming widespread loss of privacy in society. Depending on adversary’s resources, opponents may infer private user information from public data available on the internet and unsecured/public servers. A motivating example of privacy loss is the potential use of data from smart electrical meters by criminals, advertising agencies, and governments, for monitoring the presence and activities of occupants [1, 2]. Other examples are privacy loss caused by information sharing in distributed control systems and cloud computing [3]; the use of travel data for traffic estimation in intelligent transportation systems [4]; and data collection and sharing by the Internet-of-Things (IoT) [5], which is, most of the time, done without the user’s informed consent. These privacy concerns show that there is an acute need for privacy preserving mechanisms capable of handling the new privacy challenges induced by an interconnected world.

In this manuscript, we consider the problem of hiding private information $X$ of users (modeled as discrete random vectors) within datasets when publicly sharing requested queries $Y(X)$ from the same source. In particular, the aim of our privacy scheme is to respond to queries with distorted queries of the form $Z=Y(X)+V$ such that, when releasing $Z$, the private $X$ is “hidden”. Realizations of the vector $Z$ are transmitted over a public (unsecured) communication channel to a remote station. Then, if we do not distort $Y(X)$ before transmission, information about $X$ is directly accessible through the public channel. The first problem that we address is the design of the probability distribution of $V$ to maximize privacy, i.e., the distribution of $V$ must be constructed so that $Z=Y(X)+V$ carries as little information about $X$ as possible. Here, we follow an information-theoretic approach to privacy. We use the mutual information between private information $X$ and distorted queries $Y(X)+V$, $I[X;Y(X)+V]$, as privacy metric. The design of the discrete additive vector is casted as an optimization problem where we minimize $I[X;Y(X)+V]$ using the probability mass function of $V$, $p_{V}(v)$, as optimization variables. That is, the optimal distribution, $p^{*}_{V}(v)$, is given by $p^{*}_{V}(v):=\operatorname*{arg\,min}_{p_{V}(v)}I[X;Y(X)+V]$, where $p_{V}(v)$ is taken over a class of probability mass functions. Contrary to related work [6]-[11], we do not consider any sort of privacy-distortion trade-off in our formulation. We actually aim at making $I[X;Y(X)+V]$ as small as possible regardless of the distortion between $Y(X)$ and $Y(X)+V$ induced by $V$. Distortion is not an issue because we seek to generate exactly the same realization of $V$ at the remote station; then, we could recover the query by simply subtracting this realization from the one of $Z=Y(X)+V$. In order to accomplish this, we propose an algorithm to generate pseudorandom realizations from $p^{*}_{V}(v)$ at both sides of the channel using trajectories of two synchronized chaotic oscillators.

There are a number of requirements that the oscillators must satisfy for our algorithm to work: 1) trajectories of the oscillators must be bounded and chaotic; 2) they must be synchronized, i.e., we need them to generate exactly the same trajectories on both sides of the channel synchronously in time; and 3) the synchronous solution, regarded as a random process, must be stationary. Before giving the algorithm, we provide general guidelines for selecting the dynamics of the oscillators so that all the aforementioned requirements are satisfied. In particular, we use a range of well-known results in the literature to provide a synthesis procedure that allows to choose suitable oscillators. For boundedness, we use the notion of Input-to-State-Stability (ISS); for chaos, we employ standard largest Lyapunov exponent methods [12] and the (0-1) test [13]; for synchronization, we introduce the notion of convergent systems [14]; and for stationarity, we use hyperbolicity of the chaotic trajectories [15].

To generate equal realizations, our algorithm needs trajectories of the two chaotic oscillators (one at each side of the channel) to be synchronized. We force the oscillators into exponential synchronization using a driving signal. Exponential synchronization implies that trajectories of the oscillators converge to each other exponentially for all admissible initial conditions and are perfectly synchronized in the limit only. Therefore, in finite time, there is always a “small” difference between their trajectories. However, because oscillators synchronize exponentially fast, and it is often possible in practice to select initial conditions from a known compact set (known to both sides of the channel), it is safe to assume that the interconnected systems have been operating for sufficiently large time such that oscillators are practically synchronized, i.e., the synchronization error is so small that trajectories can be assumed to be equal. This is a standard assumption that is made in most, if not all, of the existing work on chaotic encryption based on synchronization [16]-[20]. Here, we give sufficient conditions for exponential synchronization to occur, provide tools for selecting the oscillators such that these conditions are satisfied, and assume that, after transients have settled down, trajectories are perfectly synchronized to some chaotic trajectory, say $\phi(t)\in{\mathds{R}}^{n_{\zeta}}$, $\zeta\in{\mathds{N}}$. If $n_{\zeta}>1$, our algorithm uses any entry $\phi^{s}(t)\in\mathcal{S}\subset{\mathds{R}}$ of $\phi(t)$ to generate realizations from $p^{*}_{V}(v)$, where $\mathcal{S}$ denotes some compact set that characterizes the support of $\phi^{s}(t)$. Because oscillators are selected such that $\phi(t)$, regarded as a random process, is stationary, samples from $\phi^{s}(t)$ follow a stationary probability density function. We obtain this density through Monte Carlo simulations [21] and divide its support $\mathcal{S}$ into a finite set of cells $C=\{c^{1},\ldots,c^{M}\}$ such that the probability that $\phi^{s}(t)$ lies in these cells equals the optimal probability distribution $p_{V}^{*}(v)$. That is, we generate pseudorandom realizations from $p_{V}^{*}(v)$ by properly selecting $C$ and evaluating if $\phi^{s}(t)$ lies in $C$ at the sampling instants.

The use of additive noise to preserve privacy is common practice. There are mainly two classes of privacy metrics considered in the literature; namely, differential privacy [22]-[23] and information-theoretic metrics, e.g., mutual information, conditional entropy, Kullback-Leibler divergence, and Fisher information [24]-[28]. In differential privacy, because it provides certain privacy guarantees, Laplace noise is usually used [29]. However, when maximal privacy is desired, Laplace noise is generally not the optimal solution. This raises the fundamental question: what is the noise distribution achieving maximal privacy? This question has many possible answers depending on the particular privacy metric being considered and the system configuration, see, e.g., [6]-[8],[11], for differential privacy based results, and [24]-[28], for information theoretic results. In general, if the data to be kept private follows continuous distributions, the problem of finding the optimal additive noise to maximize privacy is hard to solve. If a close-form solution for the distribution is desired, the problem amounts to solving a set of nonlinear partial differential equations which, in general, might not have a solution, and even if they do have a solution, it is hard to find [24]. This problem has been addressed by imposing some particular structure on the considered distributions or assuming the data to be kept private is deterministic [24],[7],[8]. The authors in [7],[8] consider deterministic input data sets and treat optimal distributions as distributions that concentrate probability around zero as much as possible while ensuring differential privacy. Under this framework, they obtain a family of piecewise constant probability density functions that achieve minimal distortion for a given level of privacy. In [24], the authors consider the problem of preserving the privacy of deterministic databases using additive continuous noise with constrained support. They use the Fisher information and the Cramer-Rao bound to construct a privacy metric between deterministic data and the one with the additive noise, and find the probability density function that minimizes it. Moreover, they prove that, in the unconstrained support case, the optimal noise distribution minimizing the Fisher information is Gaussian. This observation has been also made in [30] when using mutual information as a measure of privacy. We remark that most of the aforementioned papers consider privacy-distortion trade-offs when designing their distorting mechanisms. We do not consider this trade-off here because, at the end of the channel, we remove the distortion that we induce using our synchronization based formulation.

Existing work on chaotic encryption based on synchronization [16]-[20] directly uses the states of the chaotic oscillators to mask private information. That is, standard algorithms do not use chaotic trajectories to generate pseudorandom realization from probability distributions (as we do here); instead, they simply add the value of the sampled chaotic trajectory (or functions of it) to private messages. Although the latter succeeds in masking messages, it does not give any privacy guarantees (neither information-theoretic nor in a differential privacy sense) on the private information, and it is not optimal in any sense. Hence, the contributions of our scheme with respect to existing work on chaotic encryption [16]-[20] are the treatment of fully stochastic datasets, the information-theoretic privacy guarantees that our framework provides, and the optimal performance of the designed distorting additive vector (optimal in the sense of minimizing the mutual information $I[X;Y(X)+V]$). The work here is inspired by the experimental results presented in [31], where the authors propose a framework similar to ours for deterministic data using a electronic circuit implementation of the Mackey-Glass chaotic oscillator [32]. The contribution of our work with respect to [31] is threefold: 1) we consider fully stochastic data, which makes the privacy scheme fundamentally very different; 2) we provide a general formulation that encompasses a large class of chaotic systems, not only the electronic circuit implementation of the Mackey-Glass oscillator; and 3) we generate realizations from optimal distorting distributions, in [31], they consider uniform distributions only which is not optimal for stochastic data.

Next, we summarize the main contributions of the chapter.

**Contributions:

**

1. We provide a general information-theoretic privacy framework based on optimal additive distorting random vectors and synchronization of chaotic oscillators; 2) We prove that the synthesis of the probability mass function $p_{V}(v)$ of the distorting random vector $V$ can be posed as a convex program in $p_{V}(v)$ over a class of probability mass functions; 3) We provide an algorithm to generate pseudorandom realizations from this distribution using trajectories of chaotic oscillators; 4) Using off-the-shelf results in the literature, we provide a synthesis procedure for selecting the dynamics of the oscillators so that our algorithm is guaranteed to work.

The remainder of the paper is organized as follows. In Section 3, we present some preliminaries results needed for the subsequent sections. We introduce the notion of convergent systems and the concept of mutual information. The general formulation and the specific problems to be addressed are given in Section 4. In Section 5, we pose the synthesis of the probability distribution of the optimal distorting vector. General guidelines for selecting the chaotic oscillators are given in Section 6. The algorithm for generating pseudorandom realizations from the optimal distribution is presented in Section 7. Simulation results are given in Section 8 and concluding remarks in Section 9.

3 Notation and Preliminaries

The symbol ${\mathds{R}}$ stands for the real numbers, ${\mathds{R}}_{>0}$(${\mathds{R}}_{\geq 0}$) denotes the set of positive (non-negative) real numbers. The symbol ${\mathds{N}}$ stands for the set of natural numbers. The Euclidian norm in ${\mathds{R}}^{n}$ is denoted simply as $|\cdot|$, $|x|^{2}=x^{\top}x$, where ⊤ defines transposition. For a given measurable function $u(t)$, $t\in{\mathds{R}}_{\geq 0}$, we denote its $\mathcal{L}_{\infty}$ norm as $||u||_{\infty}:=\operatorname*{ess\,sup}_{t\geq 0}|u(t)|$, where $\operatorname*{ess\,sup}$ denotes essential supremum. Matrices composed of only ones and only zeros of dimension $n\times m$ are denoted by $\mathbf{1}_{n\times m}$ and $\mathbf{0}_{n\times m}$, respectively, or simply $\mathbf{1}$ and $\mathbf{0}$ when their dimensions are clear. For square matrices $A\in{\mathds{R}}^{n\times n}$, $\rho[A]$ denotes the spectral radius of $A$. A continuous function $\gamma:[0,a)\rightarrow[0,\infty)$ is said to belong to class $\mathcal{K}$ if it strictly increasing and $\gamma(0)=0$. Similarly, a continuous function $\beta:[0,a)\times[0,\infty)\rightarrow[0,\infty)$ belongs to class $\mathcal{KL}$ if, for fixed $s$, $\beta(r,s)$ belongs to class $\mathcal{K}$ with respect to $r$ and, for fixed $r$, $\beta(r,s)$ is decreasing with respect to $s$ and $\lim_{s\rightarrow\infty}\beta(r,s)=0$. Consider a discrete random vector $X$ with alphabet $\mathcal{X}=\{x_{1},\ldots,x_{N}\}$, $x_{i}\in{\mathds{R}}^{m}$, $m\in{\mathds{N}}$, $i\in\{1,\ldots,N\}$, and probability mass function $p_{X}(x)=\text{Pr}[X=x]$, $x\in\mathcal{X}$, where $\text{Pr}[B]$ denotes probability of event $B$. Similarly, for two random vectors $X$ and $Y$, taking values in the alphabets $\mathcal{X}$ and $\mathcal{Y}$, respectively, their joint probability mass function is denoted by $p_{X,Y}(x,y)$, the marginal distribution of $X$ is given by $p_{X}(x)=\sum_{y\in\mathcal{Y}}p_{X,Y}(x,y)$, and the conditional distribution of $X$ given $Y$ as $p_{Y|X}(y|x)=p_{X,Y}(x,y)/p_{X}(x)$. Analogously, for a continuous random vector $Y$, we denote their (multivariate) probability density function as $f_{Y}(y)$. The notation $X\sim f_{X}(x)$ ($X\sim p_{X}(x)$) stands for continuous (discrete) random vectors $X$ following the probability density (mass) function $f_{X}(x)$ ($p_{X}(x)$). We denote by "Simplex" the probability simplex defined by $\sum_{x\in\mathcal{X}}p_{X}(x)=1$, $p_{X}(x)\geq 0$ for all $x\in\mathcal{X}$. The notation $E[a]$ denotes the expected value of the random vector $a$. We denote independence between two random vectors, $X$ and $Y$, as X\rotatebox[origin={c}]{90.0}{\models}Y.

3.1 Mutual Information

Definition 1**.**

Consider two random vectors, $X$ and $Y$, with joint probability mass function $p_{X,Y}(x,y)$ and marginal probability mass functions, $p_{X}(x)$ and $p_{Y}(y)$, respectively. Their mutual information $I[X;Y]$ is defined as the relative entropy between the joint distribution and the product distribution $p_{X}(x)p_{Y}(y)$, i.e.,

[TABLE]

Mutual information $I[X;Y]$ between two jointly distributed vectors, $X$ and $Y$, is a measure of the dependence between $X$ and $Y$.

3.2 Convergent Systems

Consider the dynamical system:

[TABLE]

with $t\in{\mathds{R}}_{\geq 0}$, state $x\in{\mathds{R}}^{n}$, input $u\in\mathcal{U}\subseteq{\mathds{R}}^{m}$, and vector field $r:{\mathds{R}}^{n}\times\mathcal{U}\rightarrow{\mathds{R}}^{n}$. The vector field $r(x,u)$ is continuously differentiable in $x$, and $u(t)$ is piecewise continuous in $t$ and takes values in some compact set $\mathcal{U}\subseteq{\mathds{R}}^{m}$.

Definition 2**.**

*[33]**. System (1) is said to be globally asymptotically convergent if and only if for any bounded input $u(t)$, $t\in{\mathds{R}}$, there is a unique bounded globally asymptotically stable solution $\bar{x}_{u}(t)$, $t\in{\mathds{R}}$, such that $\lim_{t\rightarrow\infty}\left|x(t)-\bar{x}_{u}(t)\right|=0$ for all initial conditions. *

For a convergent system, the limit solution is solely determined by the external excitation $u(t)$ and not by the initial conditions. A sufficient condition for convergence obtained by Demidovich [33] and later extended in [14] is presented in the following proposition.

Proposition 1**.**

[33, 14]**. If there exists a positive definite matrix $P\in{\mathds{R}}^{n\times n}$ such that all the eigenvalues $\lambda_{i}(Q)$ of the symmetric matrix

[TABLE]

are negative and separated from zero, i.e., there exists a constant $c\in{\mathds{R}}_{>0}$ such that $\lambda_{i}(Q)\leq-c<0,$ for all $i\in\{1,...,n\}$, $u\in\mathcal{U}$, and $x\in{\mathds{R}}^{n}$, then system (1) is globally exponentially convergent; and, for any pair of solutions $x_{1}(t),x_{2}(t)\in{\mathds{R}}^{n}$ of (1), the following is satisfied:

[TABLE]

*with constant $\alpha:=(c/\lambda_{\max}(P))$ and $\lambda_{\max}(P)$ being the largest eigenvalue of the symmetric matrix $P$. *

Remark 1**.**

*There are other methods to verify that trajectories of system (1) converge to a limit solution that is independent of the initial conditions and solely determined by the external excitation $u(t)$. For instance, contraction theory [34], Lyapunov function approach to incremental stability [35], the quadratic (QUAD) inequality approach (a Lipschitz-like condition) [36], and differential dissipativity [37], which are all concepts that are closely related to notion of convergent systems [14] that we use here. *

4 Problem Setup

Let $X$ be a discrete random vector that must be kept private. The alphabet and probability mass function of $X$ are denoted as $\mathcal{X}=\{x_{1},\ldots,x_{N}\}$, $x_{i}\in{\mathds{R}}^{n_{x}}$, $n_{x}\in{\mathds{N}}$, $i\in\{1,\ldots,N\}$ and $p_{X}(x)=\text{Pr}[X=x]$, $x\in\mathcal{X}$, respectively. The $n_{x}$ entries of $X$ represent, for instance, private entries of $n_{x}$ users within a dataset that is stored by a trusted server. The server admits queries of the form $Y=q(X)$, $Y\in{\mathds{R}}^{n_{y}}$, for some (stochastic or deterministic) mapping $q:{\mathds{R}}^{n_{x}}\rightarrow{\mathds{R}}^{n_{y}}$ characterized by the transition probabilities $p_{Y|X}(y|x)$, $x\in\mathcal{X}$, $y\in\mathcal{Y}$, where $\mathcal{Y}=\{y_{1},\ldots,y_{M}\}$, $y_{i}\in{\mathds{R}}^{n_{y}}$, $n_{y}\in{\mathds{N}}$. The aim of our privacy scheme is to respond to queries of the form $q(X)$ with distorted queries $Z=q(X)+V$, for some discrete random vector $V$ (with V\rotatebox[origin={c}]{90.0}{\models}Y), such that, when releasing $Z$, the individual entries of $X$ are “hidden”. Realizations of the vector $Z$ are transmitted over a public (unsecured) communication channel to a remote station, see Figure 1. Then, if we do not add $V$ to $q(X)$ before transmission, information about $X$ is directly accessible through the public channel. As a preliminary problem that we need to solve for the subsequent results, we address the design of the probability distribution of $V$ to maximize privacy, i.e., the distribution of $V$ must be constructed so that the sum, $Z=q(X)+V$, carries as little information about $X$ as possible. In this manuscript, we use the mutual information between $X$ and $Z=Y+V$, $I[X;Z]$, as privacy metric. We aim at finding the probability mass function of $V$, $p_{V}(v)$, that minimizes $I[X;Z]$ over a class of probability mass functions. That is, we cast the design of $p_{V}(v)$ as an optimization problem with cost function $I[X;Z]$, optimization variables $p_{V}(v)$, and subject to V\rotatebox[origin={c}]{90.0}{\models}Y and the usual probability simplex constraints. Note that, contrary to related work [9]-[11],[27],[28], we do not consider any sort of privacy-distortion trade-off in our formulation. We minimize $I[X;Y+V]$ regardless of the distortion between $Y$ and $Y+V$ induced by $V$. Distortion is not an issue because, we seek to generate exactly the same realization of $V$ at the remote station and then recover the query by subtracting this realization from the one of $Z=Y+V$. This is addressed in Problem 2 and Problem 3 below.

We let $V$ be a discrete random vector with alphabet $\mathcal{Y}$ and probability mass function $p_{V}(v)=\text{Pr}[V=v]$, $v\in\mathcal{Y}$, i.e., the alphabet of $V$ and the one of the query $Y=q(X)$ are equal. Having equal alphabets imposes a tractable convex structure on the cost $I[X;Z]$ and reduces the optimization variables to the probabilities of each element of the alphabet. The case with arbitrary alphabet leads to a combinatorial optimization problem where the objective changes its structure for different combinations. We do not address this case in this manuscript; it is left as a future work. In what follows, we formally present the optimization problem we seek to address.

Problem 1**.**

[Optimal Distribution of the Additive Distorting Signal]* For given $p_{X}(x)=\text{\emph{Pr}}[X=x]$ and $p_{Y|X}(y|x)=\text{\emph{Pr}}[Y=y|X=x]$, $x\in\mathcal{X}$, $y\in\mathcal{Y}$, find the probability mass function $p_{V}(v)=\text{\emph{Pr}}[V=v]$, $v\in\mathcal{Y}$ solution of the optimization problem:*

[TABLE]

Here, $p_{V}^{*}(v)$ denotes the optimal distribution solution to (3). To hide $X$, once we have obtained $p_{V}^{*}(v)$, we aim at generating realizations $v\in\mathcal{Y}$ from this distribution, add them to the required query ($Y=q(X)$), and send realizations of the sum $Z=Y+V$ to the remote station through the public channel. At the other end of the channel, we seek to generate the exact same realizations from $p_{V}^{*}(v)$ so that we can recover the query by simply subtracting $V$ from $Z$, see Figure 2. Note that, in Figure 2, we have a recovered $\hat{Y}$ at the remote station rather that the actual $Y$. This is because we want to remark that, due to practical errors in our algorithm–e.g., due to communication delays and transients–realizations of $V$ that we generate at both ends of the channel might be slightly different in practice. To generate these realizations, we use trajectories, $\phi^{\zeta}_{u,1}(t,\zeta_{1}(0),u(t))$, $t\in{\mathds{R}}_{\geq 0}$, $\zeta_{1}(0)\in{\mathds{R}}^{n_{\zeta}}$, $u(t)\in{\mathds{R}}^{n_{u}}$, of a chaotic dynamical system of the form:

[TABLE]

with state $\zeta_{1}(t)\in{\mathds{R}}^{n_{\zeta}}$, output $s_{1}(t)\in{\mathds{R}}$, continuous in $t$ input $u(t)\in\mathcal{U}\subset{\mathds{R}}^{n_{u}}$ taking values in some compact set $\mathcal{U}$, continuous function $h:{\mathds{R}}^{n_{\zeta}}\rightarrow{\mathds{R}}$, and vector field $r:{\mathds{R}}^{n_{\zeta}}\times\mathcal{U}\rightarrow{\mathds{R}}^{n_{\zeta}}$ continuously differentiable in its first argument, uniformly in its second argument. Hereafter, system (4) is referred to as responder 1. Responder 1 is placed at the side of the trusted server, see Figure 2. The input signal $u(t)$ is generated by a chaotic autonomous exosystem:

[TABLE]

with state $\xi(t)\in{\mathds{R}}^{n_{\xi}}$, output $u(t)\in\mathcal{U}\subset{\mathds{R}}^{n_{u}}$, and vector fields $d:{\mathds{R}}^{n_{\xi}}\rightarrow{\mathds{R}}^{n_{\xi}}$ and $l:{\mathds{R}}^{n_{\xi}}\rightarrow{\mathds{R}}^{n_{u}}$. The vector field $d(\xi)$ is locally Lipschitz in $\xi$ and $l(\xi)$ is continuous. We refer to (5) as the driver system. We let $u(t)$ be connected to the remote station via the public channel, see Figure 2. At the other end of the channel, driven by the same input signal $u(t)$, we have a third chaotic oscillator with the same dynamics as (4) but with potentially different initial conditions, i.e., the second oscillator is given by

[TABLE]

with state $\zeta_{2}(t)\in{\mathds{R}}^{n_{\zeta}}$ and output $s_{2}(t)\in{\mathds{R}}$. We denote trajectories of (6) as $\phi^{\zeta}_{u,2}(t,\zeta_{2}(0),u(t))$ with $t\in{\mathds{R}}_{\geq 0}$, $\zeta_{2}(0)\in{\mathds{R}}^{n_{\zeta}}$, and $u(t)\in\mathcal{U}\subset{\mathds{R}}^{n_{\zeta}}$. System (6) is referred to as responder 2. Note that if $\zeta_{1}(t)=\zeta_{2}(t)$, $t\in{\mathds{R}}_{\geq 0}$, i.e., if systems (4) and (6) are synchronized, and we use the synchronous chaotic solution, say $\phi^{\zeta}_{u}(t,u(t))$, to generate realizations from $p_{V}^{*}(v)$, we could have the same realization of $V$ at both sides of the channel.

Problem 2**.**

[Boundedness, Chaos, and Synchronization]* State sufficient conditions on the vector fields $r(\cdot)$, $h(\cdot)$, $d(\cdot)$, and $l(\cdot)$ of the coupled system (4)-(6) such that: 1) trajectories of (4)-(6) exist and are bounded and chaotic; and 2) systems (4) and (6) exponentially synchronize, i.e., $\lim_{t\rightarrow\infty}|\zeta_{1}(t)-\zeta_{2}(t)|=0$, exponentially fast. *

Remark 2**.**

*Problem 2 seeks to enforce exponential synchronization by selecting the dynamics of the oscillators. Exponential synchronization implies that trajectories of the responders converge to each other exponentially for all initial conditions and are perfectly synchronized in the limit only. It follows that, in finite time, there is always a “small” difference between their trajectories. Nevertheless, because oscillators synchronize exponentially fast, and it is often possible in practice to select initial conditions from a known compact set (known to both the trusted server and the remote station), it is safe to assume that the interconnected systems have been operating for sufficiently large time such that oscillators are practically synchronized, i.e., the synchronization error is so small that trajectories can be assumed to be equal. This is a standard assumption that is made in most, if not all, of the existing work on chaotic encryption based on synchronization [16]-[20]. *

Finally, once we have found functions solution to Problem 2, which guarantees exponential synchronization of the responders, and assuming that responders are synchronized (see Remark 2), we aim at designing a procedure to generate pseudorandom realizations from $p_{V}^{*}(v)$ using the synchronous chaotic solution $\phi^{\zeta}_{u}(t,u(t))$. Note that $\zeta_{1}(t)=\zeta_{2}(t)\Rightarrow s_{1}(t)=h(\zeta_{1}(t))=s_{2}(t)=h(\zeta_{2}(t))$, for all $t\geq 0$. Moreover, because $\zeta_{1}(t)=\zeta_{2}(t)=\phi^{\zeta}_{u}(t,u(t))$; then, $s_{1}(t)=s_{2}(t)=h(\phi^{\zeta}_{u}(t,u(t)))=:\phi^{s}_{u}(t,u(t))\in\mathcal{S}\subset{\mathds{R}}$ for some compact set $\mathcal{S}$.To reduce the complexity of the algorithm, we use the lower dimensional synchronous solution $\phi^{s}_{u}(t,u(t))$ to generate the realizations from $p_{V}^{*}(v)$.

Problem 3**.**

[Generation of Optimal Pseudorandom Numbers]* Using the lower dimensional synchronous solution, $\phi^{s}_{u}(t,u(t))$, design an algorithm to generate pseudorandom realizations from the optimal distribution $p_{V}^{*}(v)$, $v\in\mathcal{Y}$. *

5 Optimal Distribution of the Additive Distorting Signal

In this section, we prove that Problem 1 can be posed as a convex program in the probabilities $p_{V}(v)$, $v\in\mathcal{Y}$. We derive an explicit expression for the cost function $I[X;Z]$, $Z=Y+V$, in terms of the given $p_{X}(x)$ and $p_{Y|X}(y|x)$ and the variables $p_{V}(v)$, restricted to satisfy the independence constraint V\rotatebox[origin={c}]{90.0}{\models}Y.

Lemma 1**.**

$I[X;Z]$* with $Z=Y+V$, V\rotatebox[origin={c}]{90.0}{\models}Y, is a convex function of $p_{V}(v)$, $v\in\mathcal{Y}$, for given $p_{X}(x)$ and $p_{Y|X}(y|x)$, $x\in\mathcal{X}$, $y\in\mathcal{Y}$; and can be written compactly in terms of $p_{X}(x)$, $p_{Y|X}(y|x)$, and $p_{V}(v)$, as follows:*

[TABLE]

Proof: The expression on the right-hand side of (7a) follows by inspection of Definition 1 and the fact that $p_{Z,X}(z,x)=p_{X}(x)p_{Z|X}(z|x)$. By [38, Theorem 2.7.4], cost (7a) is convex in $p_{Z|X}(z|x)$ for given $p_{X}(x)$. However, our optimization variables are $p_{V}(v)$ and not $p_{Z|X}(z|x)$. Note that $X$, $Y$, and $Z$ form a Markov chain in that order [39]; therefore, $p_{X,Y,Z}(x,y,z)=p_{X}(x)p_{Y|X}(y|x)p_{Z|Y}(z|y)$. Marginalizing $p_{X,Y,Z}(x,y,z)$ with respect to $Y\in\mathcal{Y}$ and then conditioning with respect to $X$ yields $p_{X,Z}(x,z)=\sum_{y\in\mathcal{Y}}p_{X}(x)p_{Y|X}(y|x)p_{Z|Y}(z|y)$ and $p_{Z|X}(z|x)=\sum_{y\in\mathcal{Y}}p_{Y|X}(y|x)p_{Z|Y}(z|y)$, respectively. Note that $p_{Z|X}(z|x)$ is just a linear transformation of $p_{Z|Y}(z|y)$. Hence, convexity with respect to $p_{Z|X}(z|x)$ implies convexity with respect to $p_{Z|Y}(z|y)$ because convexity is preserved under affine transformations [40]. Next, consider $p_{Z|Y}(z|y)=p_{Z,Y}(z,y)/p_{Y}(y)$. By definition, $p_{Z,Y}(z,y)=\text{Pr}[Z=z,Y=y]$, $z\in\mathcal{Z}$, $y\in\mathcal{Y}$. Note that

[TABLE]

where (a) follows from independence between $V$ and $Y$. Thus,

[TABLE]

We have concluded convexity of $I[X;Z]$ with respect to $p_{Z|Y}(z|y)$ above. Hence, because $p_{Z|Y}(z|y)=p_{V}(z-y)$ and $p_{V}(z-y)$ is a linear transformation of $p_{V}(v)$ ($p_{V}(z-y)=p_{V}(v)$ for $z-y=v$ and zero otherwise), the cost $I[X;Z]$ is convex in $p_{V}(v)$. Moreover, since $p_{Z|X}(z|x)=\sum_{y\in\mathcal{Y}}p_{Y|X}(y|x)p_{Z|Y}(z|y)$ and $p_{Z|Y}(z|y)=p_{V}(z-y)$, equality (7b) holds true. It remains to prove that $p_{Z}(z)$ can be written as (7c). Because $Z=Y+V$, $p_{Z}(z)=\text{Pr}[Z=z]$, for a given $z\in\mathcal{Z}$, can be written as the sum of the probabilities of all $Y=y$ and $V=v$ that result in $z$, i.e.,

[TABLE]

*where (b) follows from independence between $V$ and $Y$. $\blacksquare$ *

By Lemma 1, the cost $I[X;Z]$, for V\rotatebox[origin={c}]{90.0}{\models}Y, is convex in $p(v)$ and parametrized by $p_{X}(x)$ and $p_{Y|X}(y|x)$. In what follows, we cast the nonlinear program for solving Problem 1.

Theorem 1**.**

Given $p_{X}(x)$ and $p_{Y|X}(y|x)$, $x\in\mathcal{X}$, $y\in\mathcal{Y}$, the mapping $p_{V}(v)$, $v\in\mathcal{Y}$, that minimizes $I[X;Z]$, $Z=V+Y$, subject to V\rotatebox[origin={c}]{90.0}{\models}Y can be found by solving the following convex program:

[TABLE]

Proof: Theorem 1 follows from Lemma 1.* *

6 Boundedness, Chaos, and Synchronization

6.1 Existence, Uniqueness, and Boundedness of Solutions

We start addressing existence, uniqueness, and boundedness of the solutions of the coupled systems (4)-(6). To be able to use synchronous solutions to generate realizations from $p_{V}^{*}(v)$, we first need these solutions to exist and be bounded. In the system description given above, we have assumed that $r(\zeta,u(t))$ is continuously differentiable in $\zeta$ uniformly in $u(t)$, $u(t)$ is continuous in $t$, and $d(\xi)$ is locally Lipschitz. These alone imply uniqueness and existence of solutions of (4)-(6) over some finite time interval $t\in[0,\tau]$, $\tau\in{\mathds{R}}_{>0}$, [41, Theorem 2.2].To conclude the latter for arbitrarily large $\tau$, besides the locally Lipschitz assumption on the functions, we need boundedness of the solutions of (4)-(6) [41, Theorem 2.4]. Note that the coupled systems (4)-(6) have a cascade structure. The driver dynamics is independent of the responders states, and its output, $u(t)$, is the input of the responders. Then, an approach to conclude boundedness of the overall system is to conclude boundedness of the driver first, and then boundedness of the responders when driven by bounded inputs. In what follows, we formally introduce the notion of boundedness that we use here.

Definition 3**.**

[41]* The solutions of (5) are bounded for a bounded set of initial conditions if there exists a positive constant $c$, independent of the initial time instant, and for every $a\in(0,c)$, there is $b=b(a)>0$, independent of the initial time instant, such that $|\xi(0)|\leq a\Rightarrow|\xi(t)|\leq b$, $\forall\hskip 2.84526ptt\geq 0$. If the latter holds for arbitrarily large $a$; then, the solutions of (5) are globally bounded. *

Remark 3**.**

*Because $l(\xi)$ is continuous, by the extreme value theorem, boundedness of $\xi(t)$ implies boundedness of $u(t)=l(\xi(t))$. *

Remark 4**.**

*We do not give conditions for boundedness of the solutions of (5). It is assumed that the vector field $d(\xi)$ is such that the solutions of the driver are globally bounded. We refer the reader to, for instance, [41, Theorem 4.18], where sufficient conditions for boundedness are given in terms of Lyapunov-like results. *

Next, for bounded solutions of the driver, we need the solutions of the responders to be bounded when driven by $u(t)$. To address this, we use the notion if Input-to-State-Stability (ISS) [42].

Definition 4**.**

[42]* System (4) (and thus system (5) as well) is said to be Input-to-State-Stable if there exist a class $\mathcal{KL}$ function $\beta(\cdot)$ and a class $\mathcal{K}$ function $\gamma(\cdot)$ such that for any initial condition $\zeta_{1}(0)$ and any bounded input $u(t)$, the solution $\zeta_{1}(t)$ exists for all $t\in{\mathds{R}}_{\geq 0}$ and satisfies: $|\zeta_{1}(t)|\leq\beta(\zeta_{1}(0),t)+\gamma\left(||u||_{\infty}\right)$. *

Remark 5**.**

ISS* of the responders with respect to $u(t)$ guarantees that, for any bounded $u(t)$, the states $\zeta_{1}(t)$ and $\zeta_{2}(t)$ are bounded. Moreover, as $t$ increases, $|\zeta_{1}(t)|$ and $|\zeta_{2}(t)|$ are ultimately bounded [41] by $\gamma\left(||u||_{\infty}\right)$, see [42] for further details. *

Remark 6**.**

*Sufficient conditions for the responders to be ISS with input $u(t)$ are not provided here. We assume that the vector field $r(\zeta,u(t))$ is such that systems (4) and (5) are ISS with respect to $u(t)$. We refer the reader to, for instance, [41, Theorem 4.19], where sufficient conditions for ISS are given in terms of ISS-Lyapunov functions. *

Remark 7**.**

*The weaker property of integral Input-to-State-Stable (iISS) [43] could be used to conclude boundedness of the responder’s trajectories when driven by “sufficiently small” inputs. We refer the reader to [44], where sufficient conditions for *iISS * and related boundedness results are given. *

6.2 Synchronization

Next, we give sufficient conditions on $r(\zeta,u(t))$ such that $\lim_{t\rightarrow\infty}|\zeta_{1}(t)-\zeta_{2}(t)|=0$, i.e., the responders exponentially synchronize. We assume that solutions of the coupled systems (4)-(6) exist and are bounded, i.e., vector fields $r(\cdot)$, $d(\cdot)$, and $l(\cdot)$ satisfy the conditions stated in the previous subsection. Then, for bounded $u(t)$, a sufficient condition for the responders to exponentially synchronize is that systems (4) and (6) are convergent systems in the sense of Definition 2. The latter implies that, because both responders are driven by the input $u(t)$ and their dynamics are described by the same set of differential equations, trajectories of (4) and (6) converge to the same the limit solution, $\phi^{\zeta}_{u}(t,u(t))$, and this solution is solely determined by $u(t)$ and not by the initial conditions. In the following corollary of Proposition 1, we give a sufficient condition for the responders to be exponentially convergent (and thus to exponentially synchronize).

Corollary 1**.**

Consider the responders (4) and (6). If there exists a positive definite matrix $P\in{\mathds{R}}^{n_{\zeta}\times n_{\zeta}}$ such that, for all $u\in{\mathds{R}}^{n_{u}}$ and $\zeta\in{\mathds{R}}^{n_{\zeta}}$, all the eigenvalues of the symmetric matrix:

[TABLE]

*are negative and separated from zero; then, responders (4) and (6) are globally exponentially convergent, and thus $\lim_{t\rightarrow\infty}|\zeta_{1}(t)-\zeta_{2}(t)|=0$, exponentially fast. *

Remark 8**.**

*If the driver’s output $u(t)$ is to be sent over a network and quantization (or some sort of coding) is required, we would need to drive responders by the same quantized $u(t)$, say $u_{Q}(t)$, to achieve exponential synchronization. That is, if we quantize $u(t)$ to obtain $u_{Q}(t)$, and we drive both responders by $u_{Q}(t)$ (with, e.g., a Zero-Order-Hold (ZOH)), they would also exponentially synchronize. They would synchronize to a different trajectory than when driven by $u(t)$, but they would synchronize exponentially fast. *

Besides the notion of convergent systems, there are other methods available in the literature that can be used to verify that trajectories of responders asymptotically synchronize to a limit solution that is independent of the initial conditions. See Remark 1 for details.

6.3 Chaotic Dynamics

There are mainly two branches of methods to identify chaotic dynamics; namely, standard largest Lyapunov exponent methods [12], and the more recent (0-1) test [13]. Both methods use trajectories (numerical or experimental) of the systems under study to decide whether they are chaotic or not. In general, there are no sufficient conditions directly on the differential equations (the vector fields $r(\cdot)$ and $d(\cdot)$) such that chaotic trajectories are guaranteed to occur. There are, however, many well known systems in the literature known to exhibit chaotic trajectories. For instance, the Lorenz system [45], Duffing [46] and van der Pol [47] oscillators, the Rössler [48] and Chua [49] systems, and neural oscillators [50] (e.g., the Hodgkin-Huxley, Morris-Lecar, Hindmarsh-Rose, and FitzHugh-Nagumo oscillators). We can use any of these chaotic systems (if they satisfy all the required extra conditions, see Section 6.4) as driver and then select a pair of responders with convergent dynamics. Indeed, we need to verify that the responders that we choose produce chaotic trajectories when driven by the chaotic driver. Moreover, to generate the pseudorandom realizations from $p_{V}^{*}(v)$ (this is addressed in the next section), we need the chaotic trajectories of the responders, regarded as a random process, to be stationary, i.e., after transients have settled down, trajectories must follow a stationary probability distribution [39] which is independent of the initial conditions. The latter is a strong condition that is not satisfied for all chaotic systems. The existence of stationary distributions for chaotic trajectories has been proven for hyperbolic and quasi-hyperbolic (also called singular-hyperbolic) chaotic systems [15]. The definition of (quasi) hyperbolic dynamical systems [15, 51] is technical and not needed for the subsequent results. It requires concepts from differential topology that we prefer to omit here for readability of the manuscript. It suffices to know that the chaotic system that we use for the driver must lead to stationary distributions of the responders. This can be tested numerically by Monte Carlo simulations [21]. Moreover, there are many well-known chaotic systems with (quasi) hyperbolic dynamics in the literature, e.g., the Lorenz and Chua systems [52], neural oscillators [53], the many predator-pray like systems given in [54, 55], and some mechanical nonlinear oscillators [56]. In the next subsection, we provide a synthesis procedure to choose the functions of the coupled systems (4)-(6) such that all the required conditions mentioned above are satisfied.

6.4 General Guidelines

**Synthesis Procedure:

** 1) Select a driver dynamics (5) (i.e., the vector field $d(\xi)$) known to be chaotic and (quasi) hyperbolic (e.g., systems in [52]-[56]).

2) Verify that the corresponding $d(\xi)$ is locally Lipschitz and the trajectories of the driver are globally bounded, in the sense of Definition 3, using, e.g., [41, Theorem 4.18].

3) In (5), let $\xi=(\xi^{1},\ldots,\xi^{n_{\xi}})^{\top}\in{\mathds{R}}^{n_{\xi}}$, $\xi^{i}\in{\mathds{R}}$, and $u(t)=l(\xi(t))=\xi^{j}(t)$, $i,j\in\{1,\ldots,n_{\xi}\}$, i.e., fix the output of the driver to be any state of (5). In doing this, we ensure that $u(t)$ is continuous, bounded, chaotic, and (quasi) hyperbolic.

4) For the responders (4) and (6), select any continuously differentiable vector field $r(\zeta,u)$ (with respect to $\zeta$) leading to ISS dynamics, see Remark 6, and satisfying the conditions for convergence in Corollary 1, e.g., $r(\zeta,u)=A\zeta+\psi(u)$, for any matrix $A\in{\mathds{R}}^{n_{\zeta}\times n_{\zeta}}$ with spectral radius $\rho[A]<1$ and differentiable vector field $\psi:{\mathds{R}}^{n_{u}}\rightarrow{\mathds{R}}^{n_{\zeta}}$. Then, we ensure that the responders have bounded trajectories and exponentially synchronize.

5) Verify that the trajectories of the responders, when driven by the chaotic driver, are chaotic (using Lyapunov exponents or the (0-1) test) and, after transients have settled down, lead to a stationary probability distribution independent of the initial conditions. See Section 6.3 for details.

6) In (4) (and respectively in (6)), let $\zeta_{1}=(\zeta_{1}^{1},\ldots,\zeta_{1}^{n_{\zeta}})^{\top}\in{\mathds{R}}^{n_{\zeta}}$, $\zeta_{1}^{i}\in{\mathds{R}}$, and $s_{1}(t)=l(\zeta_{1}(t))=\zeta_{1}^{j}(t)$, $i,j\in\{1,\ldots,n_{\xi}\}$, i.e., fix the output of the responders to be any state of (4) and (6), respectively. Indeed, we need the same $j$ for both responders, i.e., $s_{1}(t)=\zeta_{1}^{j}(t)$ and $s_{2}(t)=\zeta_{2}^{j}(t)$. In doing this, we ensure that $s_{1}(t)$ and $s_{2}(t)$ are continuous, bounded, chaotic, and lead to stationary probability distributions.

7 Generation of Optimal Pseudorandom Numbers

In this section, we assume that the driver and the responders dynamics have been designed following the general guidelines in Section 6.4. Then, for sufficiently large $t$, the chaotic trajectories of the responders are practically synchronized, i.e., for any finite $t^{*}\in{\mathds{R}}_{>0}$, there is $\epsilon_{t^{*}}\in{\mathds{R}}_{>0}$, such that $|s_{1}(t)-\phi^{s}_{u}(t,u(t))|\leq\epsilon_{t^{*}}$ and $|s_{2}(t)-\phi^{s}_{u}(t,u(t))|\leq\epsilon_{t^{*}}$, for all $t\geq t^{*}$, where $\phi^{s}_{u}(t,u(t))\in\mathcal{S}\subset{\mathds{R}}$ denotes the asymptotic synchronous solution for some compact set $\mathcal{S}$; and samples from $\phi^{s}_{u}(t,u(t))$ follow a stationary probability distribution. Here, we assume that the responders have been operating for sufficiently large time such that the synchronization error, $|s_{1}(t)-s_{2}(t)|$, is so small that trajectories of the responders can be assumed to be equal to $\phi^{s}_{u}(t,u(t))$ (see Remark 2), i.e., $t^{*}$ is sufficiently large so that $\epsilon_{t^{*}}$ is practically zero. In Section 7.1, we quantify the worst-case distortion induced by assuming $s_{1}(t)=s_{2}(t)=\phi^{s}_{u}(t,u(t))$ in finite time. In particular, we give an upper bound on the mean squared error $E[|Y-\hat{Y}|^{2}]$, where $\hat{Y}$ denotes the estimate of realizations of $Y$ using $s_{1}(t)$, $s_{2}(t)$, and the algorithm provided below. In the remainder of this section, we assume $s_{1}(t)=s_{2}(t)=\phi^{s}_{u}(t,u(t))$. Note that the sample space of $\phi^{s}_{u}(t,u(t))$, regarded as a random process, is some compact set $\mathcal{S}\subset{\mathds{R}}$, i.e., the sample space is a subset of the real line and thus samples from $\phi^{s}_{u}(t,u(t))$ follow some stationary probability density function (pdf), say $f_{S}(s)$, for some virtual continuous random variable $S$. That is, for $s(t):=\phi^{s}_{u}(t,u(t))$, define the sampled sequence $s_{k}:=s(t_{k})$ for sampling time-instants $t_{k}\in{\mathds{R}}_{>0}$, $t_{k}:=\Delta k$, $k\in{\mathds{N}}$, and sampling period $\Delta\in{\mathds{R}}_{>0}$; then, $s_{k}\sim f(s)$ for all $k$. Because we know the dynamics (4)-(6), we can obtain $f_{S}(s)$ by Monte Carlo simulations [21]. If we know $f_{S}(s)$, we can always find a set of cells $C:=\{c^{1},\ldots,c^{M}\}$, $M\in{\mathds{N}}$, $j\in\{1,\ldots,M\}$, such that $\bigcup_{j}c^{j}={\mathds{R}}$, $\bigcap_{j}c^{j}=\emptyset$, and $\text{Pr}[s_{k}\in c]=\text{Pr}[V=v]=p_{V}^{*}(v)$ for $v\in\mathcal{Y}$ and $c\in C$. In other words, using the pdf $f_{S}(s)$, we can select the cells $C$ so that the probability that $s_{k}$ lies in the cells equals the optimal probability distribution $p_{V}^{*}(v)$. It follows that we can generate pseudorandom realizations from $p_{V}^{*}(v)$ by properly selecting $C$. Note that, because realizations are being generated by a deterministic process, there would be high correlation between consecutive realizations for small sampling period $\Delta$. However, because the $s_{k}$ is a stationary process (see Section 6.3), the larger the $\Delta$, the smaller the correlation between $s_{k}$ and $s_{k+1}$ for all $k\in{\mathds{N}}$. Indeed, large $\Delta$ would introduce large time-delays for generating realizations. There is a trade-off between correlation and time-delay that should be taken into account in practice. One way to deal with this trade-off is to compute the normalized autocorrelation function [15, 20] of $s_{k}$. Then, we select the smallest time-delay $\tau\in{\mathds{N}}$ that leads to a desired correlation between $s_{k}$ and $s_{k+\tau}$, $k\in{\mathds{N}}$, and use the delayed sequence $s^{\tau}(\cdot):=\{s_{k},s_{k+\tau},s_{k+2\tau},\ldots\}$ to generate realizations from $p_{V}^{*}(v)$. In the following algorithm, we summarize the ideas introduced above.

**Algorithm 1: Pseudorandom Number Generation:

** 1) Consider the probability mass function $p_{V}^{*}(v)=\text{Pr}[V=v]$, $v\in\mathcal{Y}=\{y_{1},\ldots,y_{M}\}$, solution to Problem 1; and the synchronous solution $s(t)=\phi^{s}_{u}(t,u(t))$ of the responders.

2) Fix the sampling period $\Delta\in{\mathds{R}}_{>0}$ and obtain, by Monte Carlo simulations [21], the probability density function $f_{S}(s_{k})$ of the sampled sequence $s_{k}=s(t_{k})$, $t_{k}=\Delta k$, $k\in{\mathds{N}}$.

3) Select a finite set of cells $C=\{c^{1},\ldots,c^{M}\}$, $M\in{\mathds{N}}$, $j\in\{1,\ldots,M\}$, such that $\bigcup_{j}c^{j}={\mathds{R}}$, $\bigcap_{j}c^{j}=\emptyset$, and $\text{Pr}[s_{k}\in c^{j}]=\text{Pr}[V=y_{j}]$ for all $y_{j}\in\mathcal{Y}$.

4) Generate realization from $p_{V}^{*}(v)$ using the piecewise function:

[TABLE]

7.1 Distortion Induced by Synchronization Errors

Algorithm 1 in Section 7 is constructed under the assumption that responders are perfectly synchronized. However, because we only have exponential synchronization, in finite time, there is always a “small” difference between $s_{1}(t)$ and $s_{2}(t)$ due to potentially different initial conditions. It follows that there is also a difference between realizations generated using $s_{1}(t_{k})$, denoted as $v_{k}^{1}\in\mathcal{Y}$, and realizations $v_{k}^{2}\in\mathcal{Y}$ generated through $s_{2}(t_{k})$, where $\mathcal{Y}=\{y_{1},\ldots,y_{M}\}$. Exponential synchronization implies that for any finite $t^{*}\in{\mathds{R}}_{>0}$, there is $\delta(t^{*},|s_{1}(0)-s_{2}(0)|)\in{\mathds{R}}_{>0}$ (denoted as $\delta_{t^{*}}$ for simplicity), parametrized by $t^{*}$ and the initial synchronization error $|s_{1}(0)-s_{2}(0)|$, such that $|s_{1}(t_{k})-s_{2}(t_{k})|\leq\delta_{t^{*}}$ for all $t_{k}\geq t^{*}_{k}$, and $\lim_{k\rightarrow\infty}|s_{1}(t_{k})-s_{2}(t_{k})|=0$. Consider the cell $c^{j}$, $c^{j}\in C$, with end points $c^{j}_{1}$ and $c^{j}_{2}$, $c^{j}_{1}<c^{j}_{2}$, the length of $c^{j}$ is defined as $l(c^{j}):=c^{j}_{2}-c^{j}_{1}$. If $c^{j}_{1}=\pm\infty$ (or $c^{j}_{2}=\pm\infty$), $l(c^{j})=\infty$. Without loss of generality, let $l(c^{2})\leq l(c^{3})\leq\ldots\leq l(c^{M-1})$, $l(c^{1})=\infty$, and $l(c^{M})=\infty$. Note that, if $\delta_{t^{*}}\leq l(c^{2})$, $v_{k}^{1}$ and $v_{k}^{2}$ are at most one level apart from each other, e.g., if $v_{k}^{1}=y_{1}$, then either $v_{k}^{2}=y_{1}$ or $v_{k}^{2}=y_{2}$; and if $v_{k}^{1}=y_{3}$, then $v_{k}^{2}=y_{2}$, $v_{k}^{2}=y_{3}$, or $v_{k}^{2}=y_{4}$. It follows that $p_{\hat{Y}|Y}(\hat{y}|y)$, $y,\hat{y}\in\mathcal{Y}$, is of the form depicted in Figure 3, where $\hat{Y}$ denotes the estimate of realizations of $Y$ using $s_{1}(t_{k})$, $s_{2}(t_{k})$, and Algorithm 1. Similarly, if $l(c^{2})<\delta_{t^{*}}\leq l(c^{3})$, $v_{k}^{1}$ and $v_{k}^{2}$ are at most two levels apart from each other and thus lead to a different structure of the transition probabilities. Here, we only consider the case where $\delta_{t^{*}}\leq l(c^{2})$. Distortion induced by larger synchronization errors can be estimated following the same methods. Note that, because responders synchronize exponentially, as $\delta_{t^{*}}\rightarrow 0$ ($t^{*}\rightarrow\infty$), $p_{\hat{Y}|Y}(\hat{y}|y)\rightarrow 1$ for $\hat{y}=y$, and $p_{\hat{Y}|Y}(\hat{y}|y)\rightarrow 0$, for $\hat{y}\neq y$, for all $y,\hat{y}\in\mathcal{Y}$. That is, distortion due to synchronization errors disappears exponentially fast. The actual value of the transition probabilities depend on the responders and driver dynamics, the initial conditions, and the cells $C$. However, we do not need these probabilities, only the structure of $p_{\hat{Y}|Y}(\hat{y}|y)$ depicted in Figure 3 is used to derive an upper bound on the expected distortion. Let $\mathcal{V}_{\delta}\subseteq\mathcal{Y}\times\mathcal{Y}$ denote the set of pairs $(y_{j},y_{i})$ for which there is a nonzero transition probability $p_{\hat{Y}|Y}(y_{j}|y_{i})$ between $Y=y_{j}$ and $\hat{Y}=y_{i}$, $y_{j},y_{i}\in\mathcal{Y}$, as depicted in Figure 3. The set $\mathcal{V}_{\delta}$ is parametrized by the upper bound on the synchronization error $|s_{1}(t_{k})-s_{2}(t_{k})|\leq\delta_{t^{*}}\leq l(c^{2})$. Define the distortion function $d(Y,\hat{Y}):=|Y-\hat{Y}|^{2}$. The function $d(Y,\hat{Y})$ is a deterministic function of two jointly distributed random vectors, $Y$ and $\hat{Y}$, with joint distribution $p_{Y,\hat{Y}}(y,\hat{y})=p_{Y}(y)p_{\hat{Y}|Y}(\hat{y}|y)$. Hence, see [39] for details, we can write the expected distortion as follows

[TABLE]

where the left-hand side of (11) follows from the definition of $\mathcal{V}_{\delta}$ above, and the last inequality from the fact that $p_{\hat{Y}|Y}(\hat{y}|y)\leq 1$ for all $y,\hat{y}\in\mathcal{Y}$. The constant $\bar{d}_{\delta}\in{\mathds{R}}_{>0}$ provides an upper bound on the worst-case distortion induced by a $\delta_{t^{*}}$ synchronization error. Moreover, as $\delta_{t^{*}}\rightarrow 0$, $\mathcal{V}_{\delta}\rightarrow\{(y_{1},y_{1}),(y_{2},y_{2}),\ldots,(y_{M},y_{M})\}$; therefore, $\lim_{\delta_{t^{*}}\rightarrow 0}\bar{d}_{\delta}=0$. That is, distortion due to synchronization errors is bounded by $\bar{d}_{\delta}$ and vanishes exponentially fast.

8 Simulation Results

We next present an evaluation of our algorithms on real data. We use the adult-dataset, available from the UCI Machine Learning Repository [57], which contains census data. Each attribute within the dataset has $3.9\times 10^{4}$ entries. We use three of these attributes: race, sex, and income, which take values on finite discrete sets. We let race and sex be the private information, $X$, and use income as the information requested by the query, $Y$. The probability mass functions of $X$ and $Y$, and part of the one of $(X,Y)$ are given in Table 1.In Figure 4, we depict $p_{X}(x)$, $p_{Y}(y)$, and $p_{X,Y}(x,y)$ with mass points indexed in the order given in Table 1.We first compute the optimal distribution $p_{V}^{*}(v)$ of the distorting additive noise $V$. We solve the convex program (8) in Theorem 1. The optimal distribution is depicted in Figure 5 and the corresponding numerical values are given in Table 2. This $p_{V}^{*}(v)$ leads to $I[X;Y+V]=0.0024$ while the mutual information without distortion is $I[X;Y]=0.0251$, i.e., according to our metric, by optimally distorting the query, we leak about ten times less information. To generate realization from this distribution at both sides of the channel, we use trajectories of two chaotic responders as introduced in Section 2. We use the synthesis procedure in Section 6.4 to select suitable driver and responders. As driver (5), we use the Lorenz system:

[TABLE]

with states $\xi_{1},\xi_{2},\xi_{3}\in{\mathds{R}}$ and driving signal $u\in{\mathds{R}}$. The Lorenz system produces bounded trajectories [58], and is known to be chaotic and quasi-hyperbolic [52]. For the responders $\eqref{2}$ and $\eqref{4}$, we let $r(\zeta,u)=A\zeta+\psi(u)$, with $A=\text{diag}[-1,-2.5]$ and $\psi(u)=(-5u^{2},50\sin(u))^{\top}$. Because $A$ is diagonal and has negative eigenvalues, responders satisfy the conditions of Corollary 1 with $P=I_{2}$; hence, they are convergent systems and thus exponentially synchronize when driven by the same input $u(t)$. Moreover, since responders are linear in $\zeta$ and $A$ is Hurwitz, systems can be proved to be ISS with input $\psi(u)$ [42]. Because $u$ is bounded and $\psi(u)$ is continuous, by the extreme value theorem, $\psi(u)$ is bounded, which, together with ISS, imply boundedness of the responders’ trajectories [42]. We let the outputs of the responders be $s_{1}(t)=\zeta_{1}^{2}$ and $s_{2}(t)=\zeta_{2}^{2}$ (their second state). In Figure 6, we show traces of the chaotic driver and responders trajectories obtained by computer simulations (using Matlab from Mathworks), and in Figure 7, we plot the synchronization error between the outputs of the responders. We initialized the responders in antiphase $\zeta_{1}(0)=-\zeta_{2}(0)=(150,150)^{\top}$, and far from the limit trajectory. Note, in Figure 7, that responders synchronize exponentially and are practically synchronized for $t\geq 5$. Moreover, after $t\geq 14$, the synchronization error is within Matlab’s precision ($10^{-12}$). Because the Lorenz system is quasi-hyperbolic, samples from the driving signal $u(t)$ follow a stationary distribution that is independent of the initial conditions of the driver, see Section 6.3. Then, according to the synthesis procedure in Section 6.4, we next verify, using Monte Carlo simulations, that samples $s_{k}=s(t_{k})$ (see Section 7), from the synchronous trajectory, $s_{1}(t)=s_{1}(t)=s(t)$, are also stationary. To do so, we compute the probability density function $f_{S}(s)$, $s_{k}\sim f_{S}(s)$, for different initial conditions and verify that all of them lead to the same density. In Figure 8, we depict probability densities of $s_{k}$ for twenty different initial conditions, sampling instants $t_{k}=\Delta k$, $\Delta=0.001$, and $t\in[0,4000]$. Note that they all lead to the same density $f_{S}(s)$. The support (obtained numerically) of $f_{S}(s)$ is given $\mathcal{S}=[-10.8585,10.8683]$. Finally, we use the piecewise function (10) to generate realizations from $p_{V}^{*}(v)$ using samples, $s_{k}$, from the synchronous trajectory. Following the algorithm given in Section 7, we have to divide the support $\mathcal{S}$ of $f_{S}(s)$ into a set of partitions $C=\{c^{1},\ldots,c^{M}\}$, such that the probability that $s_{k}$ lies in the cells equals the optimal probability distribution $p_{V}^{*}(v)$. This can be done using the empirical Cumulative Distribution Function (CDF), $F_{S}(s)$, corresponding to $f_{S}(s)$. We depict this CDF in Figure 9. Then, we simply select the cells $C$ such that $p_{V}^{*}(y_{i})=\text{Pr}[V=y_{i}]=\text{Pr}[c^{i}\leq S\leq c^{i+1}]=F_{S}(c^{i+1})-F_{S}(c^{i})$ for all $i\in\{1,\ldots,M-1\}$, $M=9$ (the cardinality of the alphabet of $Y$). For this CDF and $p_{V}^{*}(v)$ in Table 2, we obtain the following cells:

[TABLE]

In Figure 10, we show realizations generated by the piecewise function (10) at both sides of the channel, and the corresponding probability mass functions. To generate this realizations, at the trusted server, we use samples from $s_{1}(t)$ and, at the remote station, we sample $s_{2}(t)$. Note that, as expected, all samples are perfectly synchronized and their probability mass functions are equal to $p_{V}^{*}(v)$ in Figure 5.

9 Conclusions

Using an information-theoretic privacy metric (mutual information), we have provided a general privacy framework based on additive distorting random vectors and exponential synchronization of chaotic systems. The synthesis of the optimal probability distribution, $p^{*}_{V}(v)$, of the additive distorting vector $V$ has been posed as a convex program in $p_{V}(v)$. We have provided an algorithm for generating pseudorandom realizations from this distribution using trajectories of chaotic oscillators. To generate equal realizations at both sides of the channel, we have induced exponential synchronization on two chaotic oscillators (one at each side of the channel), and use their trajectories and the proposed algorithm to generate realizations. However, exponential synchronization implies that, in finite time, there is always a small error between trajectories (and thus also between realizations). We have derived an upper bound on the worst-case distortion induced by finite-time synchronization errors and showed that this distortion disappears exponentially fast. Using off-the-shelf results in the literature, we have provided general guidelines for selecting the dynamics of the responders and driver so that our algorithm for generating synchronized realizations from $p^{*}_{V}(v)$ is guaranteed to work. We have presented simulation results to illustrate our results.