# Privacy-Preserving Anomaly Detection in Stochastic Dynamical Systems:
Synthesis of Optimal Gaussian Mechanisms

Haleh Hayati, Carlos Murguia, Nathan van de Wouw Haleh Hayati, Carlos Murguia,
and Nathan van de Wouw are with the Department of Mechanical Engineering,
Dynamics and Control Group, Eindhoven University of Technology, The
Netherlands. Emails: &<EMAIL_ADDRESS>&<EMAIL_ADDRESS>&
<EMAIL_ADDRESS>
###### Abstract

We present a framework for the design of distorting mechanisms that allow
remotely operating anomaly detectors in a privacy-preserving fashion. We
consider the problem setting in which a remote station seeks to identify
anomalies using system input-output signals transmitted over communication
networks. However, in such a networked setting, it is not desired to disclose
true data of the system operation as it can be used to infer private
information – modeled here as a system private output. To prevent accurate
estimation of private outputs by adversaries, we pass original signals through
distorting (privacy-preserving) mechanisms and send the distorted data to the
remote station (which inevitably leads to degraded monitoring performance).
The design of these mechanisms is formulated as a privacy-utility (tradeoff)
problem where system utility is characterized by anomaly detection
performance, and privacy is quantified using information-theoretic metrics
(mutual information and differential entropy). We cast the synthesis of
dependent Gaussian mechanisms as the solution of a convex program (log-
determinant cost with linear matrix inequality constraints) where we seek to
maximize privacy over a finite window of realizations while guaranteeing a
bound on monitoring performance degradation. We provide simulation results to
illustrate the performance of the developed tools.

##  I Introduction

Scientific and technological advancements in today’s hyperconnected world have
resulted in a massive amount of user data being collected and processed by
hundreds of companies over public networks. Companies use this data to provide
personalized services and targeted advertising. However, these technologies
have come at the cost of a significant loss of privacy in society. Depending
on their resources, adversaries can infer sensitive (private) information
about users/systems from public data available on the internet and/or
unsecured networks/servers. An example of this privacy loss is the potential
use of data from smart electrical meters by criminals, advertising agencies,
and governments to monitor occupants’ activities [1, 2]. Another example is
the privacy loss caused by information sharing in distributed control systems
and cloud computing [3]. This is why researchers from different fields (e.g.,
information theory, computer science, and control engineering) have been drawn
to the broad research subject of privacy and security of Cyber-Physical
Systems (CPSs), [4]-[13].

In most engineering applications, information about the state of systems, say
$x$, is obtained through sensor measurements and then sent to a remote station
through communication networks for signal processing and decision-making
purposes. One such application is remotely operating anomaly detection, where
a remote station seeks to identify anomalies based on system data transmitted
over communication networks. If the communication network is public/unsecured
and/or the remote station is not trustworthy, adversaries might access and
estimate the system state, reference signals, product specifications, etc. To
avoid an accurate estimation, before transmission, we randomize sensor data
and input signals using additive dependent Gaussian random vectors and send
the corrupted data to the remote station instead. These vectors are designed
to hide the private part of the state – a private output, say $s$, modeled as
some linear function of the system state, $s=Dx$, for some deterministic
matrix $D$.

The aim of this work is to devise synthesis tools that allow designing such
distorting mechanisms that protect private information but, at the same time,
allow for the implementation of remote anomaly detection procedures – here, we
consider standard Kalman filter-based chi-squared change detection procedures
[14]. Note that by injecting privacy-preserving noises, the uncertainty level
in input-output data increases for both the anomaly monitor at the remote
station and eavesdropping adversaries since the privacy noise enters the
system dynamics in a similar fashion to disturbances. The impact of
uncertainty on the detectability of model-based fault detectors is well
studied in the literature [15]. It follows that when designing privacy-
preserving noises, we want to maximize uncertainty in the direction of private
outputs and minimize it in that of the anomaly detection procedure. In
general, stochastic additive noise increases the false alarm rate of the chi-
squared procedure. Therefore, there is a trade-off between privacy and
monitoring performance that must be taken into account in the design. We
remark that we protect against both eavesdropping adversaries and the
(potentially untrustworthy) remote station. We aim to develop a privacy scheme
that considers no trusted party and seeks to destroy information in the
direction of private outputs as much as possible while keeping monitoring
performance degradation bounded.

We follow an information-theoretic approach to privacy. As _privacy metric_ ,
we propose a combination of _mutual information_ and _entropy_ [16] between
disclosed and private data. To characterize anomaly detection performance, we
consider the false alarm rate of the chi-squared procedure as _performance
metric_. Hence, the synthesis algorithm seeks to optimize our privacy metric
by selecting the stochastic properties of the injected dependent noise (their
joint covariance matrices over a finite window of realizations) while not
degrading the false alarm rate of the monitor by more than desired (allowed)
amount. We prove that this synthesis problem can be posed as a constrained
convex optimization problem – log-determinant cost with linear matrix
inequality constraints.

Using additive random noise is common practice to enforce privacy of sensitive
data. Differential Privacy (DP) is a popular method in the context of privacy
of databases [7, 17]. In differential privacy, because it provides certain
(not optimal) privacy guarantees, Laplacian noise is commonly used [18].
However, most results using differential privacy do not consider the trade-off
between privacy and utility, i.e., they do not account for the negative effect
of the privacy mechanism on the application performance (which leads to
conservative results when data utility is of high importance). In general,
when considering privacy-utility trade-offs, optimal additive noise
distributions are determined by the particular privacy and distortion metrics
being considered, and the system configuration [19]-[21].

There are also techniques addressing privacy in dynamical systems from an
information-theoretic perspective [1, 4, 5, 22, 9, 23]. In this line of work,
privacy is characterized using information-theoretic metrics – e.g., mutual
information, entropy, Kullback-Leibler divergence, and Fisher information.
However, independently of the metric being used, if the data to be kept
private follows continuous probability distributions, the problem of finding
the optimal additive noise to maximize privacy is difficult to solve [4, 24].
This issue has been addressed by assuming the data to be kept private is
deterministic [4, 20, 10]. However, in a Cyber-Physical-Systems context, the
inherent system dynamics and unavoidable system and sensor noise lead to
stochastic non-stationary data, and thus, existing tools do not fit this
setting. Recently, the authors in [13, 25] have proposed a framework for
synthesizing optimal transition probabilities to maximize privacy for a class
of quantized CPSs with discrete multivariate probability distributions.
Although this framework is general and leads to distorting mechanisms with
arbitrary distributions, the computational complexity induced by exploring all
the possible transition probabilities from private to disclosed distorted data
is high and increases exponentially with the alphabet of the private data. In
[4], the authors neglect quantization and work with dynamical systems driven
by continuous (Gaussian) noises. They prove that the optimal noise
distribution minimizing the Fisher information is Gaussian for unconstrained
additive noise cases (so data distortion/utility is not considered). This
observation has also been made in [24] when using mutual information as a
privacy metric. In [26], we have started exploring these ideas assuming
Gaussianity in all system disturbances and without considering any particular
application at the remote station – we mainly focused on formulating the
privacy-utility problem over a finite horizon and characterized utility
degradation by the expected squared distortion between distorted and original
data.

Motivated by these results, in this manuscript, we propose a general class of
Gaussian mechanisms and seek the optimal parametrization of the additive
dependent noise distribution (its joint covariance matrix over a window of
realizations) given the system dynamics and the anomaly detection procedure at
the remote station. Compared to [26], the contributions of this manuscript are
the following: 1) We propose an optimization-based framework for synthesizing
optimal dependent Gaussian mechanisms to hold a desired trade-off between
privacy and anomaly detection performance at the remote station over finite
horizons. 2) We implement the proposed mechanisms for an arbitrarily large
horizon by solving finite horizon problems repeatedly and then, sequentially
implementing the obtained mechanisms, and prove that the amount of information
leakage in this large horizon is upper bounded by the sum of the information
leaked in each finite time horizon (which allows for the synthesis of
arbitrary large horizon sub-optimal distorting mechanisms by repeatedly
solving the associated finite horizon optimization problems). 3) We prove that
we can allow for prior distributions in the system dynamics (the joint
probability distributions between private and disclosed data) that are not
necessarily Gaussian, as long as they are log-concave [27]. The latter further
generalizes the class of stochastic dynamical systems that our framework can
deal with. To the best of the authors’ knowledge, there is only one related
result in the literature. For remote anomaly detection in stochastic dynamical
systems, the authors in [28] consider differential privacy and provide tools
to select the variance of i.i.d. additive noise so that disclosed data
guarantees a certain level of differential privacy. Their results consider
privacy of input signals (so no private state-dependent outputs), and at every
time-step, i.e., they do not consider information disclosed over a history of
observations. While this indeed gives some privacy guarantees, it is well
known that for dynamical systems, the more data is collected, the more
accurate the input/state estimation can be. For instance, in the deterministic
case, if the system is observable and of dimension $n$, after $n$ time-steps,
the state can be completely reconstructed without error. Therefore, if the
privacy scheme is to avoid input/state estimation, giving guarantees at every
time-step is certainly not sufficient. Moreover, making noise i.i.d. at every
time-step is also conservative as there is an obvious correlation among
measurements and inputs in dynamical systems. Taking into account these
points, our contributions with respect to [28] are the following: 1)
mechanisms that provide optimal trade-offs between privacy and anomaly
detection performance (not only worst-case guarantees based on differential
privacy); 2) an optimization-based synthesis framework that casts the design
of privacy-preserving mechanisms as the solution of convex programs where
privacy is maximized for a desired level of anomaly detection; 3) results that
consider privacy over a history of observation (finite horizon privacy) which
is key for dynamical systems; 4) dependent noise distributions over the
history of observations, so we design mechanisms that take into account the
correlation of data induced by the system dynamics; and 5) the notion of
private outputs which provide freedom to choose the specific part of the
system that must be kept private.

##  II Notation and Preliminaries

###  II-A Notation

The symbol ${\mathds{R}}$ stands for the real numbers, ${\mathds{R}}^{+}$
denotes the set of positive real numbers. The symbol ${\mathds{N}}$ stands for
the set of natural numbers. The Euclidian norm in ${\mathds{R}}^{n}$ is
denoted by $||X||$, $||X||^{2}=X^{\top}X$, where ⊤ denotes transposition. The
$n\times n$ identity matrix is denoted by $I_{n}$ or simply $I$ if $n$ is
clear from the context. Similarly, $n\times m$ matrices composed of only ones
and only zeros 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. The notation
$\operatorname{col}\left(x_{1},\ldots,x_{n}\right)$ stands for the column
vector composed of the elements $x_{1},\ldots,x_{n}$. This notation will also
be used in case the components $x_{i}$ are vectors. For positive definite
(semidefinite) matrices, we use the notation $P>0$ ($P\geq 0$). The notation
$X\sim\mathcal{N}[\mu,\Sigma^{X}]$ means $X\in{\mathds{R}}^{n}$ is a normally
distributed random vector with mean $E[X]=\mu\in{\mathds{R}}^{n}$ and
covariance matrix $E[(X-\mu)(X-\mu)^{T}]=\Sigma^{X}\in{\mathds{R}}^{n\times
n}$, where $E[a]$ denotes the expected value of the random vector $a$. The
notation $Y\sim\Gamma(K,\theta)$ means $Y$ follows a gamma distribution with
shape parameter $K$ and scale parameter $\theta$. Finite sequences of vectors
are written as
$X^{K}:=(X(1)^{\top},\ldots,X(K)^{\top})^{\top}\in{\mathds{R}}^{Kn}$,
$X(i)\in{\mathds{R}}^{n}$, $i\in\\{1,\ldots,K\\}$, and $n,K\in{\mathds{N}}$.
To avoid confusion, we denote powers of matrices as $(A)^{K}=A\cdots A$ ($K$
times) for $K>0$, $(A)^{0}=I$, and $(A)^{K}=\mathbf{0}$ for $K<0$. The
operators $\log(\cdot)$, $\det(\cdot)$, and $\text{tr}(\cdot)$ are logarithm,
determinant, and trace, respectively.

###  II-B Preliminaries

In this section, we present some definitions and results needed in subsequent
sections.

######  Definition 1 (MMSE Estimator [29]).

Let $X$ and $Y$ be two jointly Gaussian random vectors. The Minimum Mean
Square Error (MMSE) estimate of $X$ given $Y$ is given by _:_

$\displaystyle\hat{X}=\Sigma^{XY}{\Sigma^{Y}}^{-1}\left(Y-E(Y)\right)+E(X),$
(1)

where $\Sigma^{XY}$ denotes the cross-covariance matrix between $X$ and $Y$
and $\Sigma^{Y}$ is the covariance matrix of $Y$.

######  Definition 2 (Differential Entropy [16]).

Consider the random vector, $S\in{\mathds{R}}^{n}$,
$S\sim\mathcal{N}[\mu^{S},\Sigma^{S}]$. Its differential entropy, $h[S]$, can
be written in terms of its covariance matrix as _:_

$\displaystyle
h[S]=\frac{1}{2}\log\det\left(\Sigma^{S}\right)+\frac{n}{2}+\frac{n}{2}\log(2\pi).$
(2)

Differential entropy is a measure of the average uncertainty in a random
vector. We use base two $\log(\cdot)$, so entropy is given in bits.

######  Definition 3 (Mutual Information [16]).

Consider two jointly distributed continuous random vectors, $S$ and $Z$, with
joint entropy $h[S,Z]$ and marginal entropies $h[S]$ and $h[Z]$. Their mutual
information, $I[S;Z]$, is given as _:_

$\displaystyle I[S;Z]=h[S]+h[Z]-h[S,Z].$ (3)

The mutual information between two jointly distributed vectors is a measure of
the statistical dependence between them.

##  III System Description and Anomaly Detection

###  III-A System Description

We consider discrete-time stochastic systems of the form:

$\displaystyle\left\\{\begin{aligned}
x_{k+1}&=Ax_{k}+Bu_{k}+t_{k}+G\delta_{k},\\\
y_{k}&=Cx_{k}+w_{k}+H\delta_{k},\\\ s_{k}&=Dx_{k},\end{aligned}\right.$ (4)

with time-index $k\in{\mathds{N}}$, state $x_{k}\in{\mathbb{R}^{{n_{x}}}}$,
measurable output $y_{k}\in{\mathbb{R}^{{n_{y}}}}$, _known_ input
$u_{k}\in{\mathbb{R}^{{n_{u}}}}$, private performance output
$s_{k}\in{\mathbb{R}^{{n_{s}}}}$, and matrices $(A,B,C,D,G,H)$ of appropriate
dimensions. Vector $\delta_{k}\in{\mathbb{R}^{{n_{\delta}}}}$ models changes
in the system dynamics due to faults or anomalies (thus, in the absence of
anomalies $\delta_{k}=\mathbf{0}$). Matrix $D$ is full row rank, and the pair
$(A,C)$ is detectable. The state disturbance $t$ and the output disturbance
$w$ are multivariate i.i.d. Gaussian processes with zero mean and positive
definite covariance matrices ${\Sigma^{t}}$ and ${\Sigma^{w}}$, respectively.
The initial state $x_{1}$ is assumed to be a Gaussian random vector with
$E[x_{1}]=\mu^{x}_{1}\in\mathbb{R}^{n_{x}}$ and covariance matrix
$\Sigma^{x}_{1}\in\mathbb{R}^{n_{x}\times n_{x}}$, $\Sigma^{x}_{1}>0$.
Disturbances $t_{k}$ and $w_{k}$ and the initial condition $x_{1}$ are
mutually independent. We assume that matrices (vectors) $(A,$
$B,C,D,\Sigma^{x}_{1},\mu^{x}_{1},\Sigma^{t},\Sigma^{w})$ and input $u_{k}$
are known for all $k$.

###  III-B Anomaly Detection

The aim of the anomaly detection algorithm is to identify anomalies (i.e.,
when $\delta_{k}\neq\mathbf{0}$ in (4)) using system input-output signals.
Here, we consider standard Kalman filter-based chi-squared change detection
procedures [14]. The main idea is to use an estimator to anticipate the system
evolution in the absence of anomalies. This prediction is subsequently
compared with actual measurements from sensors. If the difference between what
is measured and the estimation (commonly referred to as residuals) is larger
than expected, there might be an anomaly in the system. We use steady-state
Kalman filters as the state estimator.

####  III-B1 Steady-state Kalman Filter (anomaly-free case)

Consider the following one step-ahead Kalman filter [30] for (4):

$\hat{x}_{k+1}=A\hat{x}_{k}+Bu_{k}+L(y_{k}-C\hat{x}_{k}),$ (5)

with estimated state $\hat{x}_{k}\in{\mathbb{R}^{{n_{x}}}}$,
$\hat{x}_{1}=E[x_{1}]$, and output injection gain matrix
$L\in{\mathbb{R}^{{n_{x}}\times{n_{y}}}}$. Define the estimation error
$e_{k}:=x_{k}-\hat{x}_{k}$. The optimal filter gain $L$ minimizing the trace
of the asymptotic estimation error covariance matrix
$P:=\lim_{k\rightarrow\infty}E[e_{k}e_{k}^{\mathrm{T}}]$ in the absence of
anomalies is given by

$L:=\left(APC^{\top}\right)\left(\Sigma^{w}+CPC^{\top}\right)^{-1},$ (6)

where $P$ is the solution of the Riccati equation:

$APA^{\top}-P+\Sigma^{t}=APC^{\top}\left(\Sigma^{w}+CPC^{\top}\right)^{-1}CPA^{\top}.$
(7)

Equation (7) always has a unique solution $P$ under the assumption of
detectability of the pair $(A,C)$ [30].

Figure 1: System configuration.

####  III-B2 Residuals and Hypothesis Testing

Consider the process dynamics (4) and the steady-state Kalman filter (5), and
define the residual sequence $r_{k}:=y_{k}-C\hat{x}_{k}$, which can be shown
to evolve according to the difference equation:

$\displaystyle\left\\{\begin{aligned}
e_{k+1}&=(A-LC)e_{k}+t_{k}-Lw_{k}+(G-LH)\delta_{k},\\\
r_{k}&=Ce_{k}+w_{k}+H\delta_{k}.\end{aligned}\right.$ (8)

If there are no anomalies ($\delta_{k}=\mathbf{0}$), the mean and the
covariance matrix of the residual are as follows:

$\displaystyle\left\\{\begin{aligned}
&E[r_{k}]=CE[e_{k}]+E[w_{k}]=\mathbf{0}_{n_{y}\times 1},\\\\[5.69054pt]
&E[r_{k}r^{\top}_{k}]=CPC^{\top}+\Sigma^{w}=:\Sigma\in{\mathbb{R}^{{n_{y}}\times{n_{y}}}}.\end{aligned}\right.$
(9)

For this residual, we identify two hypotheses to be tested: $\mathcal{H}_{0}$
the normal mode (no anomalies) and $\mathcal{H}_{1}$ the faulty mode (with
anomalies). Under the normal mode, $\mathcal{H}_{0}$, the statistics of the
residual is given by (9). However, in the presence of anomalies
($\mathcal{H}_{1}$), we expect that the statistics of the residual are
different from those in the normal mode, i.e.,

$\mathcal{H}_{1}:\left\\{\begin{aligned} &E[r_{k}]\neq\mathbf{0}_{n_{y}\times
1},\text{ and/or}\\\\[5.69054pt]
&E[r_{k}r^{\top}_{k}]\neq\Sigma.\end{aligned}\right.$ (10)

There exist many well-known hypothesis testing techniques, that could be used
to examine the residual and subsequently detect anomalies/faults, for
instance, SPRT [31, 32], CUSUM [33, 34], generalized likelihood ratio testing
[35], compound scalar testing (CST) [15], etc. Each technique has its own
advantages and disadvantages depending on the scenario. In this manuscript, we
consider a particular case of CST, namely the so-called chi-squared change
detection procedure [36, 37].

####  III-B3 Distance Measure and Chi-Squared Procedure

The input to the chi-squared procedure is a distance measure
$z_{k}\in{\mathbb{R}}$, i.e., a measure of how deviated the estimator is from
the sensor measurements. We use the quadratic distance measure

$z_{k}=r^{\top}_{k}\Sigma^{-1}r_{k},$ (11)

where $r_{k}$ and $\Sigma$ are the residual sequence and its covariance matrix
introduced in (9). If there are no anomalies ($\delta_{k}=\mathbf{0}$), i.e.
$E[r_{k}]=\mathbf{0}$ and $E[r_{k}r^{\top}_{k}]=\Sigma$, it follows that

$\left\\{\begin{array}[]{ll}E[z_{k}]=\text{tr}[\Sigma^{-1}\Sigma]+E[r_{k}]^{\top}\Sigma^{-1}E[r_{k}]=n_{y},\\\\[5.69054pt]
\begin{aligned}
\text{var}[z_{k}]&=2\text{tr}[\Sigma^{-1}\Sigma\Sigma^{-1}\Sigma]+4E[r_{k}]^{\top}\Sigma^{-1}\Sigma\Sigma^{-1}E[r_{k}]\\\\[2.84526pt]
&=2n_{y},\end{aligned}\end{array}\right.$ (12)

see, e.g. [38] for details. Because $r_{k}\sim\mathcal{N}(\mathbf{0},\Sigma)$,
$z_{k}$ follows a chi-squared distribution with $n_{y}$ degrees of freedom.
The chi-squared procedure is characterized by $z_{k}$ and its cumulative
distribution as follows:
Chi-Squared procedure:
If $z_{k}=r_{k}^{\top}\Sigma^{-1}r_{k}>\alpha$, $\tilde{k}=k$.
Design parameter: threshold $\alpha\in\mathbb{R}_{+}$.
Output: alarm time(s) $\tilde{k}$.

The idea is that if $z_{k}$ exceeds the threshold $\alpha$, alarms are
triggered, and the time instant $k$ is identified as an alarm time (denoted by
$\tilde{k}$). The parameter $\alpha$ is selected to satisfy a desired false
alarm rate $\mathcal{A}^{*}$. Assume that there are no anomalies and consider
the chi-squared procedure with threshold $\alpha\in\mathbb{R^{+}}$,
$r_{k}\sim\mathcal{N}(\mathbf{0},\Sigma)$. Let
$\alpha=\alpha^{*}:=2{P}^{-1}\left(\frac{n_{y}}{2},1-\mathcal{A}^{*}\right)$,
where ${P}^{-1}(\cdot,\cdot)$ denotes the inverse regularised lower incomplete
gamma function, then the false alarm rate $\mathcal{A}$, induced by the
threshold $\alpha$, satisfies $\mathcal{A}=\mathcal{A}^{*}$ (see [37] for
proof and details).

###  III-C Privacy Mechanism and Adversarial Capabilities

We aim to prevent adversaries from accurately estimating the private output
$s_{k}$ (see (4)). To this end, we randomize measurements $y_{k}$ and input
signals $u_{k}$ before transmission and send the corrupted data to the remote
station instead. The idea is to randomize $y_{k}$ and $u_{k}$ through a
stochastic mapping $M(k,y_{k},u_{k})$ of the form:

$M(k,y_{k},u_{k}):=\left\\{\begin{array}[]{ll}\tilde{y}_{k}=y_{k}+v_{k},\\\\[5.69054pt]
\tilde{u}_{k}=u_{k}+j_{k},\end{array}\right.$ (13)

for some dependent Gaussian processes,
$v_{k}\sim\mathcal{N}[\mathbf{\mathbf{0}},\Sigma^{v}_{k}]$ and
$j_{k}\sim\mathcal{N}[\mathbf{\mathbf{0}},\Sigma^{j}_{k}]$ with zero mean and
time-varying covariance matrices $\Sigma^{v}_{k}$ and $\Sigma^{j}_{k}$ to be
designed for each time instant. The randomized vectors $\tilde{y}_{k}$ and
$\tilde{u}_{k}$ are disclosed to the remote station to be used by the remote
anomaly detector to detect faults/attacks and trigger alarms, see Figure 1.
Note that, from the remote station point of view (as they can only use the
received signals $(\tilde{y}_{k},\tilde{u}_{k})$), the system dynamics have
the following structure in the presence of anomalies:

$\displaystyle\left\\{\begin{aligned}
x_{k+1}&=Ax_{k}+B\tilde{u}_{k}+t_{k}+G\delta_{k}\\\
&=Ax_{k}+Bu_{k}+Bj_{k}+t_{k}+G\delta_{k},\\\
\tilde{y}_{k}&=y_{k}+v_{k}=Cx_{k}+w_{k}+H\delta_{k}+v_{k}.\end{aligned}\right.$
(14)

We consider worst-case adversaries that eavesdrop data at the communication
network and/or the remote station. They do not only have access to all
distorted sensor measurements $\tilde{y}_{k}$ and distorted input signals
$\tilde{u}_{k}$, but also have prior knowledge of the dynamics and the
stochastic properties of the system, i.e., matrices
$(A,B,C,D,\Sigma^{x}_{1},\mu^{x}_{1},\Sigma^{t},\Sigma^{w})$ are known by the
adversary. Moreover, the adversary also knows the means and covariance
matrices of $\tilde{y}_{k}$, $(\mu^{\tilde{y}}_{k},\Sigma^{\tilde{y}}_{k})$,
and $\tilde{u}_{k}$, $(\mu^{\tilde{u}}_{k},\Sigma^{\tilde{u}}_{k})$, as these
can be estimated from the disclosed data $(\tilde{y}_{k},\tilde{u}_{k})$. We
assume that the adversary uses a linear MMSE estimator (see Definition 1) to
reconstruct the private output $s_{k}$, which, for jointly Gaussian vectors,
produces the best estimation performance among all unbiased estimators [29].
In practice, actual adversaries would typically not have all the capabilities
that we assume here. However, if we maximize privacy under such worst-case
adversaries, we ensure that adversaries with less capabilities perform even
worse (or equal at most).

###  III-D Performance Degradation Due to Privacy Mechanisms

Note that whenever the privacy mechanism in (13) is in place, the anomaly
detector at the remote station uses the distorted disclosed data
$(\tilde{y}_{k},\tilde{u}_{k})$ instead of $(y_{k},u_{k})$ to drive its Kalman
filter and subsequently detect anomalies. That is, the remote filter takes the
following form:

$\hat{x}_{k+1}=A\hat{x}_{k}+B\tilde{u}_{k}+L(\tilde{y}_{k}-C\hat{x}_{k}).$
(15)

Because we assume that the gain matrix $L$ has been designed a priori without
considering the privacy-preserving mechanism (so considering the stochastic
properties of $(y_{k},u_{k})$ and not those of
$(\tilde{y}_{k},\tilde{u}_{k})$), the use of distorting mechanisms in (13)
will inevitably lead to anomaly detection performance degradation.

Consider the system (14) and the distorted Kalman filter (15), and define the
distorted residual sequence $\tilde{r}_{k}:=\tilde{y}_{k}-C\hat{x}_{k}$, which
can be shown to evolve according to the difference equation:

$\left\\{\begin{array}[]{ll}\begin{aligned}
\tilde{e}_{k+1}&=(A-LC)\tilde{e}_{k}+t_{k}-Bj_{k}-Lv_{k}-Lw_{k}\\\
&+(G-LH)\delta_{k},\end{aligned}\\\\[5.69054pt]
\tilde{r}_{k}=C\tilde{e}_{k}+w_{k}+v_{k}+H\delta_{k}.\end{array}\right.$ (16)

Hence, if we apply the privacy filter (13), the statistics of the distorted
residual $\tilde{r}_{k}$ in (16) in the anomaly-free case
($\delta_{k}=\mathbf{0}$) are:

$\left\\{\begin{array}[]{ll}E[\tilde{r}_{k}]=CE[\tilde{e}_{k}]+E[w_{k}]+E[v_{k}]=\mathbf{0}_{n_{y}\times
1},\\\\[5.69054pt] \begin{aligned}
E[\tilde{r}_{k}\tilde{r}^{\top}_{k}]&=CE[\tilde{e}_{k}\tilde{e}^{\top}_{k}]C^{\top}+\Sigma^{w}+\Sigma^{v}_{k}\\\\[5.69054pt]
&=C(P+B\Sigma^{j}_{k}B^{\top}+L\Sigma^{v}_{k}L^{\top})C^{\top}+\Sigma^{w}+\Sigma^{v}_{k}\\\\[5.69054pt]
&=\Sigma+\Sigma^{v}_{k}+CL\Sigma^{v}_{k}L^{\top}C^{\top}+CB\Sigma^{j}_{k}B^{\top}C^{\top}\\\\[5.69054pt]
&=:\tilde{\Sigma}_{k}.\end{aligned}\end{array}\right.$ (17)

Then, by applying the privacy filter to the measurement and input, the
residual becomes an i.i.d. multivariate Gaussian process with covariance
matrix $\tilde{\Sigma}_{k}$ as in (17). It follows that the distance measure
used for anomaly detection takes the form
$\tilde{z}_{k}=\tilde{r}^{\top}_{k}\Sigma^{-1}\tilde{r}_{k}$ in the distorted
case. Note that because we assume the detector is designed without considering
the privacy mechanism, the distorted distance measure $\tilde{z}_{k}$ is still
normalized by the undistorted residual covariance metric $\Sigma$ introduced
in (9). This results in $\tilde{z}_{k}$ not following a chi-squared
distribution anymore. Therefore, the detection threshold $\alpha$ (which was
computed for the undistorted measure $z_{k}$) will always lead to performance
degradation when applied to the distorted $\tilde{z}_{k}$ for detection.

##  IV Privacy-utility trade-off Formulated as an Optimization Problem

As discussed in previous sections, adding the privacy mechanism changes the
residual properties, see (17), which results in increasing the false alarm
rate. The privacy noise in (13) acts as an additional source of uncertainty
from the detection point of view, which may reduce the performance of the
anomaly detection scheme. We aim to formulate an optimization problem to
characterize a desired trade-off between privacy and the performance of the
anomaly detection algorithm. For a given time horizon
$\mathcal{K}:=\\{1,\ldots,K\\}$, $K\in\mathbb{N}$, the aim of the proposed
privacy scheme is to make inference of the sequence of private output vectors,
${s^{K}}=\operatorname{col}(s_{1},\ldots,s_{K})$, from the distorted disclosed
sequences,
${{\tilde{y}}^{K}}=\operatorname{col}({\tilde{y}}_{1},\ldots,{\tilde{y}}_{K})$
and
${{\tilde{u}}^{K}}=\operatorname{col}({\tilde{u}}_{1},\ldots,{\tilde{u}}_{K})$,
by designing privacy noises sequences,
${v^{K}}=\operatorname{col}(v_{1},\ldots,v_{K})$ and
${j^{K}}=\operatorname{col}(j_{1},\ldots,j_{K})$, as difficult as possible
without degrading the anomaly detection performance excessively.
As privacy metric, an intuitive candidate to use is the mutual information
between private and disclosed data, i.e.,
$I[s^{K};{\tilde{y}}^{K},{\tilde{u}}^{K}]$. However, because the input
sequence $u^{K}$ is deterministic and $v^{K}$ and $j^{K}$ are independent, it
is easy to verify that
$I[s^{K};\tilde{y}^{K},\tilde{u}^{K}]=I[s^{K};\tilde{y}^{K}]$ (see[16] for
proof and details). That is, in the proposed setting, the randomized input
data $\tilde{u}^{K}$ does not affect $I[s^{K};\tilde{y}^{K},\tilde{u}^{K}]$ at
all. To overcome this obstacle, we add to $I[s^{K};\tilde{y}^{K}]$ the
negative differential entropy of the difference of $u^{K}$ and
$\tilde{u}^{K}$, $-h[u^{K}-\tilde{u}^{K}]=-h[j^{K}]$, to capture the
uncertainty between original and disclosed input data. That is, we propose
$I[s^{K};\tilde{y}^{K}]-h[j^{K}]$ as privacy metric to be minimized.

To capture performance degradation in the anomaly detection algorithm, we
consider the change in the original false alarm rate $\mathcal{A}^{*}$ induced
by the privacy mechanism. Let $\alpha$ denote the threshold designed in the
undistorted case to ensure a false alarm rate of $\mathcal{A}^{*}$ for the
chi-squared detector. Then, we seek privacy mechanisms that do not change the
false alarm rate by more than an allowed $\epsilon$ when the original $\alpha$
and the distorted residuals $\tilde{r}_{k}$ are used to run the detector,
i.e., we seek to satisfy the following probabilistic constraint:

$\displaystyle\operatorname{Pr}[\tilde{z}_{k}>\alpha]<\mathcal{A}^{*}+\epsilon,$
(18)

which is equivalent to

$\displaystyle\operatorname{Pr}[\tilde{z}_{k}\leq\alpha]>1-\mathcal{A}^{*}-\epsilon,$
(19)

where $\epsilon$ is the false alarm distortion level. Using the definition of
Cumulative Distribution Function (CDF), we can write (19) as

$F_{\tilde{z}_{k}}(\alpha)>1-\mathcal{A}^{*}-\epsilon,$ (20)

where $F_{\tilde{z}_{k}}(\alpha)$ is the CDF of $\tilde{z}_{k}$, i.e.,
$F_{\tilde{z}_{k}}(\alpha)=\operatorname{Pr}[\tilde{z}_{k}\leq\alpha]$.
Therefore, the optimization problem for maximizing privacy while considering
the above constraint on the false alarm rate can be written as follows.

Problem 1 Given the system dynamics (4), $k\in\mathcal{K}$ with
$\mathcal{K}=\\{1,\ldots,K\\}$, $K\in\mathbb{N}$, desired false alarm rate
$\mathcal{A}^{*}$, maximum distortion level ${\epsilon}\in{\mathbb{R}^{+}}$,
the chi-squared procedure threshold $\alpha$, the distorting mechanism (13),
and the remote Kalman filter (15), find the sequence of covariance matrices
${\Sigma^{v}_{k}}$ and ${\Sigma^{j}_{k}}$ in (13) as the solution to the
following optimization problem:

$\left\\{\begin{aligned}
&\min_{{\Sigma^{v}_{k}},{\Sigma^{j}_{k}},k\in\mathcal{K}}\
I[{s^{K}};{\tilde{y}^{K}}]-h[{j}^{K}],\\\\[2.84526pt] &\hskip
22.76219pt\text{s.t.
}F_{\tilde{z}_{k}}(\alpha)>1-\mathcal{A}^{*}-\epsilon,\\\\[2.84526pt] &\hskip
22.76219pt\text{and }(v^{K},j^{K},y^{K})\text{ mutually
independent}.\end{aligned}\right.$ (21)

##  V Solution to Problem 1

To solve the optimization problem (21), we first need to write the cost
function and constraints in terms of our design variables $\Sigma^{v}_{k}$ and
$\Sigma^{j}_{k}$, $k\in\mathcal{K}$.

###  V-A Cost Function: Formulation and Convexity

We first formulate the cost function, $I[{s^{K}};{\tilde{y}^{K}}]-h[{j}^{K}]$,
in terms of the design variables $\Sigma^{v}_{k}$ and $\Sigma^{j}_{k}$. The
differential entropy $h[{j}^{K}]$ is fully characterized by the covariance
matrix of ${j}^{K}$, $\Sigma^{j}_{K}:=E[j^{K}{j^{K}}^{\top}]$. By Definition
$2$, $h[{j}^{K}]$ is given by

$h[{j}^{K}]=\frac{1}{2}\log\det\left(\Sigma^{j}_{K}\right)+\frac{Kn_{u}}{2}+\frac{Kn_{u}}{2}\log(2\pi).$
(22)

The mutual information $I\left[{{s^{K}};{\tilde{y}^{K}}}\right]$ can be
written in terms of differential entropies as
$I[{s^{K}};{\tilde{y}^{K}}]=h[s^{K}]+h[\tilde{y}^{K}]-h[s^{K},\tilde{y}^{K}]$,
see Definition 3. Moreover, these entropies are fully characterized by the
covariance matrices of the corresponding random vectors, see Definition 2. So,
to write $I\left[{{s^{K}};{\tilde{y}^{K}}}\right]$ in terms of
$\Sigma^{v}_{k}$, we need to write the covariance matrices of $s^{K}$ and
$\tilde{y}^{K}$, and their joint covariance in terms of $\Sigma^{v}_{k}$. By
lifting the system dynamics (4) over $\mathcal{K}$, the stacked vector
$((\tilde{y}^{K})^{\top},(s^{K})^{\top})^{\top}\in\mathbb{R}^{K(n_{y}+n_{s})}$
can be written as

$\displaystyle\begin{bmatrix}\tilde{y}^{K}\\\
s^{K}\end{bmatrix}=\begin{bmatrix}\tilde{C}_{K}\\\
\tilde{D}_{K}\end{bmatrix}F_{K}x_{1}+\begin{bmatrix}\tilde{C}_{K}\\\
\tilde{D}_{K}\end{bmatrix}J_{K}t^{K-1}$ (23) $\displaystyle\hskip
42.67912pt+\begin{bmatrix}\tilde{C}_{K}\\\
\tilde{D}_{K}\end{bmatrix}N_{K}u^{K-1}+\begin{bmatrix}I\\\
\mathbf{0}\end{bmatrix}(w^{K}+v^{K}),$

with sequence of inputs
${u^{K-1}}=(u_{1}^{\top},\ldots,u_{K-1}^{\top})^{\top}$, stacked matrices
$N_{K}:=J_{K}(I_{K-1}\otimes B)$, $\tilde{C}_{K}:=I_{K}\otimes C$,
$\tilde{D}_{K}:=I_{K}\otimes D$, and

$\left\\{\begin{aligned}
F_{K}&:=\begin{bmatrix}I&A^{\top}&\ldots&(A^{\top})^{K-1}\end{bmatrix}^{\top},\\\
J_{K}&:=\begin{bmatrix}\mathbf{0}&\mathbf{0}&\mathbf{0}&\cdots&\mathbf{0}\\\
I&\mathbf{0}&\mathbf{0}&\cdots&\mathbf{0}\\\
A&I&\mathbf{0}&\cdots&\mathbf{0}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\\[2.84526pt]
(A)^{K-2}&(A)^{K-3}&(A)^{K-4}&\cdots&I\end{bmatrix}.\end{aligned}\right.$ (24)

Let ${{\Sigma}^{v}_{K}}\in\mathbb{R}^{Kn_{y}\times Kn_{y}}$ and
${{\Sigma}^{j}_{K}}\in\mathbb{R}^{Kn_{u}\times Kn_{u}}$ denote the non-
diagonal covariance matrix of the stacked additive dependent vectors $v^{K}$
and $j^{K}$. In the following lemma, we give a closed-form expression of the
joint density of $(\tilde{y}^{K},s^{K})$.

######  Lemma 1.

$\begin{pmatrix}\tilde{y}^{K}\\\
s^{K}\end{pmatrix}\sim\mathcal{N}\left[\mu^{\tilde{y},s}_{K},\Sigma^{\tilde{y},s}_{K}\right],$

with mean $\mu^{\tilde{y},s}_{K}\in\mathbb{R}^{K(n_{s}+n_{y})}$ and covariance
matrix $\Sigma^{\tilde{y},s}_{K}\in\mathbb{R}^{K(n_{s}+n_{y})\times
K(n_{s}+n_{y})}$, $\Sigma^{\tilde{y},s}_{K}>0$_:_

$\displaystyle\mu^{\tilde{y},s}_{K}$
$\displaystyle:=\begin{bmatrix}\tilde{C}_{K}\\\
\tilde{D}_{K}\end{bmatrix}F_{K}\mu^{x}_{1}+\begin{bmatrix}\tilde{C}_{K}\\\
\tilde{D}_{K}\end{bmatrix}N_{K}u^{K-1},$ (25)
$\displaystyle\Sigma^{\tilde{y},s}_{K}$ $\displaystyle:=\begin{bmatrix}I\\\
\mathbf{0}\end{bmatrix}(I_{K}\otimes\Sigma^{w})\begin{bmatrix}I\\\
\mathbf{0}\end{bmatrix}^{\top}+\begin{bmatrix}I\\\
\mathbf{0}\end{bmatrix}{{\Sigma}^{v}_{K}}\begin{bmatrix}I\\\
\mathbf{0}\end{bmatrix}^{\top}$ (26)
$\displaystyle+\begin{bmatrix}\tilde{C}_{K}\\\
\tilde{D}_{K}\end{bmatrix}F_{K}\Sigma^{x}_{1}F_{K}^{\top}\begin{bmatrix}\tilde{C}_{K}\\\
\tilde{D}_{K}\end{bmatrix}^{\top}$
$\displaystyle+\begin{bmatrix}\tilde{C}_{K}\\\
\tilde{D}_{K}\end{bmatrix}J_{K}(I_{K-1}\otimes\Sigma^{t})J_{K}^{\top}\begin{bmatrix}\tilde{C}_{K}\\\
\tilde{D}_{K}\end{bmatrix}^{\top}.$

_Proof_ : To simplify notation, we introduce the stacked vector
${\Theta^{K}}:={({{{({{\tilde{y}^{K}}})}^{\top}},{{({{s^{K}}})}^{\top}}})^{\top}}$.
By assumption, the initial condition $x_{1}$, and the processes, $t_{k}$,
$w_{k}$, and $v_{k}$, $k\in\mathbb{N}$, are mutually independent, and
$x_{1}\sim\mathcal{N}[{\mu_{1}^{x},\Sigma_{1}^{x}}]$,
$t_{k}\sim\mathcal{N}[{\mathbf{0},\Sigma^{t}}]$,
$w_{k}\sim\mathcal{N}[{\mathbf{0},\Sigma^{w}}]$,
$v_{k}\sim\mathcal{N}[{\mathbf{0},\Sigma^{v}}]$, for some positive definite
covariance matrices $\Sigma_{1}^{x}$, $\Sigma^{t}$, $\Sigma^{w}$, and
$\Sigma^{v}$. Then, see [39] for details, we have
${L_{1}}x\left(1\right)\sim\mathcal{N}[{{L_{1}}\mu_{1}^{x},{L_{1}}\Sigma_{1}^{x}L_{1}^{\top}}]$,
${L_{2}}t^{K-1}\sim\mathcal{N}[\mathbf{0},{{L_{2}}({I_{K-1}}\otimes{\Sigma^{t}}){L_{2}}^{T}}]$,
${L_{3}}w^{K}\sim\mathcal{N}[\mathbf{0},{{L_{3}}({I_{K}}\otimes{\Sigma^{w}}){L_{3}}^{\top}}]$,
${L_{4}}v^{K}\sim\mathcal{N}[\mathbf{0},{{L_{4}}{{\Sigma}^{v}_{K}}{L_{4}}^{\top}}]$,
for any deterministic matrices $L_{j}$, $j=1,2,3,4$, of appropriate
dimensions. It follows that the stacked vector ${\Theta^{K}}$ given in (23) is
the sum of a deterministic vector,
${({\tilde{C}_{K}^{\top},\tilde{D}_{K}^{T}})^{\top}}{N_{K}}{u^{K-1}}$, and
four independent normally distributed vectors. Therefore, ${\Theta^{K}}$
follows a multivariate normal distribution with
$E[{{\Theta^{K}}}]=\mu^{\tilde{y},s}_{K}$ as in (25). By inspection, using the
expression of ${\Theta^{K}}$ in (23), mutual independence among $x_{1}$,
$t_{k}$, $w_{k}$, and $v_{k}$, $k\in N$, and the definition of
$\Sigma_{1}^{x}$,
$\Sigma_{1}^{x}=E[{\left({x_{1}-\mu_{1}^{x}}\right){{\left({x_{1}-\mu_{1}^{x}}\right)}^{\top}}}]$,
it can be verified that the covariance matrix of ${\Theta^{K}}$,
$E[{({{\Theta^{K}}-E({{\Theta^{K}}})}){{({{\Theta^{K}}-E({{\Theta^{K}}})})}^{\top}}}]$,
is given by $\Sigma^{\tilde{y},s}_{K}$ in (26). It remains to prove that the
distribution of ${\Theta^{K}}$ is not degenerate, that is,
$\Sigma^{\tilde{y},s}_{K}>0$. Note that $\Sigma^{\tilde{y},s}_{K}$ in (26) can
be written as

$\displaystyle\Sigma_{K}^{\tilde{y},s}$
$\displaystyle=\left[{\begin{array}[]{*{20}{c}}\Sigma_{K}^{\tilde{y}}&{\tilde{C}_{K}Q\tilde{D}_{K}^{\top}}\\\
{\tilde{D}_{K}Q\tilde{C}_{K}^{\top}}&{\tilde{D}_{K}Q\tilde{D}_{K}^{\top}}\end{array}}\right],$
(29) $\displaystyle\Sigma_{K}^{\tilde{y}}$
$\displaystyle={\tilde{C}_{K}Q\tilde{C}_{K}^{\top}+(I_{K}\otimes\Sigma^{w})+{{\Sigma}^{v}_{K}}},$
(30)

with
$Q={F_{K}}\Sigma_{1}^{x}F_{K}^{\top}+{J_{K}}({{I_{K-1}}\otimes{\Sigma^{t}}})J_{K}^{\top}$.
Necessary and sufficient conditions for $\Sigma_{K}^{\tilde{y},s}>0$ are its
right-lower block is positive definite
(${\tilde{D}_{K}Q\tilde{D}_{K}^{\top}}>0$) and that the Schur complement of
block ${\tilde{D}_{K}Q\tilde{D}_{K}^{\top}}$ of $\Sigma_{K}^{\tilde{y},s}$,
denoted as $\Sigma_{K}^{\tilde{y},s}/({\tilde{D}_{K}Q\tilde{D}_{K}^{\top}})$,
is positive definite ([40], Theorem 1.12). The right-lower block is positive
definite if $\tilde{D}_{K}$ is full row rank and $Q$ is positive definite.
Because $D$ is full row rank by assumption, matrix
$\tilde{D}_{K}=\left({{I_{K}}\otimes D}\right)$ is also full row rank ([41],
Theorem 4.2.15). Note that $Q$ can be factored as follows

$\displaystyle
Q=\left[{\begin{array}[]{*{20}{c}}{{F_{K}}}&{{J_{K}}}\end{array}}\right]\underbrace{\left[{\begin{array}[]{*{20}{c}}{\Sigma_{1}^{x}}&\mathbf{0}\\\
\mathbf{0}&{{I_{K-1}}\otimes{\Sigma^{t}}}\end{array}}\right]}_{Q^{\prime}}\left[{\begin{array}[]{*{20}{c}}{{F_{K}}}\\\
{{J_{K}}}\end{array}}\right].$ (36)

That is, $Q$ is a linear transformation of the block diagonal matrix
${Q^{\prime}}$ above. By inspection, it can be verified that matrix
$P=[F_{K}\hskip 2.84526ptJ_{K}]$, see (24), is lower triangular with identity
matrices on the diagonal; thus, $P$ is invertible. It follows that
$Q=PQ^{\prime}P^{\top}$ is a congruence transformation of ${Q^{\prime}}$ [42].
Hence, $Q$ is positive definite if and only if the block diagonal matrices of
${Q^{\prime}}$ are positive definite [42]. Matrices ${\Sigma_{1}^{x}}$ and
${\Sigma^{t}}$ are positive definite by assumption (which implies
${{I_{K-1}}\otimes{\Sigma^{t}}}>0$), and thus we can conclude that $Q>0$,
which implies ${\tilde{D}_{K}Q\tilde{D}_{K}^{\top}}>0$, because
$\tilde{D}_{K}$ is full row rank. Then, we need to prove that the Schur
complement $\Sigma_{K}^{\tilde{y},s}/({\tilde{D}_{K}Q\tilde{D}_{K}^{\top}})$,
is positive definite. This Schur complement is given by

$\displaystyle\Sigma_{K}^{\tilde{y},s}/({\tilde{D}_{K}Q\tilde{D}_{K}^{\top}})=I_{K}\otimes\Sigma^{w}+{{\Sigma}^{v}_{K}}+$
(37)
$\displaystyle\tilde{C}_{K}\left({Q-Q\tilde{D}_{K}^{\top}\left({\tilde{D}_{K}Q\tilde{D}_{K}^{\top}}\right)^{-1}\tilde{D}_{K}Q}\right)\tilde{C}_{K}^{\top}.$

Since $(I_{K}\otimes\Sigma^{w})+{{\Sigma}^{v}_{K}}$ is positive definite, a
sufficient condition for
$\Sigma_{K}^{\tilde{y},s}/({\tilde{D}_{K}Q\tilde{D}_{K}^{\top}})>\mathbf{0}$
is

$\displaystyle
Q^{\prime\prime}:=\left({Q-Q\tilde{D}_{K}^{\top}\left({\tilde{D}_{K}Q\tilde{D}_{K}^{\top}}\right)^{-1}\tilde{D}_{K}Q}\right)>\mathbf{0}.$
(38)

Regarding ${Q^{\prime\prime}}$ as the Schur complement of a higher dimensional
matrix ${Q^{\prime\prime\prime}}$, we can conclude that:

$\displaystyle\begin{array}[]{l}Q^{\prime\prime}>\mathbf{0}\Leftrightarrow
Q^{\prime\prime\prime}=\left[{\begin{array}[]{*{20}{c}}Q&{Q\tilde{D}_{K}^{\top}}\\\
{\tilde{D}_{K}Q}&{\tilde{D}_{K}Q\tilde{D}_{K}^{\top}}\end{array}}\right]\\\
=\left[{\begin{array}[]{*{20}{c}}Q\\\
{\tilde{D}_{K}Q}\end{array}}\right]{Q^{-1}}\left[{\begin{array}[]{*{20}{c}}Q&{Q\tilde{D}_{K}^{\top}}\end{array}}\right]>\mathbf{0},\end{array}$
(46)

which is trivially true because $Q^{-1}$ is positive definite since
$Q>\mathbf{0}$. Hence, ${\tilde{D}_{K}Q\tilde{D}_{K}^{\top}}$ and
$\Sigma_{K}^{\tilde{y},s}/({\tilde{D}_{K}Q\tilde{D}_{K}^{\top}})$ are both
positive definite, and thus $\Sigma_{K}^{\tilde{y},s}>\mathbf{0}$.
$\blacksquare$
To obtain the densities of $\tilde{y}^{K}$ and $s^{K}$ and their cross-
covariance $\Sigma^{s\tilde{y}}_{K}$, we just marginalize (see [39]) their
joint density in Lemma $1$ over $s_{K}$ and $\tilde{y}_{K}$, see the following
corollary.

######  Corollary 1.

$\tilde{y}^{K}\sim\mathcal{N}[\mu^{\tilde{y}}_{K},\Sigma^{\tilde{y}}_{K}]$ and
$s^{K}\sim\mathcal{N}[\mu^{s}_{K},\Sigma^{s}_{K}]$_:_

$\displaystyle\mu^{\tilde{y}}_{K}=\tilde{C}_{K}F_{K}\mu^{x}_{1}+\tilde{C}_{K}N_{K}u^{K-1},$
$\displaystyle\mu^{s}_{K}=\tilde{D}_{K}F_{K}\mu^{x}_{1}+\tilde{D}_{K}N_{K}u^{K-1},$
$\displaystyle\Sigma^{\tilde{y}}_{K}=\begin{pmatrix}I_{n_{y}}\\\
\mathbf{0}\end{pmatrix}^{\top}\Sigma^{\tilde{y},s}_{K}\begin{pmatrix}I_{n_{y}}\\\
\mathbf{0}\end{pmatrix}\in\mathbb{R}^{Kn_{y}\times Kn_{y}},$
$\displaystyle\Sigma^{s}_{K}=\begin{pmatrix}\mathbf{0}\\\
I_{n_{s}}\end{pmatrix}^{\top}\Sigma^{\tilde{y},s}_{K}\begin{pmatrix}\mathbf{0}\\\
I_{n_{s}}\end{pmatrix}\in\mathbb{R}^{Kn_{s}\times Kn_{s}},$
$\displaystyle\Sigma^{s\tilde{y}}_{K}=\begin{pmatrix}\mathbf{0}\\\
I_{n_{s}}\end{pmatrix}^{\top}\Sigma^{\tilde{y},s}_{K}\begin{pmatrix}I_{n_{y}}\\\
\mathbf{0}\end{pmatrix}\in\mathbb{R}^{Kn_{s}\times Kn_{y}}.$ (47)

At this point, we have the three covariance matrices
$(\Sigma^{\tilde{y}}_{K},\Sigma^{s}_{K},\Sigma^{s\tilde{y}}_{K})$ that we need
to compute the mutual information $I\left[{{s^{K}};{\tilde{y}^{K}}}\right]$.
Note that $\Sigma^{s}_{K}$ and $\Sigma_{K}^{s\tilde{y}}$ are independent of
the design variables $\Sigma_{K}^{v}$ and $\Sigma_{K}^{j}$, and
$\Sigma_{K}^{\tilde{y}}$ depends on $\Sigma_{K}^{v}$:

$\displaystyle\Sigma_{K}^{\tilde{y},s}=\begin{pmatrix}\Sigma_{K}^{\tilde{y}}&{\Sigma_{K}^{s\tilde{y}}}^{\top}\\\
\Sigma_{K}^{s\tilde{y}}&\Sigma_{K}^{s}\end{pmatrix},$ (48)
$\displaystyle\Sigma_{K}^{\tilde{y}}={\tilde{C}_{K}Q\tilde{C}_{K}^{\top}+(I_{K}\otimes\Sigma^{w})+{{\Sigma}^{v}_{K}}},$
(49) $\displaystyle\Sigma^{s}_{K}={\tilde{D}_{K}Q\tilde{D}_{K}^{\top}},$ (50)
$\displaystyle\Sigma_{K}^{s\tilde{y}}={\tilde{D}_{K}Q\tilde{C}_{K}^{\top}},$
(51)

with $Q$ as in (36) independent of $\Sigma^{v}_{K}$. Then, given
$(\Sigma_{K}^{\tilde{y},s},\Sigma^{\tilde{y}}_{K},\Sigma^{s}_{K},\Sigma^{s\tilde{y}}_{K})$
in (48)-(51), by Definition 2, we can write $I[{s^{K}};{\tilde{y}^{K}}]$ as
follows:

$\displaystyle
I[{s^{K}};{\tilde{y}^{K}}]=h\left[s^{K}\right]+h\left[\tilde{y}^{K}\right]-h[\tilde{y}^{K},s^{K}]$
(52a)
$\displaystyle=\frac{1}{2}\log{\det\left({\Sigma_{K}^{s}}\right)}-\frac{1}{2}\log{\det{\left({\Sigma}_{K}^{s}-{{\Sigma}_{K}^{s\tilde{y}}}{{\Sigma}_{K}^{\tilde{y}}}^{-1}{{\Sigma}_{K}^{s\tilde{y}}}^{\top}\right),}}$
(52b)

where (52b) follows from standard determinant and logarithm formulas. In the
following lemma, we prove that minimizing the cost function
$I[s^{K};\tilde{y}^{K}]-h[j^{K}]$ using
$({{\Sigma}^{v}_{K}},{{\Sigma}^{j}_{K}})$ as optimization variables is
equivalent to solving a convex program subject to Linear Matrix Inequalities
(LMI) constraints.

######  Lemma 2.

Minimizing $I[s^{K};\tilde{y}^{K}]-h[j^{K}]$ is equivalent to solving the
following convex program:

$\displaystyle\left\\{\begin{aligned}
&\min_{\Sigma^{j}_{K},{\Pi_{K}},\Sigma_{K}^{v}}\
-\log{\det\left({\Pi_{K}}\right)}-\log\det\left(\Sigma^{j}_{K}\right)\\\\[2.84526pt]
&\hskip 11.38109pt\text{\emph{s.t.
}}\begin{bmatrix}\Sigma_{K}^{s}-\Pi_{K}&\Sigma_{K}^{s\tilde{y}}\\\\[5.69054pt]
{\Sigma_{K}^{s\tilde{y}}}^{\top}&\Sigma_{K}^{\tilde{y}}\end{bmatrix}\geq\mathbf{0},\Pi_{K}>\mathbf{0}.\end{aligned}\right.$
(53)

_Proof_ : (a) For any positive definite matrix $\Sigma$, the function
$-\log\det(\Sigma)$ is convex in $\Sigma$ [43]. It follows that $-h[j^{K}]$ in
(22) is a convex function of $\Sigma^{j}_{K}$. Minimizing $-h[j^{K}]$ amounts
to minimizing $-\log\det\left(\Sigma^{j}_{K}\right)$ as all other terms in
(22) are constants.
(b) Next, consider the expression for $I[\tilde{y}^{K};s^{K}]$ in (52b). Due
to the monotonicity of the determinant function and the fact that
$\Sigma_{K}^{s}$ is independent of the design variables, minimizing (52b) is
equivalent to

$\displaystyle\left\\{\begin{aligned} &\min_{{\Pi_{K}},\Sigma_{K}^{v}}\
-\log{\det\left({\Pi_{K}}\right)}\\\\[2.84526pt] &\hskip 11.38109pt\text{s.t.
}\mathbf{0}<{\Pi_{K}}\leq{\left({\Sigma}_{K}^{s}-{{\Sigma}_{K}^{s\tilde{y}}}{{\Sigma}_{K}^{\tilde{y}}}^{-1}{{\Sigma}_{K}^{s\tilde{y}}}^{\top}\right).}\end{aligned}\right.$
(54)

(c) The inequality term in (54) can be rewritten using Schur complement
properties [40] as

$\displaystyle\left[{\begin{array}[]{*{20}{c}}{\Sigma_{K}^{s}}-{\Pi_{K}}&{{\Sigma}_{K}^{s\tilde{y}}}\\\
{{\Sigma}_{K}^{s\tilde{y}}}^{\top}&{{\Sigma}_{K}^{\tilde{y}}}\end{array}}\right]\geq\mathbf{0}\,{,{\Pi_{K}}>\mathbf{0}}.$
(57)

Since $\Sigma_{K}^{s\tilde{y}}$ is constant in the design variables (see
(51)), and $\Sigma_{K}^{\tilde{y}}$ is convex in $\Sigma_{K}^{v}$, we can
conclude that minimizing $I[s^{K};\tilde{y}^{K}]-h[j^{K}]$ is equivalent to
solving the convex program in (53). $\blacksquare$

By Lemma 2, minimizing the cost in (21) is equivalent to solving the convex
program in (53). Then, if the distortion constraint in (21) is convex in the
design variables, we can find optimal distorting mechanisms efficiently using
off-the-shelf optimization algorithms.

###  V-B Detection Performance Constraints: Formulation and Convexity

We consider a constraint on the false alarm rate in each time step, see (21),
as a metric of performance degradation of the detector. This constraint is
fully characterized by the CDF of the distorted distance measure
$\tilde{z}_{k}$. Unfortunately, no closed-form expression is available for
this CDF. To address this, we seek a lower bound on the CDF of $\tilde{z}_{k}$
to use it to enforce constraint (20). Define a standard Gaussian vector
$m_{k}\sim\mathcal{N}(\mathbf{0},I_{n_{y}})$. Since
$\tilde{r}_{k}\sim\mathcal{N}(\mathbf{0},\tilde{\Sigma}_{k})$, it is easy to
verify that $\tilde{r}_{k}$ can be written in terms of $m_{k}$ as follows:

$\tilde{r}_{k}=\left(\tilde{\Sigma}_{k}\right)^{1/2}m_{k},$ (58)

where $\left(\tilde{\Sigma}_{k}\right)^{1/2}$ is the square root of
$\tilde{\Sigma}_{k}$. Hence, $\tilde{z}_{k}$ can be written as:

$\displaystyle\tilde{z}_{k}$
$\displaystyle=\tilde{r}^{\top}_{k}\Sigma^{-1}\tilde{r}_{k}=m_{k}^{\top}(\tilde{\Sigma}^{1/2}_{k})^{\top}\Sigma^{-1}\tilde{\Sigma}^{1/2}_{k}m_{k}.$
(59)

Define the sequence $q_{k}:=\beta m_{k}^{\top}m_{k}$,
$\beta\in\mathbb{R}^{+}$, $k\in\mathcal{K}$, and let
$\Sigma^{\prime}_{k}:=\left((\tilde{\Sigma}_{k})^{1/2}\right)^{\top}\Sigma^{-1}(\tilde{\Sigma}_{k})^{1/2}$;
then, $\tilde{z}_{k}\leq q_{k}$ for all $k\in\mathcal{K}$ implies
$m_{k}^{\top}\Sigma^{\prime}_{k}m_{k}\leq\beta m_{k}^{\top}m_{k}$. It follows
that $\tilde{z}_{k}\leq q_{k}$ if and only if

$\displaystyle\Sigma^{\prime}_{k}\leq\beta{I}_{n_{y}}.$ (60)

To be able to find the constraints of Problem $1$ in terms of the CDF of the
upper bound $q_{k}$, we need to find a relation between the CDFs of
$\tilde{z}_{k}$ and $q_{k}$. In the next lemma, we show such a relation.

######  Lemma 3.

$\tilde{z}_{k}\leq q_{k}$, $k\in\mathcal{K}$, implies _:_

$F_{\tilde{z}_{k}}(\alpha)\geq F_{q_{k}}(\alpha),$ (61)

where $F_{\tilde{z}_{k}}(\alpha)$ and $F_{q_{k}}(\alpha)$ denote the CDFs of
$\tilde{z}_{k}$ and $q_{k}$, respectively.

_Proof_ : The proof is given in Appendix A. $\blacksquare$

From (60), Lemma 3, and (20), we can conclude that the constraints of the
optimization problem (21) hold if $\Sigma^{\prime}_{k}\leq\beta{I}_{n_{y}}$
and

$F_{\beta m_{k}^{\top}m_{k}}(\alpha)>1-\mathcal{A}^{*}-\epsilon\,\,\text{for
all }\,k\in\mathcal{K}.$ (62)

Therefore, if we choose $\beta$ to satisfy (62), the false alarm constraint in
(21) translate to $\Sigma^{\prime}_{k}\leq\beta{I}_{n_{y}}$. Given a desired
distortion level $\epsilon$ and desired false alarm rate $\mathcal{A}^{*}$, a
lower bound on $F_{\beta m_{k}^{\top}m_{k}}(\alpha)$ can be computed. Since
$m_{k}^{\top}m_{k}$ follows a chi-squared distribution with $n_{y}$ degree of
freedom, the CDF $F_{q_{k}}(\alpha)$ can be computed as follows (see [44]):

$F_{q_{k}}(\alpha)=P\left({n_{y}}/{2},{\alpha}/{2\beta}\right),$ (63)

where $P(\cdot,\cdot)$ is the regularized gamma function [44]. Substituting
(63) into (62) yields:

$P\left({n_{y}}/{2},{\alpha}/{2\beta}\right)>1-\mathcal{A}^{*}-\epsilon.$ (64)

Because the regularized gamma function is an increasing function of
$\frac{\alpha}{2\beta}$ [45], the largest $\beta$ satisfying (64) for given
$\alpha$, $\epsilon$, and $\mathcal{A}^{*}$ can be calculated using the
inverse of the lower incomplete gamma function [45] as follows:

$\beta<\beta^{*}:=\frac{\alpha}{2P^{-1}\left(\frac{n_{y}}{2},1-\mathcal{A}^{*}-\epsilon\right)}.$
(65)

Then, once we have $\beta^{*}$, the false alarm constraint in Problem $1$
amounts to $\Sigma^{\prime}_{k}\leq\beta{I}_{n_{y}}<\beta^{*}{I}_{n_{y}}$,
which is a characterization in terms of the design variables (as
$\Sigma^{\prime}_{k}$ depends on $\Sigma^{v}_{k}$ and $\Sigma^{j}_{k}$). We
state this explicitly in the following lemma.

######  Lemma 4.

Inequality
$\Sigma^{\prime}_{k}=(\tilde{\Sigma}^{1/2}_{k})^{\top}\Sigma^{-1}\tilde{\Sigma}^{1/2}_{k}<\beta^{*}{I}_{n_{y}}$
can be written in terms of $\Sigma^{v}_{k}$ and $\Sigma^{j}_{k}$ as:

$\tilde{\Sigma}_{k}<\beta^{*}\Sigma,$ (66)

with $\beta^{*}$ defined in (65) and
$\tilde{\Sigma}_{k}=\Sigma+\Sigma^{v}_{k}+CL\Sigma^{v}_{k}L^{\top}C^{\top}+CB\Sigma^{j}_{k}B^{\top}C^{\top}$.
_Proof_ : The proof is given in Appendix B. $\blacksquare$

By Lemma 2 and Lemma 4 the cost $I[s^{K};\tilde{y}^{K}]-h[j^{K}]$ and false
alarm constraint $F_{\tilde{z}_{k}}(\alpha)>1-\mathcal{A}^{*}-\epsilon$ are
written linearly in terms of the optimization variables ${{\Sigma}^{v}_{K}}$
and ${{\Sigma}^{j}_{K}}$. In what follows, we pose the complete nonlinear
convex program to solve Problem $1$.

######  Theorem 1.

Consider the system dynamics (4), distorting mechanism (13), $k\in\mathcal{K}$
with $\mathcal{K}=\\{1,\ldots,K\\}$ and $K\in\mathbb{N}$, desired false alarm
rate $\mathcal{A}^{*}$, maximum distortion level
${\epsilon}\in{\mathbb{R}^{+}}$, the chi-squared procedure threshold $\alpha$,
the anomaly detector Kalman filter (15), and covariance matrices
$\Sigma^{\tilde{y}}_{K}$, $\Sigma_{K}^{s}$, and $\Sigma_{K}^{s\tilde{y}}$ _(_
given in (48)-(51)_)_. Then, the covariance matrices ${{\Sigma}^{v}_{K}}>0$
and ${{\Sigma}^{j}_{K}}>0$ that minimize $I[s^{K};\tilde{y}^{K}]-h[j^{K}]$
subject to $F_{\tilde{z}_{k}}(\alpha)>1-\mathcal{A}^{*}-\epsilon$ can be found
by solving the convex program in (67).

_Proof:_ The expressions for the cost and constraints and their convexity
follow from Lemma 2 and Lemma 4. $\blacksquare$

$\displaystyle\min_{{\Pi_{K}},\Sigma_{K}^{v},{{\Sigma}^{j}_{K}}}-\log\det[\Sigma^{j}_{k}]-\log\det[\Pi_{K}],$
(67) $\displaystyle\text{ s.t. }\left\\{\begin{aligned}
&{\Pi_{K}}>\mathbf{0},\\\
&\begin{bmatrix}\Sigma_{K}^{s}-\Pi_{K}&\Sigma_{K}^{s\tilde{y}}\\\
{\Sigma_{K}^{s\tilde{y}}}^{\top}&\Sigma_{K}^{\tilde{y}}\end{bmatrix}\geq\mathbf{0},\\\\[2.84526pt]
&\Sigma_{K}^{v}>\mathbf{0},\\\\[2.84526pt]
&\Sigma_{K}^{j}>\mathbf{0},\\\\[2.84526pt]
&\Sigma+\Sigma^{v}_{k}+CL\Sigma^{v}_{k}L^{\top}C^{\top}+CB\Sigma^{j}_{k}B^{\top}C^{\top}<\beta^{*}\Sigma,\\\\[2.84526pt]
&\beta^{*}=\frac{\alpha}{2P^{-1}\left(\frac{n_{y}}{2},1-\mathcal{A}^{*}-\epsilon\right)}.\end{aligned}\right.$

Theorem $1$ provides a tool to find optimal privacy-preserving mechanisms, as
in (13), to achieve a trade-off between privacy and anomaly detector
performance when the prior distributions are jointly Gaussian. In the
following, we significantly generalize the class of stochastic dynamical
systems that Theorem $1$ can deal with.

######  Remark 1.

In the following proposition, we prove that we can employ Theorem $1$ to find
optimal privacy-preserving mechanisms for prior distributions that are not
necessarily Gaussian, as long as they are log-concave _[27]_. In particular,
we prove that if $s^{K}$ and $\tilde{y}^{K}$ follow non-Gaussian log-concave
distributions, their mutual information is upper-bounded by an affine function
of the mutual information between Gaussian random vectors with the same mean
and covariance matrices as $s^{K}$ and $\tilde{y}^{K}$, respectively. Examples
of log-concave distributions are normal, exponential, uniform, Laplace, chi,
beta, and logistic distributions. We refer the reader to _[46, 47]_ for a
comprehensive list and properties of log-concave distributions.

######  Proposition 1.

Let $s^{K}$ and $\tilde{y}^{K}$ be jointly distributed following log-concave
multivariate probability distribution. The mutual information
$I[s^{K};\tilde{y}^{K}]$ can be upper-bounded as follows:

$\displaystyle I[s^{K};\tilde{y}^{K}]\leq
I[s^{K}_{\mathcal{N}};\tilde{y}^{K}_{\mathcal{N}}]+C_{n},$ (68)

where $C_{n}=\frac{1}{2}\log(2\pi ec(n))^{n}$, $n=K(n_{s}+n_{y})$,
$c(n)=\frac{e^{2}n^{2}}{4\sqrt{2}(n+2)}$, and $s^{K}_{\mathcal{N}}$ and
$\tilde{y}^{K}_{\mathcal{N}}$ denote Gaussian random vectors with the same
covariance matrices as $s^{K}$ and $\tilde{y}^{K}$.

_Proof_ : The proof is given in Appendix C. $\blacksquare$

######  Remark 2.

Proposition 1 implies that by minimizing
$I[s^{K}_{\mathcal{N}};\tilde{y}^{K}_{\mathcal{N}}]-h[j^{K}]$ (based on the
optimization problem in Theorem 1), we can effectively decrease the
information leakage $I[s^{K};\tilde{y}^{K}]-h[j^{K}]$ for any
$(s^{K},\tilde{y}^{K})$ following log-concave distributions.

A tighter bound for the mutual information of log-concave random vectors can
be achieved based on their maximum density. The maximum density of a random
vector is calculated by its $L_{\infty}$ norm [48]. In the following
proposition, we prove that the mutual information between log-concave random
vectors $s^{K}$ and $\tilde{y}^{K}$ can be upper-bounded by an affine function
of the mutual information of Gaussian random vectors with maximum density
being the same as those of $s^{K}$ and $\tilde{y}^{K}$.

######  Proposition 2.

Let $s^{K}$ and $\tilde{y}^{K}$ be jointly distributed following a log-concave
multivariate joint probability distribution with mutual information
$I[s^{K};\tilde{y}^{K}]$, then:

$\displaystyle I[s^{K};\tilde{y}^{K}]\leq
I[s^{K}_{{\mathcal{N}}^{*}};\tilde{y}^{K}_{{\mathcal{N}}^{*}}]+n,$ (69)

where $n=K(n_{s}+n_{y})$, and $s^{K}_{{\mathcal{N}}^{*}}$ and
$\tilde{y}^{K}_{{\mathcal{N}}^{*}}$ denote multivariate normally distributed
random vectors with maximum density being the same as those of $s^{K}$ and
$\tilde{y}^{K}$, respectively.

_Proof_ : The proof is given in Appendix D. $\blacksquare$

######  Remark 3.

In Proposition 2 we prove that the difference between $I[s^{K};\tilde{y}^{K}]$
and $I[s^{K}_{{\mathcal{N}}^{*}};\tilde{y}^{K}_{{\mathcal{N}}^{*}}]$ can be
upper-bounded by $n$. This upper bound is still conservative since it needs to
cover all log-concave distributed random vectors. However, the authors in [48]
observe that every log-concave random vector has approximately the same
entropy per coordinate as a related Gaussian vector (which resulted in
approximately the same mutual information), and log-concave distributions
resemble Gaussian distributions. Therefore, by minimizing
$I[s^{K}_{{\mathcal{N}}^{*}};\tilde{y}^{K}_{{\mathcal{N}}^{*}}]$ (based on the
optimization problem in Theorem 1), we can effectively decrease the
information leakage $I[s^{K};\tilde{y}^{K}]-h[j^{K}]$ for any
$(s^{K},\tilde{y}^{K})$ following a log-concave distribution.

So far, we have addressed the problem of designing optimal probabilistic
mappings to maximize privacy for a finite time horizon. In the next section,
we prove that the proposed mechanisms can be implemented for an arbitrarily
large time horizon by solving finite horizon problems repeatedly and
sequentially implementing the obtained mechanisms.

##  VI Arbitrary Large Time Horizon

For the setting given in Theorem 1, we have addressed the problem of designing
optimal probabilistic privacy-preserving mechanisms to maximize privacy while
considering a constraint on the false alarm rate for a finite time horizon
$\mathcal{K}$. In this section, we seek to prove that the arbitrary large
horizon sub-optimal mappings can be designed by repeatedly solving finite time
problems of the form (67). To this goal, at first, we aim to evaluate the
effect of splitting time horizon $\mathcal{K}$ into subwindows of time in the
optimization problem on the amount of information leakage, which is equivalent
to the optimal cost function in optimization problem (67).

######  Lemma 5.

If we split a time horizon $\mathcal{K}:=\\{1,\ldots,K\\}$, $K\in\mathbb{N}$,
into $N$ subwindows of time as $\mathcal{K}_{1}:=\\{1,\ldots,K_{1}\\}$,
$\mathcal{K}_{2}:=\\{K_{1}+1,\ldots,K_{2}\\}$,…,$\mathcal{K}_{N}:=\\{K_{N-1}+1,\ldots,K_{N}=K\\}$,
and then, solve the optimization problem (67) in each subwindow of time, the
optimal information leakage will be more than optimal information leakage
without splitting time _:_

$\gamma_{\mathcal{K}}\leq\gamma_{\mathcal{K}_{1}}+\gamma_{\mathcal{K}_{2}}+\ldots+\gamma_{\mathcal{K}_{N}},$
(70)

where
$\gamma_{\mathcal{K}},\gamma_{\mathcal{K}_{1}},\ldots,\gamma_{\mathcal{K}_{N}}$
are the optimal information leakage in
$\mathcal{K},\mathcal{K}_{1},\ldots,\mathcal{K}_{N}$, respectively.

_Proof_ : The proof is given in Appendix E. $\blacksquare$

######  Remark 4.

From Lemma 5, we can conclude that the amount of information leakage in a time
horizon is upper bounded by the summation of the amount of information leakage
in each subwindow of time. Therefore, arbitrary large horizon sub-optimal
mappings can be designed by repeatedly solving splitted time problems of form
(67).

######  Remark 5.

In this work, we consider the problem of solving an optimization problem to
hold a trade-off between privacy and the performance of the anomaly detection
algorithm in order to design optimal privacy-distorting mechanisms. For
designing the privacy mechanisms in the general case, i.e., for utilities of
data other than anomaly detection, another optimization problem to hold a
trade-off between _privacy_ and the amount of _distortion_ is described in
[26, 49]. This generalizes the applicability of the results in the paper to a
much wider range of applications (beyond anomaly detection) in which a balance
between privacy and performance (challenged by the distortion of data) is
needed.

##  VII On Anomaly Detection Performance Metrics

In order to assess the cost of the proposed privacy-preserving mechanisms in
terms of detectability loss, three different metrics are employed.

The first metric for evaluating the performance of the anomaly detector is the
false alarm rate. In the absence of privacy mechanisms, the false alarm rate
is $\mathcal{A}^{*}$. After applying the privacy mechanisms, based on
distortion constraint
$\operatorname{Pr}[\tilde{z}_{k}>\alpha]<\mathcal{A}^{*}+\epsilon$ in
optimization problem (21), the false alarm rate
$\mathcal{A}=\operatorname{Pr}[\tilde{z}_{k}>\alpha]$ can be increased by
$\epsilon$. Therefore, the privacy distorting mechanism degrades the
performance of the anomaly detector in terms of the false alarm rate. We
calculate the false alarm rate based on the CDF of $\tilde{z}_{k}$ as follows:

$\displaystyle\mathcal{A}=\operatorname{Pr}\left[\tilde{z}_{k}>\alpha|\delta_{k}=\mathbf{0}\right]=1-F_{\tilde{z}_{k}}\left(\alpha|\delta_{k}=\mathbf{0}\right).$
(71)

Since $\tilde{z}_{k}$ approximately follows a gamma distribution (see Appendix
F), its CDF can be determined by the inverse of the lower incomplete gamma
function [44].

Another metric is the detection rate of the anomaly detector for deterministic
additive fault $\delta_{k}$. The detection rate to detect faulty systems with
distance measure ${z}_{k}$ is defined as
$\mathcal{B}=\operatorname{Pr}\left[{z}_{k}>\alpha|\delta_{k}\neq\mathbf{0}\right]$.
We start with the case when there is no privacy distortion in the system but
$\delta_{k}\neq\mathbf{0}$. In this case, the faulty residual sequence is
determined by (16). According to the superposition principle, the faulty
residual sequence can be written as:

$\displaystyle r_{k}=r_{k}^{\delta}+r_{k}^{n},$ (72)

where $r_{k}^{\delta}$ is a signal representing the change in the residual due
to the fault $\delta_{k}$, and $r_{k}^{n}$ represents the change in the
residual due to system and measurement noises. In (9), we proved that if there
is no anomaly and no privacy distortion in the system, the residual follows
Gaussian distribution ${r}_{k}\sim\mathcal{N}(\mathbf{0},\Sigma)$. Besides,
since $\delta_{k}$ is a deterministic additive fault, $r_{k}^{\delta}$ will
only appear in the mean of the residual. Therefore, the residual follows
Gaussian distribution ${r}_{k}\sim\mathcal{N}(r_{k}^{\delta},\Sigma)$ when
$\delta_{k}\neq\mathbf{0}$. It can be proved that based on this new residual,
the new distance measure is distributed as a non-central chi-squared with
$n_{y}$ degrees of freedom and non-central parameter
$\lambda_{k}=\left\|\Sigma^{(-1/2)}r_{k}^{\delta}\right\|^{2}$ [50]. Note that
for a fixed $\Sigma$, $\lambda_{k}$ is an increasing function of the fault.
Let $\mathcal{B}^{*}$ be the detection rate without privacy distortion. This
detection rate can be calculated based on the CDF of $z_{k}$ as

$\displaystyle\mathcal{B}^{*}=1-F_{{z}_{k}}\left(\alpha|\delta_{k}\neq\mathbf{0}\right).$
(73)

Since ${z}_{k}$ is a non-central chi-squared variable for
$\delta_{k}\neq\mathbf{0}$, and due to the fact that the CDF of a non-central
chi-squared variable is a decreasing function of its non-central parameter
(see [51]), the CDF of ${z}_{k}$,
$F_{{z}_{k}}\left(\alpha|\delta_{k}\neq\mathbf{0}\right)$, is a decreasing
function of $\lambda_{k}$. As a result, by increasing $\lambda_{k}$,
$\mathcal{B}^{*}=1-F_{{z}_{k}}\left(\alpha|\delta_{k}\neq\mathbf{0}\right)$
increases. Therefore, given the fact that $\mathcal{B}^{*}$ is an increasing
function of $\lambda_{k}$ and $\lambda_{k}$ is an increasing function of
$\delta_{k}$, we can conclude that for large fault $\delta_{k}$,
$\mathcal{B}^{*}$ will be increased which amounts better detectability.

To calculate the detection rate for the case when there is both privacy
distortion and fault in the system, based on the superposition principle, we
need to consider another term in the residual for the additive noises for
privacy compared with (72). Therefore, the distorted faulty residual is given
by:

$\displaystyle\tilde{r}_{k}=r_{k}^{\delta}+r_{k}^{n}+r_{k}^{p},$ (74)

where $r_{k}^{p}$ is a signal representing the change in the residual due to
privacy distortions. In Equation (17), we showed that the distorted residual
follows $\tilde{r}_{k}\sim\mathcal{N}(\mathbf{0},\tilde{\Sigma}_{k})$.
Consequently, the distorted faulty residual follows Gaussian distribution as
$\tilde{r}_{k}\sim\mathcal{N}(r_{k}^{\delta},\tilde{\Sigma}_{k})$. It can be
proved that based on the new residual, the new distance measure
$\tilde{z}_{k}={\tilde{r}}^{\top}_{k}\Sigma^{-1}\tilde{r}_{k}$ is distributed
as a generalized chi-squared [52]. Then, the detection rate $\mathcal{B}$ in
the presence of privacy mechanisms can be calculated as follows:

$\displaystyle\mathcal{B}=1-F_{\tilde{z}_{k}}\left(\alpha|\delta_{k}\neq\mathbf{0}\right).$
(75)

Since $\tilde{z}_{k}$ follows generalized chi-squared, its PDF and CDF,
unfortunately, cannot be expressed in closed form, although approximations are
available [52]. Therefore, in this case, the detection rate needs to be
calculated numerically. That is why this metric cannot effectively be employed
to define the detection performance constraint in Problem $1$.

To determine the effect of privacy distortion on the detection rate, the
function of the detection rate in terms of $\delta_{k}$ in with and without
privacy distortion cases can be compared. At the start point, when
$\delta_{k}=\mathbf{0}$, the detection rate equals the false alarm rate, and
therefore, the detection rate is higher when privacy distortions are added.
Besides, as we mentioned before, when there is no privacy distortion, the
detection rate is monotonically increasing for an increasing amount of
$\delta_{k}$. In addition, when $\delta_{k}$ goes to infinity, the detection
rate goes to $1$ for both with and without privacy distortions cases. Hence,
by adding privacy distortions, the function of detection rate in terms of
fault starts from a higher point (compared to without privacy distortions
case) and reaches the same point for infinitely large fault where
$\mathcal{B}=1$. Hence, the privacy distortion flattens the detection rate
function in terms of fault magnitude (see Fig. 6 for our implemented numerical
example of this occurrence).

To investigate the overall effect of privacy mechanisms on detection
performances, the Receiver Operating Characteristic (ROC) is considered. The
ROC curves show how the detection rate changes as the false-alarm rate
increases. In (71) and (75), the false alarm rate and the detection rate are
written in terms of $\alpha$, where $\alpha$ is selected to satisfy the
desired false alarm rate $\mathcal{A}^{*}$. Therefore, calculating the false
alarm rate and detection rate for different $\mathcal{A}^{*}$ will allow
characterizing the ROC of the proposed detection scheme. The ROC plot can be
employed to explore the net effect of privacy mechanisms on detection
performances in terms of false alarm rate and detection rate.

##  VIII Illustrative case study

The authors of [53, 37] studied the anomaly detection problem for a well-
stirred chemical reactor with a heat exchanger. We use this case study to
demonstrate our results. The state, input, and output vectors of the
considered reactor are:

$\displaystyle\left\\{\begin{array}[]{ll}x_{k}=\begin{pmatrix}C_{0}\\\
T_{0}\\\ T_{w}\\\ T_{m}\end{pmatrix},u_{k}=C_{u},y_{k}=\begin{pmatrix}C_{0}\\\
T_{0}\end{pmatrix},\end{array}\right.$

where

$\left\\{\begin{array}[]{ll}\begin{aligned} &C_{0}:\text{Concentration of the
chemical product},\\\ &T_{0}:\text{Temperature of the product},\\\
&T_{w}:\text{Temperature of the jacket water of the heat exchanger},\\\
&T_{m}:\text{Coolant temperature},\\\ &C_{u}:\text{Inlet concentration of the
reactant}.\end{aligned}\end{array}\right.$

Figure 2: Evolution of the optimal cost function (information leakage) based
on increasing $\epsilon$ for horizon $K=60$.

Figure 3: Comparison between the first element of the measurement vector
$y_{k}$ and the distorted measurements $\tilde{y}_{k}$ for different
distortion levels $\epsilon\in\\{0.1,0.75,0.9\\}$, and $\mathcal{A}^{*}=0.1$.

$\left\\{\begin{array}[]{l}A=\left(\begin{array}[]{cccc}0.8353&0&0&0\\\
0&0.8324&0&0.0031\\\ 0&0.0001&0.1633&0\\\
0&0.0280&0.0172&0.9320\end{array}\right),\quad
B=\left(\begin{array}[]{ccc}0.0458\\\ 0\\\ 0\\\ 0\end{array}\right),\quad
C=\left(\begin{array}[]{cccc}1&0&0&0\\\
0&1&0&0\end{array}\right),\quad\Sigma^{w}=0.01I_{n_{y}},\\\
\Sigma^{t}=\left(\begin{array}[]{cccc}0.1&0&0&0\\\ 0&0.2&0&0\\\ 0&0&0.3&0\\\
0&0&0&0.4\end{array}\right),\quad L=\left(\begin{array}[]{ccc}0.8271&0\\\
0&0.8243\\\ 0&0.0002\\\
0&0.0481\end{array}\right),\quad\Sigma=\left(\begin{array}[]{ccc}1.0169&0\\\
0&1.0169\end{array}\right).\end{array}\right.$ (76)

We use the discrete-time dynamics of the reactor introduced in [25] for this
illustrative simulation study with matrices as given in (76), normally
distributed initial condition
$x_{1}\sim\mathcal{N}[(6.94;13.76;1;1)^{\top},I_{4}]$, and reference signal
$u_{k}=50\cos[0.5k]^{2}$. As private output, we use the concentration of the
chemical product; then, matrix $D$ in (4) is given by the full row rank matrix
$D=(1,0,0,0)$. Therefore, the aim of the privacy scheme is to hide the private
state $C_{0}$, which is the concentration of the reactant, as much as possible
without distorting anomaly detection performance excessively.

First, in Figure 2, we show the effect of false alarm distortion level
$\epsilon$ on the amount of optimal information leakage, which is shown by the
optimal cost $I[s^{K};\tilde{y}^{K}]-h[j^{K}]$. This figure depicts the
evolution of the optimal cost for increasing $\epsilon$ with time horizon
$K=60$ and desired false alarm rate $\mathcal{A}^{*}=0.1$ in the absence of
fault. As expected, by increasing the distortion level $\epsilon$, which
increases the amount of distortion on the false alarm rate as we discussed in
Section 7 (i.e., leading to poorer anomaly detection performance), the
information leakage decreases monotonically (i.e., leading to improved
privacy). This figure illustrates that while designing the privacy distorting
mechanisms, we need to balance the trade-off between privacy and anomaly
detector performance.

The effect of optimal distorting mechanisms is shown in Figure 3, where we
contrast actual and distorted measurement data for different distortion levels
$\epsilon\in\\{0.1,0.75,0.9\\}$. We consider $\mathcal{A}^{*}=0.1$; therefore,
$\epsilon=0.9$ means that the optimization problem in (67) is solved without
considering the distortion constraint (because in this case, the upper bound
of the false alarm rate is $\mathcal{A}^{*}+\epsilon=1$ which means that there
is no limitation on the false alarm rate). As can be seen in these figures,
the difference between distorted and actual measurements will increase by
increasing the distortion level. The same result is achieved for the
comparison of actual and distorted inputs.

Next, in Figure 4, we depict the realization of $s_{k}$, the estimated private
output without distorting $y_{k}$ and $u_{k}$ using the MMSE estimator,
${{\hat{s}}_{k}^{yu}}$, and the estimated private output given the distorted
vectors $\tilde{y}_{k}$ and ${\tilde{u}}_{k}$,
${{\hat{s}}_{k}^{\tilde{y}\tilde{u}}}$. We consider time horizon $K=60$,
distortion levels $\epsilon\in\\{0.1,0.5,0.9\\}$, and desired false alarm rate
$\mathcal{A}^{*}=0.1$ in the absence of fault. As can be seen in these
figures, by increasing the distortion level $\epsilon$, the accuracy of
estimation decreases. This figure demonstrates that by randomizing measurement
and input data, we can prevent worst-case adversaries from accurately
estimating $s_{k}$.

Figure 4: Comparison between the private output $s_{k}$, the estimated private
output without distorting $y_{k}$ and $u_{k}$, ${{\hat{s}}_{k}^{yu}}$, and the
estimated private output given distorted data $\tilde{y}_{k}$ and
$\tilde{u}_{k}$, ${{\hat{s}}_{k}^{\tilde{y}\tilde{u}}}$, for horizon $K=60$,
desired false alarm rate $\mathcal{A}^{*}=0.1$, and different distortion
levels, $\epsilon\in\\{0.1,0.5,0.9\\}$.

To evaluate the effect of privacy mechanisms on anomaly detection performance,
we apply three different metrics, as discussed in Section 7. First, in Figure
5, the effect of distortion level $\epsilon$ on the false alarm rate is shown.
As discussed in Section 7, we employ the CDF of the approximation of
$\tilde{z}_{k}$, $F_{\tilde{z}_{k}}(\alpha)$, using the lower incomplete gamma
function to determine the false alarm rate. As can be seen, by increasing
$\epsilon$, the false alarm rate increases. Therefore as we expected,
increasing false alarm distortion level $\epsilon$, which resulted in
increasing privacy distortion and decreasing the information leakage, indeed
reduces the performance of the anomaly detector in terms of false alarm rate.
Besides, this figure shows that the upper bound we designed for the false
alarm rate in distortion constraint (20), $\mathcal{A}^{*}+\epsilon$, is
respected and not so conservative.

Next, as the second detection metric, the effect of privacy distortion on the
detection rate of the chi-squared tuning procedure is depicted in Figure 6. We
consider deterministic additive constant fault $\delta_{k}=\delta$ with fault
matrices $H=(0,1)^{\top}$ and $G=\mathbf{0}$. We apply privacy distortions
with different distortion levels $\epsilon\in\\{0,0.1,0.3,0.5\\}$ and desired
false alarm rate $\mathcal{A}^{*}=0.1$, where $\epsilon=0$ indicates without
privacy distortion. As can be seen in this figure, the detection rate of the
anomaly detection algorithm $\mathcal{B}$, as defined in Section 7, is
changing shape and flattening by increasing the distortion level $\epsilon$.
For small $\delta$, the detection rate increases by increasing $\epsilon$, and
for large $\delta$, the detection rate slightly decreases by increasing
$\epsilon$. It might be mistakenly concluded from this figure that increasing
the distortion level can improve the detection performance of the anomaly
detector. On the other hand, Figure 5 illustrates that by increasing the
distortion level, the false alarm rate will increase. Therefore, to
investigate the overall effect of privacy mechanisms on detection
performances, the Receiver Operating Characteristic (ROC) is considered for
analyzing the anomaly detector performance.

The ROC shows how the detection rate changes as the false alarm rate
increases. To characterize the ROC of the anomaly detection scheme, we find
the detection rate (true positive) for a given fault $\delta$ as a function of
the false alarm rate (false positive). Generally, in ROC analysis, tests are
categorized based on the area under the ROC curve. The closer a ROC curve is
to the upper left corner, the more efficient the test is. At first, in Figure
7, we plot the ROC when there is no privacy distortion in the system for
different amounts of fault $\delta\in\\{0.1,1,2,3,4\\}$. As we expected based
on (73), by increasing the fault, the detection rate will be increased, and
the ROC curve will be closer to the upper left corner, which illustrates the
better performance of the anomaly detector. Then in Figure 8, the ROC is
depicted for the distorted faulty system with a fixed fault $\delta=2$ and
different distortion levels $\epsilon\in\\{0.01,0.3,0.5\\}$ by increasing the
desired false alarm rate. As can be seen, for the fixed false alarm rate, the
detection rate is decreased by increasing $\epsilon$. Also, for the fixed
detection rate, the false alarm rate is increased by increasing $\epsilon$. In
other words, by increasing $\epsilon$, the ROC curve will get farther from the
upper left corner, which illustrates the worse performance of the anomaly
detector. To further explore this figure, we marked with circle markers on all
plots corresponding to a specific value of the desired false alarm rate
$\mathcal{A}^{*}=0.3$. It can be seen that by increasing $\epsilon$, the
$\mathcal{A}^{*}=0.3$ point moves toward the top right corner of the ROC,
close to point $(1,1)$ which shows a severely poor detection performance,
equivalent to that of a detector that always identifies a residual as faulty,
regardless of its real value. Therefore, this figure shows that increasing the
privacy distortion will reduce the anomaly detector’s performance.

In Figure 9, we show that, as discussed in Remark 2 and 3, for non-Gaussian
random vectors $s^{K}$ and $\tilde{y}^{K}$ following joint log-concave
distributions, their mutual information can be upper-bounded by affine
functions of the mutual information of some Gaussian random vectors with the
same mean and covariance matrices. We consider $s^{K}$ and $\tilde{y}^{K}$
with joint Laplace distribution and solve the optimization problem in (67) for
$I[s^{K}_{\mathcal{N}};\tilde{y}^{K}_{\mathcal{N}}]-h[j^{K}]$ to find the
optimal distorting parameters for an increasing amount of distortion level
$\epsilon$. Since the maximum density [48] of Gaussian random vectors and
Laplace random vectors are the same for the same means and covariances, we can
conclude that for Laplace variables
$s^{K}_{{\mathcal{N}}^{*}}=s^{K}_{\mathcal{N}}$,
$\tilde{y}^{K}_{{\mathcal{N}}^{*}}=\tilde{y}^{K}_{\mathcal{N}}$, and therefore
$I[s^{K}_{{\mathcal{N}}^{*}};\tilde{y}^{K}_{{\mathcal{N}}^{*}}]-h[j^{K}]=I[s^{K}_{\mathcal{N}};\tilde{y}^{K}_{\mathcal{N}}]-h[j^{K}]$.
We calculate the optimal
$I[s^{K}_{\mathcal{N}};\tilde{y}^{K}_{\mathcal{N}}]=I[s^{K}_{{\mathcal{N}}^{*}};\tilde{y}^{K}_{{\mathcal{N}}^{*}}]$
by solving the optimization problem (67). Then, we calculate
$I[s^{K};\tilde{y}^{K}]$ (see [54] for the mutual information calculation),
where $s^{K}$ and $y^{K}$ follow Laplace distributions with the same mean and
covariance as $s^{K}_{\mathcal{N}}$ and $y^{K}_{\mathcal{N}}$, and
$\tilde{y}^{K}=y^{K}+v^{K}$ (with the same distortion variables as the
Gaussian case). As can be seen in Figure 9, $I[s^{K};\tilde{y}^{K}]$ is even
upper-bounded by $I[s^{K}_{\mathcal{N}};\tilde{y}^{K}_{\mathcal{N}}]$ without
the constant terms $n$ and $C_{n}$. It implies that by solving the
optimization problem for Gaussian variables, the information leakage for
Laplace variables decreases for increasing $\epsilon$ and is upper-bounded by
the information leakage for Gaussian variables. Furthermore, it is apparent
that even if the Gaussian upper bounds are loose for these distributions of
$s^{K}$ and $\tilde{y}^{K}$, the information leakage is still decreasing for
increasing $\epsilon$. It shows that minimizing
$I[s^{K}_{\mathcal{N}};\tilde{y}^{K}_{\mathcal{N}}]-h[j^{K}]$ or
$I[s^{K}_{{\mathcal{N}}^{*}};\tilde{y}^{K}_{{\mathcal{N}}^{*}}]-h[j^{K}]$
effectively decreases the information leakage for any log-concave $s^{K}$ and
$\tilde{y}^{K}$.

Finally, in Figure 10, we illustrate the effect of splitting time horizon into
subwindows of time in the optimization problem (67) on the amount of
information leakage. At first, we solve the optimization problem in time
horizon $\mathcal{K}=\\{1,\ldots,60\\}$ for an increasing amount of
$\epsilon$. Then, we divide $\mathcal{K}$ into three and then six time
subwindows and solve the optimization problem in them, separately, i.e., for
three divisions, time horizons are $\\{1,2,\ldots,20\\}$,
$\\{21,22,\ldots,40\\}$, and $\\{41,42,\ldots,60\\}$. As seen in this figure,
the optimal cost function $I[s^{K};\tilde{y}^{K}]-h[j^{K}]$ is increased by
increasing the number of divisions. It shows that the optimal information
leakage in $\mathcal{K}$ without splitting the time horizon is upper bounded
by suboptimal information leakage in $\mathcal{K}$ by solving the optimization
problem in subwindows of time. It deduced that the amount of information
leakage in an arbitrarily large time horizon would be upper bounded by the
summation of the amount of information leakage in each subwindow of time, and
arbitrary large horizon sub-optimal mappings can be designed by repeatedly
solving finite time problems of the form (67).

Figure 5: Comparison between expected false alarm rate
$\mathcal{A}^{*}+\epsilon$ and approximation of false alarm rate after privacy
distortions $\mathcal{A}$ for an increasing amount of $\epsilon$, horizon
$K=60$, and desired false alarm rate $\mathcal{A}^{*}=0.1$. Figure 6:
Comparison between the detection rate $\mathcal{B}$ of chi-squared tuning
algorithm for an increasing amount of fault $\delta$ for the system without
privacy-induced distortion ($\epsilon=0$) and system with privacy-induced
distortions (with $\epsilon\in\\{0.1,0.3,0.5\\}$) and $\mathcal{A}^{*}=0.1$.
Figure 7: ROC curves of the detection rate in terms of desired false alarm
rate without privacy distortion for different faults
$\delta\in\\{0.1,1,2,3,4\\}$. Figure 8: ROC curves of the detection rate in
terms of the false alarm rate for different distortion levels
$\epsilon\in\\{0.01,0.3,0.5\\}$, fixed fault $\delta=2$, and increasing
$\mathcal{A}^{*}$. The circle markers on each line correspond to a certain
value of the desired false alarm rate $\mathcal{A}^{*}=0.3$. Figure 9:
Comparison between the mutual information of Laplacian $s^{K}$ and
$\tilde{y}^{K}$, $I[s^{K};\tilde{y}^{K}]$, and the mutual information of
Gaussian vectors $s^{K}_{\mathcal{N}}$ and $\tilde{y}^{K}_{\mathcal{N}}$,
$I[s^{K}_{\mathcal{N}};\tilde{y}^{K}_{\mathcal{N}}]$, for increasing
$\epsilon$ and desired false alarm rate $\mathcal{A}^{*}=0.1$. Figure 10:
Comparison between optimal information leakage
$I[s^{K};\tilde{y}^{K}]-h[j^{K}]$ in time horizon
$\mathcal{K}=\\{1,\ldots,60\\}$, and the summation of optimal information
leakage in divided time horizons with three divisions and six divisions.

##  IX Conclusions

In this paper, we presented an anomaly detection framework to maximize privacy
for a class of stochastic dynamical systems. A detailed mathematical framework
for synthesizing distorting mechanisms is given to minimize the information
leakage induced by the use of public/unsecured communication networks. We
proposed a class of linear Gaussian distorting mechanisms to randomize sensor
measurements and input signals before transmission to prevent adversaries from
accurately estimating the private part of the system state (a private
performance output). Furthermore, for the class of systems under study, we
thoroughly characterized information-theoretic metrics (mutual information and
differential entropy) to quantify the information between private outputs and
disclosed data for a class of worst-case eavesdropping adversaries.

Anomaly detection is an important problem that has been studied in a variety
of applications. We proved that the privacy distorting mechanisms can increase
the false alarm rate in the anomaly detector, resulting in reducing the
performance of the anomaly detection scheme. Therefore, when designing the
privacy distorting mechanism, we considered the trade-off between _privacy_
and _anomaly detection performance_. We formulated the synthesis of the
distorting mechanisms as a convex program, where we minimize the information
leakage between an arbitrarily large sequence of private outputs and the
disclosed distorted data while considering a constraint on the amount of false
alarm rate in the anomaly detector to limit the effect of privacy mechanisms
on the anomaly detection algorithm. We presented simulation results to
illustrate the performance of the proposed tools.

## References

  * [1] S. R. Rajagopalan, L. Sankar, S. Mohajer, and H. V. Poor, “Smart meter privacy: A utility-privacy framework,” in _Proceedings of the IEEE International Conference on Smart Grid Communications (SmartGridComm)_ , 2011, pp. 190–195.
  * [2] O. Tan, D. Gündüz, and H. V. Poor, “Increasing smart meter privacy through energy harvesting and storage devices,” _IEEE Journal on Selected Areas in Communications_ , vol. 31, pp. 1331–1341, 2013.
  * [3] Z. Huang, Y. Wang, S. Mitra, and G. E. Dullerud, “On the cost of differential privacy in distributed control systems,” in _Proceedings of the 3rd International Conference on High Confidence Networked Systems_ , 2014, pp. 105–114.
  * [4] F. Farokhi and H. Sandberg, “Optimal privacy-preserving policy using constrained additive noise to minimize the fisher information,” in _Proceedings of the IEEE 56th Annual Conference on Decision and Control (CDC)_ , 2017.
  * [5] F. Farokhi, H. Sandberg, I. Shames, and M. Cantoni, “Quadratic Gaussian privacy games,” in _Proceedings of the 54th IEEE Conference on Decision and Control (CDC)_ , 2015, pp. 4505–4510.
  * [6] F. Miao, Q. Zhu, M. Pajic, and G. J. Pappas, “Coding sensor outputs for injection attacks detection,” in _Proceedings of the 53rd Annual Conference on Decision and Control (CDC)_ , 2014, pp. 5776–5781.
  * [7] J. L. Ny and G. J. Pappas, “Differentially private filtering,” _IEEE Transactions on Automatic Control_ , vol. 59, pp. 341–354, 2014.
  * [8] R. M. Ferrari and A. M. Teixeira, _Safety, Security and Privacy for Cyber-Physical Systems_. Springer, 2021\.
  * [9] F. du Pin Calmon and N. Fawaz, “Privacy against statistical inference,” in _Proceedings of the 50th Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ , 2012, pp. 1401–1408.
  * [10] Q. Geng and P. Viswanath, “The optimal mechanism in differential privacy,” in _Proceedings of the IEEE International Symposium on Information Theory_ , 2014, pp. 2371–2375.
  * [11] T. Tanaka, M. Skoglund, H. Sandberg, and K. H. Johansson, “Directed information as privacy measure in cloud-based control,” _arXiv preprint arXiv:1705.02802_ , 2017.
  * [12] F. F. N. D. Murguia C., Shames I., _Information-Theoretic Privacy Through Chaos Synchronization and Optimal Additive Noise. In: Farokhi F. (eds) Privacy in Dynamical Systems_. Singapore: Springer, 2020.
  * [13] C. Murguia, I. Shames, F. Farokhi, and D. Nešić, “On privacy of quantized sensor measurements through additive noise,” in _Proceedings of the 57th IEEE Conference on Decision and Control (CDC)_ , 2018.
  * [14] C. Murguia and J. Ruths, “On model-based detectors for linear time-invariant stochastic systems under sensor attacks,” _IET Control Theory & Applications_, vol. 13, no. 8, pp. 1051–1061, 2019.
  * [15] J. J. Gertler, “Survey of model-based failure detection and isolation in complex plants,” _IEEE Control systems magazine_ , vol. 8, no. 6, pp. 3–11, 1988.
  * [16] T. M. Cover and J. A. Thomas, _Elements of Information Theory_. New York, NY, USA: Wiley-Interscience, 1991.
  * [17] C. Dwork, “Differential privacy: A survey of results,” in _Theory and Applications of Models of Computation_. Berlin, Heidelberg: Springer Berlin Heidelberg, 2008, pp. 1–19.
  * [18] C. Dwork and A. Roth, “The algorithmic foundations of differential privacy,” _Foundations and Trends in Theoretical Computer Science_ , vol. 9, pp. 211–407, 2014.
  * [19] S. Han, U. Topcu, and G. J. Pappas, “Differentially private convex optimization with piecewise affine objectives,” in _Proceedings of the 53rd IEEE Conference on Decision and Control_ , 2014.
  * [20] J. Soria-Comas and J. Domingo-Ferrer, “Optimal data-independent noise for differential privacy,” _Information Sciences_ , vol. 250, pp. 200 – 214, 2013.
  * [21] Y. Wang, Z. Huang, S. Mitra, and G. E. Dullerud, “Entropy-minimizing mechanism for differential privacy of discrete-time linear feedback systems,” in _Proceedings of the 53rd IEEE Conference on Decision and Control_ , 2014, pp. 2130–2135.
  * [22] F. Farokhi and G. Nair, “Privacy-constrained communication,” _IFAC-PapersOnLine_ , vol. 49, pp. 43 – 48, 2016.
  * [23] E. V. Belmega, L. Sankar, and H. V. Poor, “Enabling data exchange in two-agent interactive systems under privacy constraints,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 9, pp. 1285–1297, 2015.
  * [24] E. Akyol, C. Langbort, and T. Basar, “Privacy constrained information processing,” in _Proceedings of the 54th IEEE Conference on Decision and Control (CDC)_ , 2015, pp. 4511–4516.
  * [25] C. Murguia, I. Shames, F. Farokhi, D. Nešić, and H. V. Poor, “On privacy of dynamical systems: An optimal probabilistic mapping approach,” _IEEE Transactions on Information Forensics and Security_ , 2021.
  * [26] H. Hayati, C. Murguia, and N. Van De Wouw, “Finite horizon privacy of stochastic dynamical systems: A synthesis framework for gaussian mechanisms,” in _2021 60th IEEE Conference on Decision and Control (CDC)_. IEEE, 2021, pp. 5607–5613.
  * [27] A. Prékopa, “Logarithmic concave measures and related topics,” in _Stochastic programming_. Citeseer, 1980.
  * [28] R. M. Ferrari, K. H. Degue, and J. Le Ny, “Differentially private anomaly detection for interconnected systems,” in _Safety, Security and Privacy for Cyber-Physical Systems_. Springer, 2021, pp. 203–230.
  * [29] M. Huemer, O. Lang, and C. Hofbauer, “Component-wise conditionally unbiased widely linear mmse estimation,” _Signal Processing_ , vol. 133, pp. 227–239, 2017.
  * [30] K. J. Aström and B. Wittenmark, _Computer-controlled Systems (3rd Ed.)_. Upper Saddle River, NJ, USA: Prentice-Hall, Inc., 1997.
  * [31] A. Wald, “Sequential tests of statistical hypotheses,” _The annals of mathematical statistics_ , vol. 16, no. 2, pp. 117–186, 1945.
  * [32] A. S. Willsky, “A survey of design methods for failure detection in dynamic systems,” _Automatica_ , vol. 12, no. 6, pp. 601–611, 1976.
  * [33] E. S. Page, “Continuous inspection schemes,” _Biometrika_ , vol. 41, no. 1/2, pp. 100–115, 1954.
  * [34] C. Murguia and J. Ruths, “Characterization of a cusum model-based sensor attack detector,” in _2016 IEEE 55th Conference on Decision and Control (CDC)_. IEEE, 2016, pp. 1303–1309.
  * [35] M. Basseville, “Detecting changes in signals and systems—a survey,” _Automatica_ , vol. 24, no. 3, pp. 309–326, 1988.
  * [36] C. Murguia and J. Ruths, “Cusum and chi-squared attack detection of compromised sensors,” in _2016 IEEE Conference on Control Applications (CCA)_. IEEE, 2016, pp. 474–480.
  * [37] C. Murguia and J. Ruths, “On model-based detectors for linear time-invariant stochastic systems under sensor attacks,” _IET Control Theory Applications_ , vol. 13, pp. 1051–1061, 2019.
  * [38] S. M. Ross, _Introduction to probability models_. Academic press, 2014.
  * [39] M. Ross, _Introduction to Probability Models, Ninth Edition_. Orlando, FL, USA: Academic Press, Inc., 2006.
  * [40] F. Zhang, _The Schur complement and its applications_. Springer Science & Business Media, 2006, vol. 4.
  * [41] R. A. Horn, _Topics in Matrix Analysis_. New York, NY, USA: Cambridge University Press, 1986.
  * [42] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan, _Linear matrix inequalities in system and control theory_. SIAM, 1994.
  * [43] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [44] E. W. Stacy, “A generalization of the gamma distribution,” _The Annals of mathematical statistics_ , pp. 1187–1192, 1962.
  * [45] E. Furman and R. Zitikis, “A monotonicity property of the composition of regularized and inverted-regularized gamma functions with applications,” _Journal of Mathematical Analysis and Applications_ , vol. 348, no. 2, pp. 971–976, 2008.
  * [46] M. Bagnoli and T. Bergstrom, “Log-concave probability and its applications,” _Economic Theory_ , vol. 26, pp. 445 – 469, 2005.
  * [47] J. Wellner, “Log-concave distributions: definitions, properties, and consequences,” _Presentation, University of Paris-Diderot_ , 2012.
  * [48] S. Bobkov and M. Madiman, “The entropy per coordinate of a random vector is highly constrained under convexity conditions,” _IEEE Transactions on Information Theory_ , vol. 57, no. 8, pp. 4940–4954, 2011.
  * [49] H. Hayati, C. Murguia, and N. Van de Wouw, “Gaussian mechanisms against statistical inference: Synthesis tools,” in _2022 European Control Conference (ECC)_. IEEE, 2022, pp. 1294–1300.
  * [50] A. Feiveson and F. Delaney, _The distribution and properties of a weighted sum of chi squares_. National Aeronautics and Space Administration, 1968, vol. 4575.
  * [51] A. Samir Kamel and A. I. Abdel-Samad, “On the computation of non-central chi-square distribution function,” _Communications in Statistics-Simulation and Computation_ , vol. 19, no. 4, pp. 1279–1291, 1990.
  * [52] J.-P. Imhof, “Computing the distribution of quadratic forms in normal variables,” _Biometrika_ , vol. 48, no. 3/4, pp. 419–426, 1961.
  * [53] J. Chen and R. J. Patton, _Robust model-based fault diagnosis for dynamic systems_. Springer Science & Business Media, 2012, vol. 3.
  * [54] J. Delpiano, “Fast mutual information of two images or signals (https://www.mathworks.com/matlabcentral/fileexchange/13289-fast-mutual-information-of-two-images-or-signals),” 2021\.
  * [55] S. Boyd and L. Vandenberghe, _Convex Optimization_. New York, NY, USA: Cambridge University Press, 2004.
  * [56] S. Dharmadhikari and K. Joag-Dev, _Unimodality, convexity, and applications_. Elsevier, 1988.
  * [57] S. Hoggar, “Chromatic polynomials and logarithmic concavity,” _Journal of Combinatorial Theory, Series B_ , vol. 16, no. 3, pp. 248–254, 1974.
  * [58] O. Johnson and C. Goldschmidt, “Preservation of log-concavity on summation,” _ESAIM: Probability and Statistics_ , vol. 10, pp. 206–215, 2006.
  * [59] A. Marsiglietti and V. Kostina, “A lower bound on the differential entropy of log-concave random vectors with applications,” _Entropy_ , vol. 20, no. 3, p. 185, 2018.
  * [60] F. E. Satterthwaite, “An approximate distribution of estimates of variance components,” _Biometrics bulletin_ , vol. 2, no. 6, pp. 110–114, 1946.
  * [61] B. L. Welch, “The generalization of ‘student’s’problem when several different population varlances are involved,” _Biometrika_ , vol. 34, no. 1-2, pp. 28–35, 1947.

##  X Appendix

###  X-A Proof of Lemma 3

$\tilde{z}_{k}$ and $q_{k}$ are defined as functions of variable $m_{k}$. For
the CDF of a function of a random variable, if $Y$ and $X$ are continuous
random variables with $Y=g(X)$, the CDF of $Y$, $F_{Y}(y)$ can be
characterized based on Probability Density Function (PDF) of $X$ as follows

$\displaystyle F_{Y}(y)$ $\displaystyle=P_{Y}(Y\leq y)=P_{Y}(g(X)\leq y)$ (77)
$\displaystyle=P_{X}(x\in\mathcal{X}:g(X)\leq
y)=\int_{\\{x\in\mathcal{X}:g(X)\leq y\\}}f_{X}(x)dx,$

where $f_{X}(x)$ and $\mathcal{X}$ are PDF and the sample space of $X$.
Therefore, CDFs of $\tilde{z}_{k}$ and $q_{k}$ can be written as functions of
the PDF of $m_{k}$ as below:

$F_{\tilde{z}_{k}}(\alpha)=\int_{\\{m_{k}:m_{k}^{\top}\Sigma^{\prime}_{k}m_{k}\leq\alpha\\}}f_{m_{k}}(\eta)d\eta,$
(78)
$F_{q_{k}}(\alpha)=\int_{\\{m_{k}:\beta_{k}m_{k}^{\top}{m_{k}}\leq\alpha\\}}f_{m_{k}}(\eta)d\eta.$
(79)

From $m_{k}^{\top}\Sigma^{\prime}_{k}m_{k}\leq\beta m_{k}^{\top}m_{k}$, we can
conclude that

${\\{m_{k}:\beta_{k}m_{k}^{\top}m_{k}\leq\alpha\\}}\subseteq{\\{m_{k}:m_{k}^{\top}\Sigma^{\prime}_{k}m_{k}\leq\alpha\\}},$
(80)

which means that the ellipsoid $\\{\beta_{k}m_{k}^{\top}m_{k}=\alpha\\}$ is a
subset of ellipsoid $\\{m_{k}^{\top}\Sigma^{\prime}_{k}m_{k}=\alpha\\}$ (see
[55] for proof and details).
From (80) and the fact that $f_{m_{k}}(\eta)$ is always non-negative, we
deduce that

$\int_{\\{m_{k}:m_{k}^{\top}\Sigma^{\prime}_{k}m_{k}\leq\alpha\\}}f_{m_{k}}(\eta)d\eta\geq\int_{\\{m_{k}:\beta_{k}m_{k}^{\top}m_{k}\leq\alpha\\}}f_{m_{k}}(\eta)d\eta.$
(81)

According to (78) and (79), (81) is equivalent to
$F_{\tilde{z}_{k}}(\alpha)\geq F_{q_{k}}(\alpha)$. $\blacksquare$

###  X-B Proof of Lemma 4

By multiplying $\tilde{\Sigma}^{-1/2}_{k}$ to the right and
$(\tilde{\Sigma}^{-1/2}_{k})^{\top}$ to the left side of
$(\tilde{\Sigma}^{1/2}_{k})^{\top}\Sigma^{-1}\tilde{\Sigma}^{1/2}_{k}\leq\beta^{*}{I}_{n_{y}}$,
it will be converted to

$\Sigma^{-1}\leq\beta^{*}(\tilde{\Sigma}_{k}^{\top}\tilde{\Sigma}_{k})^{-1/2}=\beta^{*}\tilde{\Sigma}_{k}^{-1}.$
(82)

Since $\Sigma$ and $\tilde{\Sigma}_{k}$ are positive definite matrices, (82)
can be written as (66) which is a convex function of design variables
$\Sigma^{v}_{k}$ and $\Sigma^{j}_{k}$. $\blacksquare$

###  X-C Proof of Proposition 1

First, we need to prove that for initial state $x_{1}$ and disturbances
$t_{k}$ and $w_{k}$ following non-Gaussian log-concave distributions, and
$j^{K}$ and $v^{K}$ following Gaussian distributions, $s^{K}$ and
$\tilde{y}^{K}$ follow log-concave distributions. The log-concavity of
distribution is preserved under affine transformations (see [56], Lemma 2.1).
Besides, based on Hoggar’s theorem, the sum of two independent log-concave
random variables is itself log-concave [57, 58]. Therefore, based on the
system dynamics (4), distorting mechanism (13), the independency between
initial state and disturbances, and also independency between $v^{K}$,
$j^{K}$, and $y^{K}$, we can conclude that $y^{K}$, $s^{K}$, and
$\tilde{y}^{K}$, follow log-concave distributions.

Mutual information $I\left[s^{K};\tilde{y}^{K}\right]$ can be written in terms
of differential entropies as
$I[s^{K};\tilde{y}^{K}]=h[s^{K}]+h[\tilde{y}^{K}]-h[s^{K},\tilde{y}^{K}]$.
Moreover, because Gaussian distributions maximize entropy for a given
covariance matrix [16], we can write:

$\left\\{\begin{aligned} h[s^{K}]&\leq h[s^{K}_{\mathcal{N}}],\\\
h[\tilde{y}^{K}]&\leq h[\tilde{y}^{K}_{\mathcal{N}}].\end{aligned}\right.$
(83)

The authors in [59] prove that the entropy of random vectors following log-
concave distributions is lower bounded as:

$h[s^{K},\tilde{y}^{K}]\geq\frac{n}{2}\log\frac{\left|\Sigma_{K}^{\tilde{y},s}\right|^{1/n}}{c(n)},$
(84)

where $n=K(n_{s}+n_{y})$ and $c(n)=\frac{e^{2}n^{2}}{4\sqrt{2}(n+2)}$.
Expanding the lower bound in (84), we obtain

$\displaystyle\frac{n}{2}\log\frac{\left|\Sigma_{K}^{\tilde{y},s}\right|^{1/n}}{c(n)}$
$\displaystyle=\frac{n}{2}\left(\frac{1}{n}\log\operatorname{det}\Sigma_{K}^{\tilde{y},s}-\log
c(n)\right)$
$\displaystyle=\frac{1}{2}\log\operatorname{det}\Sigma_{K}^{\tilde{y},s}-\frac{n}{2}\log
c(n).$ (85)

Let $(s^{K}_{\mathcal{N}},\tilde{y}^{K}_{\mathcal{N}})$ be jointly Gaussian
with covariance $\Sigma_{K}^{\tilde{y},s}$. Then, by Definition $2$, we can
write its differential entropy as
$h[s^{K}_{\mathcal{N}},\tilde{y}^{K}_{\mathcal{N}}]=\frac{1}{2}\log\det\left(\Sigma_{K}^{\tilde{y},s}\right)+\frac{n}{2}+\frac{n}{2}\log(2\pi)$.
Therefore, using the latter formula for
$h[s^{K}_{\mathcal{N}},\tilde{y}^{K}_{\mathcal{N}}]$ and (X-C), we can write
the following

$\frac{n}{2}\log\frac{\left|\Sigma_{K}^{s,\tilde{y}}\right|^{1/n}}{c(n)}=h[s^{K}_{\mathcal{N}},\tilde{y}^{K}_{\mathcal{N}}]-\frac{n}{2}-\frac{n}{2}\log(2\pi)-\frac{n}{2}\log
c(n).$ (86)

Combining (83), (84), and (86), we can finally conclude that

$\displaystyle I[s^{K};\tilde{y}^{K}]$ $\displaystyle\leq
h[s^{K}_{\mathcal{N}}]+h[\tilde{y}^{K}_{\mathcal{N}}]-h[s^{K}_{\mathcal{N}},\tilde{y}^{K}_{\mathcal{N}}]$
$\displaystyle+\frac{n}{2}+\frac{n}{2}\log(2\pi)+\frac{n}{2}\log c(n),$ (87)

which is equivalent to (68). $\blacksquare$

###  X-D Proof of Proposition 2

Mutual information $I\left[s^{K};\tilde{y}^{K}\right]$ can be written in terms
of differential entropies as
$I[s^{K};\tilde{y}^{K}]=h[s^{K}]+h[\tilde{y}^{K}]-h[s^{K},\tilde{y}^{K}]$.
Besides, it has been proved in Appendix C that $y^{K}$, $s^{K}$, and
$\tilde{y}^{K}$, follow log-concave distributions. Then, from the inequality
on the entropy of log-concave distributed random vectors, we can conclude that
(see Theorem $I.1$ in [48]):

$\left\\{\begin{aligned}
\frac{1}{Kn_{s}}h[s^{K}]&\leq\frac{1}{Kn_{s}}h[s^{K}_{{\mathcal{N}}^{*}}]+\frac{1}{2},\\\
\frac{1}{Kn_{y}}h[\tilde{y}^{K}]&\leq\frac{1}{Kn_{y}}h[\tilde{y}^{K}_{{\mathcal{N}}^{*}}]+\frac{1}{2},\\\
\frac{1}{n}h[s^{K},\tilde{y}^{K}]&\geq\frac{1}{n}h[s^{K}_{{\mathcal{N}}^{*}},\tilde{y}^{K}_{{\mathcal{N}}^{*}}]-\frac{1}{2}.\end{aligned}\right.$
(88)

Combining all the inequalities in (88), we have

$\displaystyle\frac{1}{n}I[s^{K};\tilde{y}^{K}]$
$\displaystyle\leq\frac{1}{n}(h[s^{K}_{{\mathcal{N}}^{*}}]+h[\tilde{y}^{K}_{{\mathcal{N}}^{*}}]-h[s^{K}_{{\mathcal{N}}^{*}},\tilde{y}^{K}_{{\mathcal{N}}^{*}}])+1,$
(89)

which is equivalent to (69). $\blacksquare$

###  X-E Proof of Lemma 5

Due to the optimization problem (21), and relations (22) and (52b), the
optimal information leakage in time window $\mathcal{K}$, which is equivalent
to optimal objective function in optimization problem (21), can be written as
follows:

$\displaystyle\gamma_{\mathcal{K}}$
$\displaystyle=\frac{1}{2}\log{\det\left({\Sigma_{K}^{s}}\right)}$ (90)
$\displaystyle-\frac{1}{2}\log{\det{\left({\Sigma}_{K}^{s}-{{\Sigma}_{K}^{sy^{\prime}}}^{\top}{{{\Sigma}_{K}^{y^{\prime}}}^{*}}^{-1}{\Sigma}_{K}^{sy^{\prime}}\right)}}$
$\displaystyle-\frac{1}{2}\log{\det\left({\Sigma_{K}^{j}}^{*}\right)}-\frac{Kn_{y}}{2}-\frac{Kn_{y}}{2}\log(2\pi),$

where ${\Sigma_{K}^{j}}^{*}$ is the optimal ${\Sigma_{K}^{j}}$ and
${{{\Sigma}_{K}^{y^{\prime}}}^{*}}={\tilde{C}_{K}Q\tilde{C}_{K}^{\top}+(I_{K}\otimes\Sigma^{w})+{{\Sigma}^{v}_{K}}^{*}}$
with ${{\Sigma}^{v}_{K}}^{*}$ as the optimal ${{\Sigma}^{v}_{K}}$. Similarly,
for the optimal information leakage in subwindows of time $\mathcal{K}_{i}$
considering $i\in{1,2,\ldots,N}$ one can write:

$\displaystyle\gamma_{\mathcal{K}_{i}}$
$\displaystyle=\frac{1}{2}\log{\det\left({\Sigma_{K_{i}}^{s}}\right)}$ (91)
$\displaystyle-\frac{1}{2}\log{\det{\left({\Sigma}_{K_{i}}^{s}-{{\Sigma}_{K_{i}}^{sy^{\prime}}}^{\top}{{{\Sigma}_{K_{i}}^{y^{\prime}}}^{*}}^{-1}{\Sigma}_{K_{i}}^{sy^{\prime}}\right)}}$
$\displaystyle-\frac{1}{2}\log{\det\left({\Sigma_{K_{i}}^{j}}^{*}\right)}-\frac{{K_{i}}n_{y}}{2}-\frac{{K_{i}}n_{y}}{2}\log(2\pi),$

where ${\Sigma_{K_{i}}^{j}}^{*}$ is the optimal ${\Sigma_{K_{i}}^{j}}$ and
${{{\Sigma}_{K_{i}}^{y^{\prime}}}^{*}}={\tilde{C}_{K}Q\tilde{C}_{K}^{\top}+(I_{K_{i}}\otimes\Sigma^{w})+{{\Sigma}^{v}_{K_{i}}}^{*}}$
with ${{\Sigma}^{v}_{K_{i}}}^{*}$ as the optimal ${{\Sigma}^{v}_{K_{i}}}$ in
subwindows of time. In the objective function (90), the logarithm function is
monotonically increasing. Therefore, to find the optimal objective function
(information leakage), we need to maximize ${{\Sigma}^{v}_{K}}$ and
${\Sigma_{K}^{j}}$. Then, ${{\Sigma}^{v}_{K}}^{*}$ and ${\Sigma_{K}^{j}}^{*}$
are the maximum ${{\Sigma}^{v}_{K}}$ and ${\Sigma_{K}^{j}}$, respectively.
Correspondingly, ${{\Sigma}^{v}_{K_{i}}}^{*}$ and ${\Sigma_{K_{i}}^{j}}^{*}$
are the maximum ${{\Sigma}^{v}_{K_{i}}}$ and ${\Sigma_{K_{i}}^{j}}$. As a
result:

${{\Sigma}^{v}_{K}}^{*}\geq{\tilde{\Sigma}^{v^{*}}_{K}}:=\operatorname{diag}({{\Sigma}^{v}_{K_{1}}}^{*},{{\Sigma}^{v}_{K_{2}}}^{*},\ldots,{{\Sigma}^{v}_{K_{N}}}^{*}).$
(92)

Because in ${\tilde{\Sigma}^{v^{*}}_{K}}$, the non-diagonal block terms are
zero, we only maximize its block diagonal terms, while in
${{\Sigma}^{v}_{K}}^{*}$, the ${{\Sigma}^{v}_{K}}$ is maximized by designing
all its blocks. Similarly,

${{\Sigma}^{j}_{K}}^{*}\geq{\tilde{\Sigma}^{j^{*}}_{K}}:=\operatorname{diag}({{\Sigma}^{j}_{K_{1}}}^{*},{{\Sigma}^{j}_{K_{2}}}^{*},\ldots,{{\Sigma}^{j}_{K_{N}}}^{*}).$
(93)

From (92) one can write:

$\displaystyle\log{\det{\left({\Sigma}_{K}^{s}-{{\Sigma}_{K}^{sy^{\prime}}}({\Sigma_{K}^{y^{\prime}_{1}}+{{\Sigma}^{v}_{K}}^{*}})^{-1}{{\Sigma}_{K}^{sy^{\prime}}}^{\top}\right)}}\geq$
(94)
$\displaystyle\log{\det{\left({\Sigma}_{K}^{s}-{{\Sigma}_{K}^{sy^{\prime}}}({\Sigma_{K}^{y^{\prime}_{1}}+{\tilde{\Sigma}^{v^{*}}_{K}}})^{-1}{{\Sigma}_{K}^{sy^{\prime}}}^{\top}\right)}},$

where
$\Sigma_{K}^{y^{\prime}_{1}}=\tilde{C}_{K}Q\tilde{C}_{K}^{\top}+(I_{K}\otimes\Sigma^{w})$.
Then, based on (94):

$\displaystyle\log{\det{\left({\Sigma}_{K}^{s}-{{\Sigma}_{K}^{sy^{\prime}}}({\Sigma_{K}^{y^{\prime}_{1}}+{{\Sigma}^{v}_{K}}^{*}})^{-1}{{\Sigma}_{K}^{sy^{\prime}}}^{\top}\right)}}\geq$
(95)
$\displaystyle\log{\det{\left({\Sigma}_{K}^{s}-{{\Sigma}_{K}^{sy^{\prime}}}({\Sigma_{K_{1}}^{y^{\prime}_{1}}+{{\Sigma}^{v}_{K_{1}}}^{*}})^{-1}{{\Sigma}_{K}^{sy^{\prime}}}^{\top}\right)}}+$
$\displaystyle\ldots+\log{\det{\left({\Sigma}_{K}^{s}-{{\Sigma}_{K}^{sy^{\prime}}}({\Sigma_{K_{N}}^{y^{\prime}_{1}}+{{\Sigma}^{v}_{K_{N}}}^{*}})^{-1}{{\Sigma}_{K}^{sy^{\prime}}}^{\top}\right)}}.$

Furthermore, from (93) and the fact that:

$\displaystyle\log{\det\left({\tilde{\Sigma}^{j^{*}}_{K}}\right)}$
$\displaystyle=\log{\det\left({{\Sigma}^{j}_{K_{1}}}^{*}\right)}+\ldots+\log{\det\left({{\Sigma}^{j}_{K_{N}}}^{*}\right)},$
(96)

we can conclude that:

$\displaystyle\log{\det\left({{\Sigma}^{j}_{K}}^{*}\right)}$
$\displaystyle\geq\log{\det\left({{\Sigma}^{j}_{K_{1}}}^{*}\right)}+\ldots+\log{\det\left({{\Sigma}^{j}_{K_{N}}}^{*}\right)}.$
(97)

Therefore, combining (90), (91), (95) and (97), we can conclude (70).
$\blacksquare$

###  X-F Proof of $\tilde{z}_{k}$ approximately following a gamma distribution

At first, we prove that defining $m_{k}$ as a standard normal variable
$m_{k}\sim\mathcal{N}(\mathbf{0},I_{n})$ and matrix $Q\in\mathbb{R}^{n\times
n}$ with $n\in\mathbb{R}^{+}$, $Q^{\prime}=m_{k}^{\top}Qm_{k}$ approximately
follows a gamma distribution. $Q^{\prime}=m_{k}^{\top}Qm_{k}$ is a weighted
sum of chi-squared distributions with one degree of freedom
($\mathscr{X}^{2}_{1}$) as
$Q^{\prime}\sim\sum\limits_{i=1}^{{n}}\lambda_{i}\mathscr{X}^{2}_{1}$ [50],
where $\lambda_{i}$, $i\in\\{1,\ldots,n\\}$, are the eigenvalues of the matrix
$Q$. For each term of the sum $\lambda_{i}\mathscr{X}^{2}_{1}$ we have [44]:

$\lambda_{i}\mathscr{X}^{2}_{1}\sim\Gamma(K_{i}=\frac{1}{2},\theta_{i}=2\lambda_{i}).$
(98)

Therefore,
$Q^{\prime}\sim\sum\limits_{i=1}^{{n}}\Gamma(K_{i}=\frac{1}{2},\theta_{i}=2\lambda_{i})$
is the sum of gamma variables. Based on the Welch-Satterthwaite Approximation
[60, 61], the summation of gamma distributions
$\sum\limits_{i=1}^{{n}}\Gamma(K_{i},\theta_{i})$ approximately follows a
gamma distribution with shape
$K_{sum}=\frac{\left(\sum_{i}\theta_{i}K_{i}\right)^{2}}{\sum_{i}\theta_{i}^{2}K_{i}}$
and scale $\theta_{sum}=\frac{\sum\theta_{i}K_{i}}{K_{sum}}$. Consequently,
$Q^{\prime}$ follows a gamma distribution $Q^{\prime}\sim\Gamma(K,\theta)$
with shape
$K=\frac{\left(\sum_{i=1}^{n}\lambda_{i}\right)^{2}}{2\sum_{i=1}^{n}\left(\lambda_{i}\right)^{2}}$
and scale
$\theta=2\frac{\sum_{i=1}^{n}\left(\lambda_{i}\right)^{2}}{\left(\sum_{i=1}^{n}\lambda_{i}\right)}$.
The accuracy of this approximation is proved in [50].

Hence, we can conclude that
$\tilde{z}_{k}=m_{k}^{\top}\Sigma^{\prime}_{k}m_{k}$ approximately follows a
gamma distribution $\tilde{z}_{k}\sim\Gamma(K_{z}^{k},\theta_{z}^{k})$ with:

$\displaystyle
K_{z}^{k}=\frac{\left(\sum_{i=1}^{n_{y}}\lambda_{i}^{k}\right)^{2}}{2\sum_{i=1}^{n_{y}}\left(\lambda_{i}^{k}\right)^{2}},$
(99)
$\displaystyle\theta_{z}^{k}=2\frac{\sum_{i=1}^{n_{y}}\left(\lambda_{i}^{k}\right)^{2}}{\left(\sum_{i=1}^{n_{y}}\lambda_{i}^{k}\right)},$
(100)

where $\lambda_{i}^{k}$, $i\in\\{1,\ldots,n_{y}\\}$ are the eigenvalues of the
matrix
$\Sigma^{\prime}_{k}=(\tilde{\Sigma}^{1/2}_{k})^{\top}\Sigma^{-1}\tilde{\Sigma}^{1/2}_{k}$.
$\blacksquare$