A Two-Stage Architecture for Differentially Private Kalman Filtering and LQG Control
Kwassi H. Degue Jerome Le Ny
cs.SY cs.CR









 




































Degue et al.: Differentially Private Kalman Filtering and LQG Control





















A Two-Stage Architecture for Differentially Private Kalman Filtering and LQG Control
    Kwassi H. Degue
        and Jerome Le Ny, Senior Member, IEEE
A preliminary version of this paper appeared in <cit.>. 
This work was supported by NSERC under Grant RGPIN-5287-2018 and RGPAS-2018-522686, 
by the Pierre Arbour Foundation doctoral scholarship and by an FRQNT doctoral scholarship.

The authors are with the Department of Electrical Engineering,
Polytechnique Montreal and GERAD, Montreal, QC H3T-1J4, Canada 
{kwassi-holali.degue, jerome.le-ny}@polymtl.ca.
    December 30, 2023
==========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================












Large-scale monitoring and control systems enabling a more intelligent infrastructure 
increasingly rely on sensitive data obtained from private agents, e.g., location 
traces collected from the users of an intelligent transportation system. 
In order to encourage the participation of these agents, it becomes then critical 
to design algorithms that process information in a privacy-preserving way. 
This article revisits the Kalman filtering and Linear Quadratic Gaussian (LQG) 
control problems, subject to privacy constraints. 

We aim to enforce differential privacy, a formal, state-of-the-art 
definition of privacy ensuring that the output of an algorithm is not too 
sensitive to the data collected from any single participating agent.
A two-stage architecture is proposed that first aggregates and combines the individual 
agent signals before adding privacy-preserving noise and post-filtering the result to 
be published. We show a significant performance improvement offered by this architecture 
over input perturbation schemes as the number of input signals increases and that 
an optimal static aggregation stage can be computed by solving a semidefinite program.
The two-stage architecture, which we develop first for Kalman filtering, is then 
adapted to the LQG control problem by leveraging the separation principle.
Numerical simulations illustrate the performance improvements over differentially 
private algorithms without first-stage signal aggregation. 




Differential privacy; Kalman filtering; Estimation; Filtering; 
LQG control; Optimal control

















§ INTRODUCTION


To monitor and control intelligent infrastructure systems such as smart grids,
smart buildings or smart cities, data needs to be continuously collected 
from the people interacting with these systems, either through sensors installed 
in the environment such as cameras and smart meters, or through personal devices 
such as smartphones. 
Hence, in contrast to more traditional control systems, the measured signals for
such systems often contain highly privacy-sensitive information, e.g., related 
to the real-time location or health of a person.
For example, the accuracy of crowd-sourced traffic maps and congestion-aware routing 
applications is increased by using data provided by smartphones and connected vehicles
<cit.>.
However, individual location data turns out to be very difficult to properly 
anonymize because individuals have highly unique mobility patterns 
<cit.>, 
and in fact individual trajectories can be reconstructed even from 
just aggregate location data <cit.>.
Similarly, fine-grained measurements of a house's electric power consumption 
collected by a smart meter can enable demand-response schemes, but can 
also be used to infer the activities of the occupants, 
by identifying the usage of individual appliances 
<cit.>.
Therefore, it is necessary to implement privacy-preserving mechanisms 
when sensitive data must be shared to improve a system's performance.


Various definitions of privacy have been proposed that are amenable to formal analysis.
While a survey of such definitions is out of the scope of this paper, we can mention 
some recent work focusing on signal processing and control problems.
Privacy is measured by a lower bound on the mutual information between published
and private signals in <cit.>, 
on the Fisher information in <cit.>, or on the error covariance 
of the estimator of a sensitive signal 
in <cit.>.
The concept of k-anonymity and its extensions has been applied to the publication 
of location traces in <cit.>.
But much of the recent research on privacy-preserving data analysis relies on the notion
of differential privacy <cit.>. 
In the standard set-up, which is also the situation considered in this article, a data
holder aims to release the results of computations based on private data. Differential
privacy is enforced by adding an appropriate amount of noise to the published results, 
in such a way that the probability distribution over the outputs does not depend too 
much on the data of any single individual.
As a result, the ability of a third party observing the outputs to make 
new inferences about a given person is roughly the same, whether or not that
person chooses to contribute its data. This guarantee can then be used to 
weigh the risks of information disclosure against the benefits of publishing 
more accurate analyses.


A large number of techniques have been developed to compute differentially private
versions of various statistics, see <cit.> for an overview.

Nevertheless, the differentially private analysis of streaming data remains
relatively less explored <cit.>, despite its 
importance for signal processing and control applications. 
Some previous work has focused on the design of differentially private dynamic estimators <cit.>, 
controllers <cit.>, 
consensus algorithms <cit.>,
or anomaly detectors 
<cit.>.
In particular, <cit.> discusses the Kalman filtering problem under 
a differential privacy constraint and compares schemes introducing noise either 
directly on the measured signals (input perturbation mechanisms) or on the 
published estimate (output perturbation mechanisms). Output perturbation provides
better performance as the number of input signals increases, but has the
drawback of leaving unfiltered noise on the output, which motivates the 
two-stage architecture that we consider here.
In <cit.>, the authors consider a multi-agent linear quadratic tracking 
problem where the trajectory tracked by each agent should remain private, while
<cit.> considers an LQG control problem where each agent wishes 
to keep its individual state private. In both cases, noise is added directly
on the individual measurements, a form of input perturbation.


In this paper we study the design of Kalman filters and LQG controllers subject 
to a differential privacy constraint on the measured signals. 
These problems, stated formally in Section <ref>, arise when a data
collector measures private signals originating from a population of agents, whose 
dynamics can be modeled as linear Gaussian systems, in order to publish in real-time 
either an estimate of an aggregate state of the agent population, or a control signal 
shared by the agents and aimed at regulating such an aggregate state.
As a motivating example, one can consider the problem of controlling the distribution
of vehicles on a road network by means of traffic messages broadcasted to all cars, with
the current density estimated from location data obtained from the smartphones of 
individual drivers.
 

Section <ref> presents the main contribution of this paper,
namely, a two-stage architecture for differentially private Kalman filtering, 
where the privacy-preserving noise is added only after an input stage appropriately 
combining the measured signals of the individual agents, while an output stage filters 
out this noise.
Such two-stage architectures were discussed in <cit.> but have not yet
been applied to the Kalman filtering problem, and we argue first in 
Section <ref> via a simple example that significant performance 
improvements can be expected compared to input perturbation mechanisms. 
We show that the optimal input stage can be computed by solving a semidefinite
program (SDP), hence, a tractable convex optimization problem. 
The fact that the input stage design problem admits an SDP formulation 
is reminiscent of other Kalman filtering problems subject to resource 
or communication constraints, see, e.g., 
<cit.>,
but the SDP capturing specifically the differential privacy constraint is new.
The design procedure is then adapted in Section <ref> 
to the LQG control problem. By exploiting the classical properties of the optimal
LQG controller (linearity and separation principle), we can view the control problem 
as the problem of estimating a certain linear combination of the agent states as 
in Section <ref>, but for a specific cost on the estimation error.


A conference version of this paper appeared in <cit.>, but contained no proof
and did not discuss the LQG control problem. 
The numerical examples discussed in Sections <ref> and <ref> 

are also new.
Finally, we introduce some notation used throughout this paper. 
We fix a generic probability triple (Ω,ℱ,ℙ), where ℱ 
is a σ-algebra on Ω and ℙ a probability measure
defined on ℱ. 
The notation X ∼𝒩(μ,Σ) means that X is a Gaussian 
random vector with mean vector μ and covariance matrix Σ.
“Independent and identically distributed” is abbreviated iid.
We denote the p-norm of a vector x∈ℝ^k by 
|x| _p:=(∑_i=1^k|x_i|^p)^1/p, for p ∈ [1,∞). 
For a matrix A, the induced 2-norm (maximum singular value of A) is denoted 
A_2 and the Frobenius norm A_F := √(Tr(A^TA)).
If A and B are symmetric matrices, A ≽ B (resp. A ≻ B) means that 
A-B is positive semi-definite (resp. positive definite).
We use the notation diag(A_1,…,A_n) to represent a block-diagonal
matrix with the matrices A_i on the diagonal. The column vector of size n with all
components equal to 1 is denoted 1_n.
Finally, for a discrete-time signal x we denote x_0:t := {x_0,…,x_t}.





§ PROBLEM STATEMENT
 
 



 §.§ Privacy-Preserving State Estimation and LQG Control for a Population of Dynamic Agents



Consider a set of n privacy-sensitive signals {y_i,t}_0 ≤ t ≤ T, 
i=1,…,n, with y_i,t∈ℝ^p_i, collected by a data aggregator, 
and which could originate from n distinct agents. Let p = ∑_i=1^n p_i.
We assume that a mathematical model capturing known dynamic and statistical properties
of these signals is publicly available, consisting of a linear system with n independent 
(vector-valued) states associated to the n measured signals

    x_i,t+1   = A_i,t  x_i,t +B_i,t  u_t+ w_i,t,    0 ≤ t ≤ T-1, 
    
    y_i,t   = C_i,t  x_i,t + v_i,t,    0 ≤ t ≤ T,

for i=1,…,n, where x_i,t, w_i,t∈ℝ^m_i.
Here w_i,t∼𝒩(0,W_i,t) and v_i,t∼𝒩(0,V_i) are 
independent sequences of iid zero-mean Gaussian random vectors with covariance matrices 
W_i,t≻ 0, V_i≻ 0, for i=1, …, n. 
In particular, assuming that the matrices W_i,t are invertible is necessary in the 
following to be able to use the “information filter” form of the Kalman filter 
equations <cit.>.
The sequence u with u_t∈ℝ^h represents a control input that is shared 
by the n individuals. This is motivated by scenarios in which a common signal is broadcast
to drive the aggregate state of a population, while individual signals can still be subject to
privacy constraints. 
The initial conditions x_i,0 are independent Gaussian random vectors that are also 
independent of the noise processes w and v, with mean x_i,0 
and covariance matrices Σ_i,0^-≻ 0. 
Let x_t := [x_1,t^T,…,x_n,t^T]^T, y_t := [y_1,t^T,…,y_n,t^T]^T, 
w_t = [w_1,t^T,…,w_n,t^T]^T and  v_t = [v_1,t^T,…,v_n,t^T]^T denote 
the global state, measurement and noise signals of (<ref>). 
Define A_t := (A_1,t,…,A_n,t), B_t = [B_1,t^T,…,B_n,t^T]^T, C_t = (C_1,t,…,C_n,t),
W_t = (W_1,t,…,W_n,t), V = (V_1,…,V_n). Then the system (<ref>)
can be rewritten more compactly as

    x_t+1    = A_t   x_t+ B_t   u_t + w_t,	    0 ≤ t ≤ T-1,	
    
    y_t    = C_t   x_t+v_t,    0 ≤ t ≤ T,

with w_t ∼𝒩(0,W_t) and v_t ∼𝒩(0,V). Throughout the paper,
the model parameters x_i,0, Σ^-_i,0, A_i,t, B_i,t, C_i,t, 
W_i,t, V_i, are assumed to be publicly known information.

In Section <ref>, we first consider a filtering problem (with u a known signal) 
where the data aggregator aims to publish at each period t a causal estimate ẑ_t of 
a linear combination z_t= L_t x_t = ∑_i=1^nL_i,tx_i,t of the individual states, 
computed from the signals y_i, with L_t := [ L_i,t, …, L_n,t ] 
some given (publicly known) matrices. 
This estimator should minimize the Mean Square Error (MSE) performance measure

    E_T := 1/T+1∑_t=0^T 𝔼[z_t - ẑ_t_2^2].

For privacy reasons, the signals y_i are not released by the data aggregator and moreover 
the publicly released estimate ẑ should also guarantee the differential privacy of the input 
signals y_i, as defined precisely in Section <ref>.
For example, the signals y_i could represent position measurements of n individuals, 
each state x_i could consist of the position and velocity of individual i, 
and the goal might be to publish only a real-time estimate of the average velocity 
of all individuals.
Note that in the absence of privacy constraint, 
the optimal estimator is ẑ_t=∑_i=1^nL_i,tx̂_i,t, with x̂_i,t 
provided by the (time-varying) Kalman filter estimating the state x_i of subsystem i from the 
signal y_i <cit.>, and in particular the estimation problem then decouples for the 
n subsystems. 


Next, in Section <ref> we build on the results obtained for
the filtering problem to study the following privacy-constrained Linear Quadratic Gaussian (LQG) 
regulation problem. The data aggregator uses the measured
signals y_i, 1 ≤ i ≤ n, to compute and broadcast a (causal) control signal u that minimizes 
the following quadratic cost for the n agents

    J_T= 1/T+1𝔼[ ∑_t=0^T-1(x_t^T Q_t x_t+u_t^T R_t u_t) + x_T^T Q_T x_T],

where Q_t ≽ 0 for 0 ≤ t ≤ T and R_t ≻ 0 for 0 ≤ t ≤ T-1 are
publicly known weight matrices.
Again, only the signal u is published (in particular, it is available to the n agents), 
and releasing u must guarantee the differential privacy of the measured signals y_i.
It is worth noting that the cost function (<ref>) can be used to
drive an aggregate value of the global population state toward 0 rather than the
individual agent states, since trying to do the latter might be in direct conflict with
the privacy requirement (which, essentially, aims to hide the value of the individual
signals y_i, and hence indirectly of the individual states x_i).
For example, we can have 
x_t Q_t x_t =  ( 1/n∑_i=1^n x_i,t)^2 if Q_t = 1/n^21_n 1_n^T,
in order to regulate the average population state.
Figure <ref> represents the estimation and control problem setups, including a basic
scheme to enforce differential privacy by injecting noise directly in the signals {y_i}_1 ≤ i ≤ n,
as described in Section <ref>.



Note that the steady-state versions of both the filtering and LQG control problems are also considered, with the 
model (<ref>)-(<ref>) in this case assumed time-invariant and the performance 
measures defined as

    E_∞ := lim_T →∞ E_T,    J_∞ := lim_T →∞ J_T.





To finish stating the above problems formally, we define in the next section the differential privacy 
constraint imposed on the signal ẑ for the filtering problem and u for the LQG problem.



 §.§ Differential Privacy Constraint

 


Differential privacy <cit.> is a property satisfied by certain randomized algorithms (also called 
mechanisms), which in an abstract setting compute outputs in a space 𝖱 based on sensitive 
data in a space 𝖣. To be differentially private, the algorithm must ensure that the probability
distribution of its randomized output is not very sensitive to certain variations in the input data, which
are specified as part of the privacy requirement.

More concretely, we equip the input space 𝖣 with a symmetric binary relation called
adjacency and denoted Adj, which captures the variations in the input datasets that we want to
make hard to detect by observing the outputs.
In our case, the input space 𝖣 is the vector space ℝ^p(T+1) of global measurement 
signals y, and a mechanism is a causal stochastic system producing an output signal (ẑ or u 
in the previous section) based on its input y.

We define two measured signals to be adjacent if the following condition holds

    Adj(y,y') iff for some  1 ≤ i ≤ n,  y_i-y_i' _2≤ρ_i, 
    and y_j=y_j' for all j≠ i,

with {ρ_i}_i=1^n∈ℝ_+^n a given set of positive numbers
and with the definition of the ℓ_2-norm v_2 := (∑_t=0^T |v_t|_2^2)^1/2 for
a vector-valued signal v.
In other words, two adjacent measurement signals can differ by the values of 
a single participant, with only ℓ_2-bounded signal deviations allowed for each individual.
For any two inputs that satisfy the adjacency relation, the following definition characterizes 
the deviation that is allowed for a differentially private mechanism's output distribution. 

 
Let 𝖣 be a space equipped with a symmetric binary relation denoted Adj 
and let (𝖱,) be a measurable space, where 
 is a given σ-algebra over 𝖱. 
Let ϵ≥ 0, 1 ≥δ≥ 0. A randomized mechanism M from 
𝖣 to 𝖱 is (ϵ,δ)-differentially private (for Adj) if 
for all d,d'∈𝖣 such that Adj(d,d'), 

    ℙ(M(d)∈ S) ≤ e^ϵ ℙ(M(d')∈ S) + δ, 
     ∀ S ∈ℛ.

  

In Definition <ref>, smaller values of ϵ and δ correspond 
to stronger privacy guarantees, i.e., distributions for M(d) and M(d') that are closer in
(<ref>). Next, we need tools that can be used to enforce the property of 
Definition <ref>. The following definition is useful for mechanisms 
such as ours that produce outputs in vector spaces.

	
Let 𝖣 be a space equipped with an adjacency relation Adj.  
Let 𝖱 be a vector space equipped with a norm ·_𝖱. 
The sensitivity of a mapping q:𝖣↦𝖱 is defined as

    q := sup_{d,d':Adj(d,d')} q(d)-q(d') _𝖱.

For 𝖱=ℝ^h(T+1) or (ℝ^h)^ℕ equipped with 
the ℓ_2-norm, this defines the ℓ_2-sensitivity of q, denoted _2q. 


The Gaussian mechanism <cit.> consists in adding Gaussian noise 
proportional to the ℓ_2-sensitivity of a mapping to enforce (ϵ,δ)-differential privacy. 
A fairly tight upper bound on the proportionality constant is provided in <cit.>.
Recall first the definition of the 𝒬-function
𝒬(x):=1/√(2π)∫_x^∞exp(-u^2/2)du, which is monotonically 
decreasing from (-∞,∞) to (0,1).
Then, for 1 > δ > 0, define
κ_δ,ϵ := 1/2ϵ( Q^-1(δ)+√((Q^-1(δ))^2+2ϵ)).
The following theorem can be found in <cit.>.

	
Let ϵ > 0, 1 > δ > 0. Let G be a dynamic system with p inputs and q outputs. 
Then the mechanism M(y) = Gy + ν, where ν is a white Gaussian noise (sequence of iid zero-mean 
Gaussian vectors) with ν_t ∼𝒩(0,κ_δ,ϵ^2 (Δ_2 G)^2 I_q),
is (ϵ,δ)-differentially private.


In other words, Theorem <ref> says that we can produce a differentially
private signal by adding white Gaussian noise at the output of a system G processing the sensitive
signal y, with covariance matrix σ^2 I_q, and σ proportional to the ℓ_2-sensitivity of G.




 §.§ A First Solution: Input Perturbation Mechanism



An important property of differential privacy is its resilience to 
post-processing <cit.>, i.e., applying further computations to an output that is 
differentially private does not degrade the differential privacy guarantee, as long as the original
dataset is not reaccessed for these computations. 
This leads immediately to a first solution for the estimation and control problems stated in the previous 
section, called the input perturbation mechanism, which consists in perturbing each measured
signal y_i directly to release differentially private versions of these signals. 

First, note that the memoryless system defined by (G y)_t = M y_t, where M is the diagonal matrix
M = diag (1/ρ_1, …, 1/ρ_n), has the sensibility bound

    Δ G := sup_Adj(y,y')My - My'_2 ≤ 1

for the adjacency relation (<ref>). Hence, by Theorem <ref>, 
releasing the signals {y_i/ρ_i + ν_i}_1 ≤ i ≤ n, where each signal ν_i is 
a white Gaussian noise with covariance matrix κ_δ,ϵ^2 I_p_i, 
is (ϵ,δ)-differentially private. Equivalently, using the resilience to 
post-processing to multiply this output by M^-1, we see that releasing the signals 
{ỹ_i := y_i + ρ_i ν_i}_1 ≤ i ≤ n is (ϵ,δ)-differentially private. 
Once these signals are released, applying further processing on them does not impact the 
differential privacy guarantee. Moreover, these signals are of the same form as the outputs of system (<ref>), except for a higher level of (still Gaussian) noise due to the addition 
of the artificial privacy-preserving noise. One can therefore produce an (ϵ,δ)-differentially 
private estimate ẑ or control signal u discussed in Section <ref> 
by applying standard Kalman filtering and LQG design techniques to the signals ỹ_i,
see Figure <ref>.
An advantage of input perturbation mechanisms is that each agent can release directly the 
differentially private signal ỹ_i, and hence does not need to trust the data aggregator
to enforce the differential privacy property. Moreover, this mechanism has the potentially useful
feature of publishing the individual signals ỹ_i, which could be used for other
purposes than the original estimation or control problem.
Nonetheless, as we discuss in the following sections, input perturbation typically leads 
to a high level of noise and hence performance degradation, which motivates the search 
for better mechanisms.









§ DIFFERENTIALLY PRIVATE KALMAN FILTERING
 






Input perturbation for the differentially private Kalman filtering problem, as discussed 
in Section <ref> and represented on Figure <ref>, 
was considered in <cit.>. 
Here, we show first in Section <ref> via a simple example
that the performance of this mechanism can be significantly improved by combining 
the individual input signals before adding the privacy preserving noise.
This leads to the two-stage mechanism of Figure <ref>,
whose systematic design is discussed in Section <ref>.





 §.§ A Scalar Example
  

Consider a scalar homogeneous version of model (<ref>) with 
A_i,t=a ∈ℝ, B_i,t = b, C_i,t=c, W_i,t=σ_w^2, V_i=σ_v^2 
and L_i,t=1 (so z_t = ∑_i=1^n x_i,t), 
and assume ρ_i = ρ 
for all i in (<ref>).

In other words, A_t = a I_n, B_t = b 1_n, C_t = c I_n, 
W_t = σ_w^2 I_n, V = σ_v^2 I_n in (<ref>).

We also let T →∞ in the problem statement and consider the steady-state MSE
E_∞, see (<ref>),
as performance measure for a given estimator ẑ of z. Moreover in this section 
we consider for simplicity the minimum mean square error (MMSE) estimate 
ẑ_t of x_t given the measurements up to time t-1 only, since the corresponding 
MSE is directly obtained by solving an algebraic Riccati equation (ARE).


Let α := κ_δ,ϵ ρ. Since 
sup_y,y':Adj(y,y')‖ y-y'‖ _2 = ρ, 
by Theorem <ref> and as explained in Section <ref>,
releasing the signal r_t=  y_t+ν_t, with white noise ν such that ν_t∼𝒩(0,α^2I_n), 
is (ϵ,δ)-differentially private for the adjacency relation (<ref>). 
The steady-state MSE of a Kalman filter estimating z for the system with dynamics as in (<ref>)
and measurements r is obtained by solving a 
scalar ARE, which leads to the following expression 

    E_∞^1 =
    n/2c^2(-β+√(β^2+4(α^2
    +σ_v^2)σ_w^2c^2)),

where β = (1-a^2) (σ_v^2+α^2)-c^2 σ_w^2.


Instead of input perturbation, we can use the architecture shown on 
Fig. <ref> with D=1_n^T, 
a 1 × n row vector of ones.
Consider the same adjacency relation (<ref>) and denote 

η_t = D y_t = ∑_i=1^ny_i,t, θ_t=∑_i=1^nw_i,t 
and λ_t=∑_i=1^nv_i,t.
We have

    z_t+1    = a z_t+θ_t,
    η_t    = c z_t+λ_t.

Since again sup_y,y':Adj(y,y')‖ Dy-Dy'‖ _2 = ρ,
releasing the scalar signal s_t= η_t + ζ_t, with ζ_t∼𝒩(0,α^2), 
is (ϵ,δ)-differentially private for the adjacency relation (<ref>). 
The MSE of a Kalman filter estimating z from this signal s, with the dynamics of the model (<ref>),
can again be obtained by solving an ARE, which leads to the following expression

    E_∞^2 =
    n/2c^2(-β_(n)+√(β_(n)^2+4(α^2/n +σ_v^2) σ_w^2 c^2)),

where β_(n) = (1-a^2) (σ_v^2+α^2/n)-c^2 σ_w^2.

Comparing (<ref>) and (<ref>), we see that the only difference is
the vanishing influence of the privacy preserving noise on E_∞^2 as n increases, with
the term α^2/n replacing α^2 in E_∞^1.
For example, if n=100, a=1, c=1, ρ=50, ϵ=ln(3), δ=0.05, σ_w^2=0.5, 
σ_v^2=0.9, we obtain E_∞^1 ≈ 6235 and E_∞^2 ≈ 650.
It is indeed desirable that as the number of agents n increases, differential privacy becomes easier to enforce
and the impact of the privacy requirement on achievable performance decreases, a feature that the architecture 
of Figure <ref> has the potential to achieve. The design of this architecture is discussed 
in the next section for the general filtering problem of Section <ref>. 






 §.§ Design of the Two-Stage Mechanism
 


Following Figure <ref>, we construct a differentially private 
estimate ẑ_t of z_t by first multiplying the global signal y with 
a constant matrix 

    D=[ D_1 D_2   … D_n ],

with the matrices D_i ∈ℝ^q× p_i, 1 ≤ i ≤ n, to be designed and q to be determined.
Then, we add white Gaussian noise ζ according to the Gaussian mechanism,
in order to make the signal s differentially private, with

    s_t = D y_t + ζ_t = 
    D C_t x_t + D v_t + ζ_t,   0 ≤ t ≤ T.

Therefore, the role of the matrix D is to combine the individual signals appropriately
before adding the privacy-preserving noise, in order to decrease the overall sensitivity
(see Definition <ref>), while preserving enough information for z to
be estimated with sufficient accuracy.
Finally, we construct a causal MMSE estimator ẑ 
of z from s, a task for which it is optimal to use a Kalman filter, since the system 
model producing s with the state dynamics of (<ref>) is still linear and Gaussian. 
This Kalman filter produces a state estimate x̂ of x and then ẑ_t = L_t x̂_t 
for all t.



Given D, for measurement signals y and y' adjacent according to (<ref>) 
and differing at index i, we have

    Dy-Dy' _2 =  D_iy_i-D_iy_i'_2≤ρ_i D_i_2,

where D_i_2 denotes the maximum singular value of the matrix D_i,
and there are adjacent signals y_i, y_i' achieving the bound.
Hence, we can bound the sensitivity of the memoryless system y ↦ Dy as follows 

    _2D     :=  sup_y,y':Adj(y,y')‖ Dy-Dy'‖ _2
        =  1 ≤ i ≤ nmax{ρ_i   D_i_2 }.
 
Therefore, from Theorem <ref>, for any matrix D, releasing 
s_t=D y_t+ζ_t, with ζ_t∼𝒩(0,(κ_δ,ϵ _2D)^2I_q),
is (ϵ,δ)-differentially private for the adjacency relation (<ref>). 
The estimate ẑ is then also (ϵ,δ)-differentially private, since it is obtained
by post-processing s, without re-accessing the sensitive signal y.





  §.§.§  Input Transformation Optimization


We can now consider the problem of optimizing the choice of matrix D. 
Let x̂_t^- = 𝔼[x_t | s_0:t-1] and x̂_t = 𝔼[x_t | s_0:t] be the 
state estimates produced by the Kalman filter of Figure <ref> after the 
prediction step and the measurement update step respectively <cit.>.
Let Σ̅_t = 𝔼[(x_t - x̂_t^-) (x_t - x̂_t^-)^T | s_0:t-1] and 
Σ_t  = 𝔼[(x_t - x̂_t) (x_t - x̂_t)^T | s_0:t] be the corresponding
error covariance matrices. 
We also denote by Σ̅_0 = diag(Σ^-_1,0, …Σ^-_n,0) 
the covariance matrix for the initial state x_0. 
For completeness, we recall here the equations of the Kalman filter.  
Given the dynamics (<ref>) and the measurement 
equation (<ref>), with ζ a Gaussian noise with covariance matrix 
(κ_δ,ϵ _2D)^2 I_q, we have for t ≥ 0 and starting
from x̂_0^- := x̅_0

    x̂_t    = x̂_t^- + K^f_t (s_t - DC_t x̂_t^-),  
    x̂_t+1^-    = A_t x̂_t + B_t u_t,     with 
    
    K^f_t    = Σ̅_t C_t^T D^T (D (C_t Σ̅_t C_t^T + V) D^T 
    + κ_δ,ϵ^2 Δ_2 D^2 I_q)^-1.

The error covariance matrices evolve for t ≥ 0 as

    Σ_t^-1 = Σ̅_t^-1+C_t^T Π C_t,   Σ̅_t+1   = A_tΣ_tA_t^T+W_t,

where 
Π = D^T(D V D^T+(κ_δ,ϵ _2D)^2I_q)^-1D.
With z = L_t x_t and its estimator ẑ_t = L_t x̂_t, we can rewrite
the MSE E_T in (<ref>) as

    E_T = 1/T+1∑_t=0^TTr (L_tΣ_t L_t^T).

As a result, a matrix D minimizing the MSE can be found by
solving the following optimization problem 


    min_q ∈ℕ, D ∈ℝ^q × p     1/T+1∑_t=0^TTr( L_tΣ_t L_t^T ) 	
    s.t. Σ_0^-1   =  Σ̅_0^-1 + C_0^T Π C_0, 	
    Σ_t+1^-1   = (A_tΣ_t A_t^T+ W_t)^-1 + C_t+1^T Π C_t+1,  
        0 ≤ t ≤ T-1, 
    Π   = D^T (D V D^T + κ_δ,ϵ^2 (Δ_2D)^2 I_q)^-1 D.


In the minimization (<ref>), we have emphasized that finding the first 
dimension q of the matrix D is part of the optimization problem.
Note also that we can write the optimization problem above equivalently as a minimization 
over the variables D, Π, and {Σ_t}_0 ≤ t ≤ T, but the variables 
other than D can be immediately eliminated using the equality 
constraints (<ref>)-(<ref>).

With our assumption V ≻ 0, we obtain an equivalent form for (<ref>)
by using the matrix inversion lemma 

    Π   =  V^-1 - V^-1( V^-1
    + D^T D/(κ_δ,ϵΔ_2 D)^2)^-1 V^-1,

or, alternatively,

    κ_δ,ϵ^2 [ ( V - V Π V )^-1 - V^-1] = 
    D^T D/Δ_2 D^2  .





  §.§.§ Semidefinite Programming-based Synthesis
     

In this section, we show that the optimization problem (<ref>)-(<ref>)
can be recast as a semidefinite program (SDP) and hence solved efficiently <cit.>,
if we impose the following additional constraints on D

    Δ_2 D = 1 = ρ_1 D_1_2 = … = ρ_n D_n_2.

First, the following Lemma shows that in fact no loss of performance occurs 
by adding the constraint (<ref>) to 
(<ref>)-(<ref>), i.e., that this constraint 
is satisfied automatically by some matrix D^* that is optimal 
for (<ref>)-(<ref>). 

	
For any feasible solution 
D of (<ref>)-(<ref>) that does not satisfy 
(<ref>), there exists a feasible solution that does satisfy this constraint 
and gives a lower or equal cost in (<ref>). 
In particular, adding the constraint (<ref>) to the problem 
(<ref>)-(<ref>) does not change the value of the
minimum nor the existence of a minimizer.



Consider a matrix D and a corresponding sequence {Σ_t}_0 ≤ t ≤ T defined by the 
iterations (<ref>)-(<ref>).
First, rescaling D to λ D for any λ≠ 0 does not impact the constraint 
(<ref>) (note that Δ_2 (λ D) = λΔ_2 D), and so we can add
the constraint Δ_2 D = 1 without changing the solution of the optimization problem 
(<ref>)-(<ref>).

Next, if the other constraints of (<ref>) are not satisfied by D, construct 
the p × p matrix M = D^T D + diag({η_i I_p_i}_1 ≤ i ≤ n),
with η_i = (Δ_2 D / ρ_i)^2 - D_i^2_2. Since η_i ≥ 0 by (<ref>),
M is positive semi-definite. The i^th diagonal block of M is
M_ii = D_i^T D_i + [(Δ_2 D / ρ_i)^2 - D_i^2_2] I_p_i, which has
maximum eigenvalue (Δ_2 D / ρ_i)^2. Define some matrix D̃ such that
D̃^T D̃ = M and group the columns of D̃ as 
D̃ = [  D̃_1 …D̃_n ] as for D, so that
D̃_i consists of p_i columns. In particular M_ii = D̃_i^T D̃_i, so D̃_i^T D̃_i
has maximum eigenvalue (Δ_2 D / ρ_i)^2 and hence D̃_i has maximum singular
value (Δ_2 D / ρ_i). In other words, D̃ satisfies (<ref>)
with a sensibility Δ_2 D = Δ_2 D̃ that is unchanged, and moreover D̃^T D̃ = M ≽ D^T D.

Therefore, when we replace D by D̃, Δ_2 D in the denominator of (<ref>) remains unchanged,
and moreover 

    ( V^-1 + D̃^TD̃/(κ_δ,ϵΔ_2 D̃)^2)^-1   ≼( V^-1 + D^T D/(κ_δ,ϵΔ_2 D)^2)^-1

hence Π̃≽Π,
where Π̃ and Π are defined according to (<ref>) or equivalently (<ref>) 
for D̃ and D respectively. Let K := Π̃- Π≽ 0.

Replacing Π by Π̃ in (<ref>), we obtain a matrix Σ̃_0
satisfying Σ̃_0^-1 = Σ_0^-1 + C_0^T K C_0≽Σ_0^-1,
so Σ̃_0≼Σ_0.
Now if we have two matrices Σ̃_t≼Σ_t, and we use these two matrices 
together with Π̃ and Π to define Σ̃_t+1, Σ_t+1 according to (<ref>), 
then immediately 

    Σ̃_t+1^-1   = (A_ tΣ̃_ t A_ t^T 
    + W_ t)^-1 + C_t+1^T Π̃C_t+1
       ≽ (A_ tΣ_ t A_ t^T + W_ t)^-1 + 
    C_t+1^T Π̃C_t+1
       = Σ_t+1^-1 + C_t+1^T K C_t+1≽Σ^-1_t+1.

Therefore, Σ̃_t+1≼Σ_t+1.

Hence, by induction, starting from D̃ we obtain a sequence {Σ̃_t}_t such 
that Σ̃_t≼Σ_t for all t ≥ 0.
This gives a smaller or equal cost

    1/T+1∑_t=0^T+1Tr (L_tΣ̃_t L_t^T) ≤1/T+1∑_t=0^T+1Tr (L_tΣ_t L_t^T),

and so the lemma is proved.


By Lemma <ref>, we can add without loss of optimality the constraints 
(<ref>) to (<ref>)-(<ref>), which allows us
in the following to recast the problem as an SDP. 
Let α_i = κ_δ,ϵρ_i, for all 1 ≤ i ≤ n.  
Denote E_i=[     0     … I_p_i     …     0 ]^T 
the p × p_i matrix whose elements are zero except for an identity matrix in 
its i^th block. The next lemma converts the constraints 
(<ref>)-(<ref>) to linear 
matrix inequalities. 

 	
If Π, D satisfy the constraints (<ref>)-(<ref>),
then Π satisfies V - V Π V ≽ 0 together with the following constraints, 
for all 1 ≤ i ≤ n,

    [ I_p_i/α_i^2 + V_i^-1                E_i^T;                  E_i              V-V Π V ]≽0, 
       [ I_p_i/α_i^2 + V_i^-1                E_i^T;                  E_i              V-V Π V ]≻0.
 
Conversely, if Π satisfies these constraints, then there exists a matrix D 
such that Π, D satisfy (<ref>)-(<ref>). 
One such D can be obtained by the factorization of

    κ_δ,ϵ^2 [ ( V - V Π V )^-1 - V^-1] = D^T D

(e.g., via singular value decomposition (SVD)) and will then satisfy Δ_2 D = 1.




V - V Π V ≽ 0 is immediate from (<ref>), since it is equal to 
( V^-1 + D^T D/(κ_δ,ϵΔ_2 D)^2)^-1.
Together with (<ref>), the right-hand side 
of (<ref>) then represents any positive 
semidefinite matrix M = D^T D such that
its diagonal blocks M_ii = D^T_i D_i have maximum eigenvalue
equal to 1/ρ_i^2, since D_i_2 = 1/ρ_i by (<ref>).
These constraints are equivalent to saying that for all 1 ≤ i ≤ n,

    E_i^T[ ( V - V Π V )^-1 - V^-1] E_i≼ I_p_i/α_i^2, 
       and not E_i^T[ ( V - V Π V )^-1 - V^-1] E_i 
    ≺ I_p_i/α_i^2.
 
Indeed, this comes from the standard fact that the maximum value λ_i, max 
of M_ii is the smallest λ satisfying M_ii≼λ I_p_i. 
The constraints given in the Lemma are obtained by 
noting that E_i^T V^-1E_i = V_i^-1
and taking Schur complements in (<ref>) and (<ref>).

Note that the fact that the left-hand side of (<ref>) is 
positive semidefinite is a simple consequence of V ≽ V - V Π V, hence adding
the constraint ( V - V Π V )^-1 - V^-1≽ 0
is unnecessary.


Next, define the information matrices Ω_t = Σ_t^-1, for 0 ≤ t ≤ T. 
If the matrices W_t are invertible, denoting
Ξ_t = W_t^-1 and using the matrix inversion lemma in (<ref>), one gets

    C_t+1^TΠ C_t+1 - Ω_t+1 + Ξ_t
          - Ξ_t A_t (Ω_t + A_t^TΞ_t A_t)^-1 A_t^TΞ_t = 0.


Replacing the equality in (<ref>) by ≽ 0 and taking a Schur complement,
together with  the inequalities of Lemma <ref>,
leads to the following SDP with variables 
Π≽ 0, {X_t≽ 0,Ω_t≻ 0}_0 ≤ t ≤ T



    min_Π≽ 0, {X_t, Ω_t}_0 ≤ t ≤ T  1/T+1∑_t=0^TTr (X_t)   s.t.
       [   X_t   L_t; L_t^T   Ω_t ]≽ 0,    0 ≤ t≤ T, 
       Ω_0 = Σ̅_0^-1 + C_0^TΠ C_0, 
       [ C_t+1^TΠ C_t+1 - Ω_t+1 + Ξ_t                      Ξ_t A_t;                     A_t^TΞ_t           Ω_t + A_t^TΞ_t A_t ]≽ 0, 
        0 ≤ t ≤ T-1, 
       [ I_p_i/α_i^2 + V_i^-1                E_i^T;                  E_i              V-V Π V ]≽0,   1 ≤ i≤ n.


Here the minimization of the cost (<ref>) has been replaced 
by the minimization of (<ref>), after introducing the 
slack variable X_t satisfying (<ref>), or equivalently 
X_t ≽ L_t Ω_t^-1 L_t^T by taking a Schur complement.
Since we replaced the equality in (<ref>) by an inequality, the SDP above
is a relaxation of the original problem (<ref>)-(<ref>).
The purpose of the next theorem is to show that this relaxation is tight.
Once an optimal solution for this SDP is obtained, we recover an optimal matrix D
from Π by the factorization (<ref>).


 
Let Π^* ≽ 0,{ X^*_t≽ 0,Ω^*_t≻ 0}_0 ≤ t ≤ T 
be an optimal solution for (<ref>)-(<ref>).
Suppose that for some 0 ≤ t ≤ T, we have L_t(Ω^*_t)^-1C_t^T ≠ 0.
Let D^* be a matrix obtained from Π^* by the factorization (<ref>).
Then D^* is an optimal solution for (<ref>)-(<ref>), which moreover
satisfies D_i^*_2 = 1/ρ_i for 1 ≤ i ≤ n, with the decomposition (<ref>).
The corresponding optimal covariance matrices {Σ_t^*}_0 ≤ t ≤ T for the
Kalman filter can be computed using the equations (<ref>)-(<ref>).
Finally, the optimal costs of (<ref>)-(<ref>) and 
(<ref>)-(<ref>) are equal, i.e., the
SDP relaxation is tight.

 

Even though the condition L_t(Ω^*_t)^-1C_t^T ≠ 0 introduced to guarantee 
the possibility of constructing the matrix D in the proof is not an explicit condition expressed 
directly in term of the problem parameters, it appears to be a weak requirement in practice. 


 

Consider Π^*, {X_t^*,Ω_t^*}_0 ≤ t ≤ T an optimal solution 
of the SDP (<ref>)-(<ref>). 
As explained in the proof of Lemma (<ref>), the constraint (<ref>)
is equivalent to

    α_i^2 E_i^T[ ( V - V Π^* V )^-1 - V^-1] E_i≼ I_p_i.

We show that we cannot have 

    α_i^2 E_i^T[ ( V - V Π^* V )^-1 - V^-1] E_i ≺ I_p_i.
 
Indeed, otherwise there exists η > 0
such that the matrix Π̃= Π^* + η I_p still satisfies (<ref>).
Using this matrix Π̃ in (<ref>), we obtain a matrix
Ω̃_0 = Ω_0^* + η C_0^T C_0 feasible for (<ref>).
Now define Ω̃_1 = Ω_1^* + η C_1^T C_1. One can immediately check that 
Π̃, Ω̃_0 and Ω̃_1 satisfy (<ref>)
for t = 0, using the fact that Ω_0^*, Ω_1^*, Π^* are feasible and that
C_0^T C_0 ≽ 0. Similarly the matrices Ω̃_t = Ω_t^* + η C_t^T C_t
are feasible in (<ref>) for all 0 ≤ t ≤ T.
Now in (<ref>), taking a Schur complement, we obtain that the
matrices X̃_t = L_t Ω̃_t^-1 L_t^T are feasible. By the matrix
inversion lemma we can write

    X̃_t = X^*_t - L_t (Ω^*_t)^-1 C_t^T K_t C_t (Ω^*_t)^-1 L_t^T.

for some matrices K_t ≻ 0.
These matrices X̃_t give a cost 
1/T+1∑_t=0^T Tr(X^*_t) - L_t (Ω^*_t)^-1 C_t^T K_t^1/2_F,
which is a strict improvement over the assumed optimal solution as soon as one matrix 
L_t (Ω^*_t)^-1 C_t^T is not zero (since the K_t's are invertible).
Hence, we have a contradiction and so we cannot have 
α_i^2 E_i^T[ ( V - V Π^* V )^-1 - V^-1] E_i ≺ I_p_i. 
We can then apply Lemma <ref> and construct a matrix D^* from Π^* 
as in (<ref>), so that the pair Π^*, D^* satisfies 
(<ref>)-(<ref>).

Let 𝒱^* be the optimum value of (<ref>)-(<ref>),
and V^* that of (<ref>)-(<ref>). First, 𝒱^* ≤ V^* since the
constraints of the original problem have been relaxed to obtain the SDP.
We now show how to construct a sequence {Σ_t^* }_0 ≤ t ≤ T, which together
with Π^* satisfy the constraints of (<ref>)-(<ref>)
and achieve the cost 𝒱^*, thereby proving the remaining claims of the theorem.
Note that since Ω^*_t + A_t^TΞ_t A_t≻ 0, (<ref>) 
is equivalent to ℛ_t(Ω^*_t, Ω^*_t+1) ≽ 0, where

    ℛ_t(Ω_t, Ω_t+1) :=     C_t+1^TΠ^* C_t+1 - Ω_t+1 + Ξ_t
       - Ξ_t A_t (Ω_t + A_t^TΞ_t A_t)^-1 A_t^TΞ_t.



First, we take Σ_0^* = (Ω_0^*)^-1. 
If ℛ_t(Ω^*_t, Ω^*_t+1) = 0 for all 0 ≤ t ≤ T-1,
then the matrices Ω^*_t satisfy (<ref>)
and we can take Σ_t^* = (Ω_t^*)^-1 for all t, since these
matrices satisfy the equivalent condition (<ref>).
Otherwise, let t̃ be the first time index such that 
ℛ_t̃(Ω^*_t̃, Ω^*_t̃+1) is not zero.
For t ≤t̃, we take Σ_t^* = (Ω_t^*)^-1 and so in particular
we have ℛ_t̃((Σ_t̃^*)^-1, Ω^*_t̃+1) ≽ 0
and not zero.
Consider the matrix 
Ω̃_t̃+1 = Ω^*_t̃+1 + ℛ_t̃(Ω^*_t̃, Ω^*_t̃+1), 
which then satisfies ℛ_t̃(Ω^*_t̃, Ω̃_t̃+1) = 0 
by definition. We set Σ_t̃+1^* = Ω̃_t̃+1^-1.

Now note that we again have 
ℛ_t̃+1((Σ_t̃+1^*)^-1, Ω̃_t̃+2) ≽ 0,
by verifying that (<ref>) is satisfied at t+1, using the fact that
(Σ_t̃+1^*)^-1 = Ω̃_t̃+1≽Ω^*_t̃+1.
From here, we can proceed by induction, assuming that Σ_0^*, …, Σ_t^* 
are set and taking

    Σ_t+1^* = (Ω̃_t+1)^-1 := (Ω_t+1^*+ℛ_t((Σ_t^*)^-1, Ω^*_t+1))^-1,

which reduces to (Ω_t+1^*)^-1 if ℛ_t((Σ_t^*)^-1, Ω^*_t+1) = 0.


The procedure above provides matrices Π^*, {Σ_t^* }_0 ≤ t ≤ T 
satisfying the constraints of the original program (<ref>)-(<ref>).
By construction, we have (Σ_t^*)^-1≽Ω_t^* and the matrices (Σ_t^*)^-1 
also satisfy (<ref>). Therefore, replacing, for each 0 ≤ t ≤ T, Ω_t^* by 
(Σ^*_t)^-1 and X_t^* by L_t Σ_t L_t^T in the solution 
of (<ref>)-(<ref>) that we started with gives a cost
𝒱≤𝒱^* for the SDP, hence 𝒱 = 𝒱^* by optimality of 𝒱^*.
But this cost 𝒱 is also equal to the cost 1/T+1∑_t=0^TTr (L_tΣ_t^* L_t^T) 
of (<ref>). Hence, we have shown that
(<ref>)-(<ref>) and
(<ref>)-(<ref>) have
the same minimum value, and constructed an optimal solution D^*, Π^*, {Σ_t^* }_0 ≤ t ≤ T
to (<ref>)-(<ref>) achieving this value.

 



 §.§ Stationary problem
 

In the stationary case with T →∞ and the model 
(<ref>)-(<ref>) now assumed time-invariant 
and detectable, we wish to find a signal aggregation matrix D followed 
by a time-invariant Kalman filter to minimize the steady-state MSE E_∞.
This can be done by solving the following SDP with variables Π≽ 0, 
X ≽ 0, Ω≻ 0


    min_Π≽ 0, X, Ω  Tr (X)   s.t.
       [   X   L; L^T   Ω ]≽ 0, 
       [ C^TΠ C - Ω + Ξ            Ξ A;           A^TΞ     Ω + A^TΞ A ]≽ 0, 
       [ I_p_i/α_i^2+ V_i^-1               E_i^T;                 E_i             V-V Π V ]≽0,   1 ≤ i≤ n.



Compared to (<ref>)-(<ref>), this SDP is of much smaller
size, due to the fact that the transient behavior is neglected in the performance measure.
The proof of the following theorem is similar to that of Theorem <ref>.
 
Let Π^* ≽ 0, X^* ≽ 0, Ω^* ≻ 0 be an optimal solution 
for (<ref>)-(<ref>).
Suppose that we have L (Ω^*)^-1 C^T ≠ 0.
Let D^* be a matrix obtained from Π^* by the factorization (<ref>).
Then D^* minimizes the steady-state MSE E_∞ among all possible matrices D 
introduced as in Figure <ref>, and the corresponding value of E_∞
is equal to the optimal value of the SDP.



Note that given the optimum matrix D^* and corresponding Π^*, an alternative way of 
computing E_∞ is by solving an ARE to obtain the steady-state prediction error 
covariance matrix Σ̅_∞ for x̂^-, then compute the steady-state error 
covariance matrix Σ_∞ = (Σ̅_∞^-1 + C^T Π^* C)^-1 for 
the estimator x̂, and finally E_∞ = (L Σ_∞ L^T).








 §.§ Syndromic Surveillance Example
 

To illustrate the differentially private filtering methodology, including issues related
to the choice of model and adjacency relation (<ref>), we discuss in this
section an example motivated by the analysis of epidemiological data.

Consider a scenario in which Public Health Services (PHS) must publish for a population
infected by a disease the number I_t of infectious people, i.e., those who have the disease 
and are able to infect others. 
PHS use privacy-sensitive data collected from n=12 hospitals, 
with each hospital i recording the number I_i,t of infectious people in its area,
as well as the number R_i,t of recovered people, i.e, those who were infected
by the disease and are now immune.
For each area i, these numbers are assumed to follow 
a discrete-time SEIR epidemiological model for the specific disease <cit.>,
written here first without process noise,
obtained by discretizing a classical continuous-time model <cit.> using 
a forward Euler discretization <cit.>

    S_i,t+1    =S_i,t-β_i S_i,tI_i,t/N_i, 
    
    E_i,t+1    =(1-τ_i)E_i,t+β_i S_i,tI_i,t/N_i,
    
    I_i,t+1    =(1-ϑ_i)I_i,t+τ_i E_i,t, 
    
    R_i,t+1    =R_i,t+ϑ_i I_i,t,

where S_i represents the number of susceptible people, i.e., those who are not infected 
but could become infected, the number of exposed people, i.e., those infected but not yet 
able to infect others, and N_i the total number of people.
The parameters τ_i, β_i and ϑ_i represent the transition rates from 
one disease stage to the next. 

For each hospital i's area, N_i in (<ref>) is assumed constant  
for the time interval of interest.
Moreover, let us make an approximation that this period  
is short enough or the disease at an early-enough stage so that
S_i,t can also be assumed approximately constant, equal to S_i,0.
Then, the remaining states 
η_i,t=[E_i,t, I_i,t, R_i,t]^T ∈ℝ^3
evolve as a linear system of the form

    η_i,t+1   = 𝒜_i η_i,t+ φ_i,t,    0 ≤ t ≤ T-1,

where 

    𝒜_i=    [         1-τ_i β_i S_i,0/N_i             0;           τ_i         1-ϑ_i             0;             0           ϑ_i             1 ],

and we introduced the white Gaussian process noise φ_i,t∼𝒩(0,Φ_i)
in the model, with covariance matrices Φ_i≻ 0, for 1 ≤ i ≤ 12. 
Now, consider the problem of choosing the level ρ_i in the adjacency
relation (<ref>) to provide a meaningful privacy guarantee to the
patients. If we assume that the measurement for hospital i in our model is 
y_i,t = [I_i,t, R_i,t]^T, then once a person becomes sick, he or she is 
counted at each period either in the signal I_i or, after recovery, in the 
signal R_i. As a result, the impact of a single individual on the measurement y_i 
could be quite large, proportional to the time horizon T, requiring in turn a large 
value of ρ_i to provide strong privacy guarantees, and hence a high level
of noise.


A practical solution to this issue is to not continuously record the same individuals, 
but simply count at each time period t the number of newly infectious individuals,
i.e., y_i,t^(1) = I_i,t-I_i,t-1, as well as the number of newly recovered individuals,
y_i,t^(2) = R_i,t-R_i,t-1. Such a measurement model is much more beneficial from
a privacy point of view. Indeed, assuming that a given individual can only become
sick once, he or she can affect y_i,t^(1) by ± 1 for at most 2 periods t
(when the person becomes sick and recovers), and y_i,t^(2) for at most one period by 1.
As a result, one can take ρ_i = √(3) in (<ref>) to provide a strong
privacy guarantee, i.e., insensitivity of the published output to the complete record of 
a single individual. The measurement noise v_i,t in the model represents counting
errors, e.g., due to people not being diagnosed by the hospital.


To obtain a dynamic model compatible with the measurements y_i,t = [y_i,t^(1), y_i,t^(2)]^T,
define for each hospital i the 4-dimensional state 
x_i,t = [I_i,t-1, R_i,t - R_i,t-1, E_i,t, I_i,t]^T.
These states x_i,t evolve as (<ref>) with  

    A_i=    [             0             0             0             1;             0             0             0           ϑ_i;             0             0         1-τ_i β_i S_i,0/N_i;             0             0           τ_i         1-ϑ_i ],   
    B_i= [ 0; 0; 0; 0 ]

and the noise w_i,t = [ϖ_i,t, φ_1,i,t, φ_2,i,t, φ_3,i,t]^T
includes the components of φ_i introduced in (<ref>)
as well as a small independent Gaussian noise ϖ_i with small variance σ^2, 
added so that the covariance matrices of W_i = diag(σ^2,Φ_i) are invertible 
as required by our algorithms (ideally ϖ_i,t would be 0 since the first line
of A_i corresponds to the delay in the model).
The measurement matrices are immediately

    C_i=[ -1  0  0  1;  0  1  0  0 ],   1 ≤ i ≤ 12,

and one can then verify that the model is observable.
Let us assume the following values for the parameters in (<ref>)

    τ_i=0.2,  β_i S_i,0/N_i = 0.5,  ϑ_i=0.1,  for   1 ≤ i ≤ 3 
       τ_i=0.3,  β_i S_i,0/N_i = 0.3,  ϑ_i=0.5,  for   4 ≤ i ≤ 6 
       τ_i=0.5,  β_i S_i,0/N_i = 0.7,  ϑ_i=0.15,   for   7 ≤ i ≤ 9 
       τ_i=0.7,  β_i S_i,0/N_i = 0.6,  ϑ_i=0.3,  for   10 ≤ i ≤ 12.

Moreover, assume for all 1 ≤ i ≤ 12

    V_i=0.4   I_2,     Φ_i = [   0.3 -0.15     0; -0.15   0.3 -0.15;     0 -0.15   0.3 ].

The goal of 
the surveillance system is to continuously release,
at each period t, an estimate of the quantity 

    z_t=    ∑_i=1^n( [ 0 0 0 1 ]× x_i,t).





Let us set the privacy parameters to δ=0.02 and ϵ=ln(3) for example
(and ρ_i to √(3) as discussed above).
We design the two-stage architecture of Figure <ref> by first 
solving the stationary optimization 
problem (<ref>)-(<ref>),
which provides an optimal matrix Π^*.
Recall that the matrix D can be then obtained from the factorization (<ref>).
The number of rows q of D is then equal to the rank of the matrix
M^* := κ_δ,ϵ^2 [ ( V - V Π^* V )^-1 - V^-1].
Due to the numerical procedures, this rank will typically be maximal 
(here equal to 24). However, we plot on Fig. <ref> the 
ratios σ_i/σ_max of the singular values of M^*, with σ_max
the maximum singular value. We see that if we select for example only the singular
values σ_i ≥ 10^-4σ_max in the SVD of M^* and set the smaller
ones manually to 0, we obtain a matrix of rank 14, hence a matrix D with 14 rows
instead of 24. We then verify (by solving an algebraic Riccati equation) that the 
performance of the steady-state Kalman filter is left virtually unchanged by this 
truncation, with a steady-state MSE of about 160, i.e., a root mean square error (RMSE) 
of 12.65 for the estimate of the number I_t of infectious people. 
In contrast, the input perturbation mechanism (i.e., taking D = I) gives a steady-state 
MSE of 777 or RMSE on I_t of 27.87. 
Reducing the number of rows of D is also beneficial for example in terms 
of processing complexity of the Kalman filter, which now has fewer inputs.




For δ=0.01, we compare on Fig. <ref> the steady-state RMSE 
of the two-stage mechanism and the input perturbation architecture for 
different values of the privacy parameter ϵ. One can see 
that by aggregating the input signals, we obtain a much better performance, 
especially in the high-privacy regime (when ϵ is small).
Other measures of performance could also be of interest for the final filtering
architecture, such as the convergence time of the estimates. 
For illustration purposes, Fig. <ref> shows sample paths 
of differentially private estimates both for the two-stage mechanism and 
for the input perturbation mechanism.













§ DIFFERENTIALLY PRIVATE LQG CONTROL

 




We now turn to the LQG control problem introduced at the end of 
Section <ref>. For 
concreteness, we assume here that the control input u_t at time t 
can depend on the measurements y_0:t up to time t. 
It is straightforward to adapt the discussion to the case where 
only y_0:t-1 are available to compute u_t. By the separation 
principle <cit.> for
the standard LQG control problem (i.e., with no privacy constraint), 
the optimal control law for the system (<ref>)-(<ref>) 
and quadratic cost (<ref>) is of the form 
u_t = K^c_t x̂_t, where: i) x̂_t = 𝔼[x_t | y_0:t] is 
the MMSE estimator, computed by the Kalman filter (<ref>) independently 
of the design of the optimal control law; and ii) K^c_t is the optimal 
gain for the deterministic linear quadratic regulator (LQR) problem, i.e., 
assuming that w = 0 in (<ref>) and C = I_n, v = 0 
in (<ref>). 
In particular, since the sequence of control gains K^c_t can be precomputed, 
the LQG control problem is similar to the filtering problem considered in the 
previous section, with the desired published output ẑ_t = L_t x̂_t 
simply replaced by u_t = K^c_t x̂_t. This motivates the architecture proposed 
on Figure <ref> for differentially private 
LQG control, which, compared to Figure <ref>, aggregates the 
measured signals y_i before adding the privacy-preserving noise. 
Essentially, the only difference with the Kalman filtering problem is that the 
performance is measured by (<ref>) instead of the 
MSE (<ref>), so that the cost function in the optimization 
problem for the matrix D needs to be changed.

The following theorem summarizes the discussion above and the classical results (see for example
<cit.>) that allow us to formulate in the following 
an efficiently solvable optimization problem for the choice of aggregation matrix D on 
Figure <ref>. 


Given a choice of matrix D for the differentially private LQG control architecture of 
Figure <ref>, the control law u_t(y_0:t), t ≥ 0, 
minimizing the cost function (<ref>) takes the form

    u_t = K^c_t x̂_t,

where x̂_t is computed by the Kalman filter (<ref>) and the gains 
K^c_t  are precomputed independently of the filtering problem as 

    K^c_t := - (R_t + B_t^T P_t+1 B_t)^-1 B_t^T P_t+1 A_t

with the matrices P_t ≽ 0 given by P_T = Q_T and the backward Riccati 
difference equation

    P_t =    A_t^T P_t+1 A_t + Q_t 
       - A_t P_t+1 B_t (R_t + B_t^T P_t+1 B_t)^-1 B_t^T P_t+1 A_t.

Moreover, the optimal objective (<ref>) corresponding to this 
control law can be written

    J_T = J^c_T + J^f_T(D)

where 

    J^c_T = 1/T+1( x̅_0^T P_0 x̅_0 + (P_0 Σ̅_0) 
    + ∑_t=1^T (P_t W_t) )

is a term independent of D and

    J_T^f(D)    = 1/T+1∑_t=0^T-1(N_t Σ_t),    with 
    
    N_t    = Q_t + A_t^T P_t+1 A_t - P_t,    0 ≤ t ≤ T-1,

and the matrices Σ_t for 0 ≤ t ≤ T-1 defined by (<ref>)-(<ref>).


Note that the dependence on D in (<ref>) is due to the fact 
that the error covariance matrices 
Σ_t depend on D via (<ref>)-(<ref>).
Moreover, from (<ref>) we see that 

N_t defined in (<ref>) is positive semidefinite. 
Hence, we can define for all 0 ≤ t ≤ T-1 matrices L_t such that N_t = L_t^T L_t, 
and L_T = 0, to rewrite the cost (<ref>) as 1/T+1∑_t=0^T(L_t Σ_t L_t^T).
Minimizing this cost over the matrices D with the relations 
(<ref>)-(<ref>) leads to an optimal 
aggregation matrix D for the architecture of Figure  <ref>. 
The reformulation of this optimization problem as an SDP then follows exactly from the 
same argument as in Section <ref>,
which led to Theorem <ref>. In other words, we have 
the following result.

	
Let L_t, 0 ≤ t ≤ T-1, be any matrices obtained from the factorization

    L_t^T L_t = N_t, 0 ≤ t ≤ T-1,

with N_t defined by (<ref>), and let L_T = 0.
Let Π^* ≽ 0,{ X^*_t≽ 0,Ω^*_t≻ 0}_0 ≤ t ≤ T 
be an optimal solution for (<ref>)-(<ref>)
with this choice of matrices L_t.
Suppose that for some 0 ≤ t ≤ T, we have L_t(Ω^*_t)^-1C_t^T ≠ 0.
Let D^* be a matrix obtained from Π^* by the factorization (<ref>).
Then D^* minimizes the LQG cost (<ref>) among all the aggregation
matrices D of Figure <ref>. This cost is equal to
J_T^c + 1/T+1∑_t=0^T (X_t^*), with J_T^c defined in (<ref>).





 §.§ Stationary Problem
 

As in Section <ref> for the filtering problem, we can consider 
the steady-state LQG problem by letting T →∞ and assuming the 
model (<ref>)-(<ref>) and the weight matrices Q and R 
in the cost (<ref>) to be time-invariant.
We assume the model to be detectable and stabilizable and the pair
(A,Q^1/2) to be detectable, in order to be able to implement a stabilizing
LQG controller.
We can take the optimal gains K^c and K^f of the controller and 
the Kalman filter respectively to be also independent of time. 
Following Theorem <ref> 
and Proposition <ref>, we then immediately have the 
following result for the design of the optimal D matrix.


	
Let P be the positive semidefinite solution of the following algebraic Riccati equation

    P =     A^T P A + Q - A^T P B(R+B^T P B)^-1 B^T PA.

Let L be any matrix obtained from the factorization

    L^T L = A^T P A + Q  - P.

Let Π^* ≽ 0, X^* ≽ 0, Ω^* ≻ 0 be an optimal solution 
for (<ref>)-(<ref>),
for this choice of matrix L. Suppose that we have L (Ω^*)^-1 C^T ≠ 0.
Let D^* be a matrix obtained from Π^* by the factorization (<ref>).
Then D^* minimizes the steady-state LQG cost J_∞ 
among all possible matrices D introduced as in Figure <ref>,
the corresponding value of the cost is J_∞ = (PW) + (X^*).







 §.§ Numerical Simulations
 
 

We illustrate the above results numerically for n = 10 independent scalar systems, 
with states x_i,t evolving as first order systems with time-invariant
dynamics (<ref>), where

    A_1    = 1.1, A_2 = 0.85, A_3 = 0.84, A_4 = 0.7, A_5 = 0.75, 
    
    A_6    = 0.9, A_7 = 0.8, A_8 = 1.05, A_9 = 0.99, A_10 = 1,

C_i = 1, W_i = 0.02 and V_i = 0.1 for all 1 ≤ i ≤ 10,
and B is a 10 × 3 matrix with B_ij = 0 except for

    B_3,1   = B_6,1 = B_9,1 = 1 
    
    B_1,2   = B_4,2 = B_7,2 = B_10,2 = 1 
    
    B_2,3   = B_5,3 = B_8,3 = 1.

In other words, the published control signal is
3-dimensional, with control input u_1 simultaneously actuating
systems 3, 6, 9, u_2 actuating systems 1, 4, 7, 10 and
u_3 actuating systems 2, 5 and 8.
We wish to regulate the sum of the states ∑_i=1^10 x_i to 0, 
hence we take 
Q to be the 10 × 10 all-ones matrix 
and R = I_3 in (<ref>). 
We set the privacy parameters to δ=0.05 and ϵ=ln(3),
and ρ_i=1 for 1 ≤ i ≤ 10.


To design the differentially private LQG controller with signal aggregation
for the stationary problem, we compute the matrix L of 
Proposition <ref> and solve the optimization 
problem (<ref>)-(<ref>).
Following the methodology discussed at the end of Section <ref>,
we find that one can take the matrix D to be a 4 × 10 matrix
at the matrix factorization stage (<ref>).
The corresponding steady-state cost J_∞ is found to be 1.37,
whereas it is 2.17 for the input perturbation mechanism (i.e., with
D = I_10). Hence, signal aggregation results in a significant
improvement.
Figure <ref> shows a comparison of the
cost J_∞ for this problem, with the two architectures, 
for different values of ϵ. 
Finally, Figure <ref> illustrates the sample
paths obtained under closed-loop control with the differentially
private controllers. 
We see in particular on Figure <ref>
that the two-stage architecture provides a much better transient
behavior for the regulated average trajectory (or sum of trajectories)
compared to the input perturbation architecture, in addition to a better
steady-state performance.
















§ CONCLUSION
 

This paper considers the Kalman filtering and LQG optimal control problems
under a differential privacy constraint. We propose  an architecture combining 
an input stage aggregating the individual signals appropriately, the Gaussian 
mechanism to enforce differential privacy and a Kalman filter to reconstruct
the desired estimate. Optimizing the parameters of this architecture can be recast
as an SDP. Examples illustrate the performance improvements 
compared to the input perturbation mechanism, which adds noise directly on the individual
signals. The methodology is then adapted to propose a similar two-stage architecture
for an LQG control problem, where the goal is to compute a shared control broadcasted
to the agent population.
Future research could consider the extension of these ideas to nonlinear systems, 
improving on the input and output perturbation mechanisms
of <cit.>.
In addition, since the size of the SDP increases rapidly with the number of agents
(and the time horizon in the non-stationary case), it would be useful to develop 
numerical methods and a problem-specific solver that take advantage of the sparsity 
of the matrices involved in the constraints, as in <cit.> for example.




































































IEEEtran




































































































