# Encrypted Distributed Lasso for Sparse Data Predictive Control

Andreea B. Alexandru  Anastasios Tsiamis  George J. Pappas The authors are
with the Department of Electrical and Systems Engineering, University of
Pennsylvania, Philadelphia, PA 19104.
<EMAIL_ADDRESS>
###### Abstract

The least squares problem with $\ell_{1}$-regularized regressors, called
_Lasso_ , is a widely used approach in optimization problems where sparsity of
the regressors is desired. This formulation is fundamental for many
applications in signal processing, machine learning and control. As a
motivating problem, we investigate a sparse data predictive control problem,
run at a cloud service to control a system with unknown model, using
$\ell_{1}$-regularization to limit the behavior complexity. The input-output
data collected for the system is privacy-sensitive, hence, we design a
privacy-preserving solution using homomorphically encrypted data. The main
challenges are the non-smoothness of the $\ell_{1}$-norm, which is difficult
to evaluate on encrypted data, as well as the iterative nature of the Lasso
problem. We use a distributed ADMM formulation that enables us to exchange
substantial local computation for little communication between multiple
servers. We first give an encrypted multi-party protocol for solving the
distributed Lasso problem, by approximating the non-smooth part with a
Chebyshev polynomial, evaluating it on encrypted data, and using a more cost
effective distributed bootstrapping operation. For the example of data
predictive control, we prefer a non-homogeneous splitting of the data for
better convergence. We give an encrypted multi-party protocol for this non-
homogeneous splitting of the Lasso problem to a non-homogeneous set of
servers: one powerful server and a few less powerful devices, added for
security reasons. Finally, we provide numerical results for our proposed
solutions.

##  I Introduction

Sparsity and compressed sensing have been widely used in signal processing,
machine learning and control applications, especially in the big-data regime
and noisy environments [1, 2]. In high-dimensional problems, it is likely that
only a subset of features affects the observations. Hence, pursuing sparse
representations reduces the model complexity, prevents overfitting and helps
with overall interpretation. Since finding the solution with the minimal
number of non-zero coefficients is NP-hard, $\ell_{1}$-regularization has been
proposed as a convex method that promotes sparsity. Perhaps one of the most
used algorithms for sparse recovery has been the celebrated Lasso algorithm
(least absolute shrinkage and selection operator), which accounts for both
sparsity and potentially noisy data, using $\ell_{1}$-regularization.

For instance, Lasso has been used in signal reconstruction for medical
imaging, wireless communication and tracking; seismology applications;
portfolio optimization; text analysis [3, 2]. Many of these applications are
large-scale or involve data coming from multiple data sources. With the recent
widespread availability and development of cloud services, it seems an
attractive and cost effective solution to outsource the computations to the
cloud, when the data owner or querier lacks the computational resources and/or
expertise to locally perform them. Given the privacy-sensitive nature of
medical data, financial data, location data, energy measurements etc., on
which such problems are computed, and how they can be used to profile users or
mount attacks on critical infrastructure, the computations should not be
performed in the clear at the cloud service.

###  I-A Contributions

In order to deal with the privacy issues, we draw on cryptographic approaches,
specifically, on homomorphic encryption, which enables polynomial computations
by the cloud over encrypted data of the client. However, encrypted Lasso
brings new challenges: evaluating non-smooth functions on encrypted data, as
well as continuing computations over multiple iterations and time steps, which
generally requires refreshing the ciphertexts.

A conventional observation is that distributing a large optimization problem
to multiple servers improves the execution time by parallelizing smaller
subproblems. Apart from this, we note that _distributing the computation
allows a streamlined execution of encrypted iterations_. In particular, using
multiple servers allows us to perform a refresh operation at a substantially
reduced cost compared to performing it only at one server. This cheaper
refresh operation enables us to continue the encrypted computations over
multiple iterations, as well as to use a high degree polynomial to approximate
the gradient of the $\ell_{1}$-norm. Specifically, we propose:

  * •

an efficient distributed encrypted solution to Lasso problems using ADMM,
offering computational privacy of all the data, including intermediate
results;

  * •

an optimized implementation of the above protocol using an efficient Chebyshev
series evaluation for polynomial approximations and reducing the number of
ciphertext levels and operations.

We apply our cloud-based sparsity framework to the problem of data-based
predictive control. Our goal is to control an unknown system using only the
privacy-sensitive input-output data that are potentially noisy. Data-driven
control is a blooming research area and methods based on the behavioral
framework have received significant renewed interest [4, 5, 6, 7, 8, 9, 10]
since their original proposal [11, 12]. The idea of such methods is that the
state representation can be replaced by a data-based representation which only
uses the trajectories of the system, bypassing the need for system
identification. In the case of noisy data, inspired by [5, 10], we reformulate
the data-based predictive control as a lasso problem.

For this use-case, we propose:

  * •

a distributed encrypted solution for $\ell_{1}$-regularized data predictive
control, using an optimized implementation.

For better convergence, we customize this solution to split the problem
heterogeneously between a powerful server and a few less powerful machines.

###  I-B Comparison to related work

The usage of ADMM for private distributed optimization is not novel, see e.g.,
[13, 14], given its convenient formulation and splitting of the objective
function and variables. (For other distributed gradient based methods, see
e.g., [15].) However, our usage of distributed ADMM substantially differs from
previous works in the following: i) we start with centralized rather than
already distributed data, so we split the centralized problem in a way that
fits our privacy and low-power requirements; ii) we assume heterogeneous
servers and we split the computations differently depending on who performs
them; iii) the data at each server is not in the clear, which complicates the
computations; iv) the servers do not learn any of the data, including
intermediate iterates and results; this requires more complex computations to
privately perform nonlinear operations; v) the $\ell_{1}$-regularization term
is non-smooth and has nonlinear gradient, leading to updates of the global
primal variable that are incompatible with the mentioned ADMM works.

Compared to [16], we consider the approximation of a non-smooth nonlinear
objective function and use multiple servers to streamline the complex
iterative computation, rather than a two-party computation, and a more
powerful homomorphic encryption scheme to partially replace the blinded
communication necessary at every iteration in [16].

The authors in [17] propose a distributed ADMM for a Lasso problem, using an
threshold additively homomorphic encryption and SPDZ [18] for computing the
nonlinear operations. In contrast to their work, the data is not distributed
in the clear to the computing servers, meaning we have less flexibility with
respect to the local computations. However, we can choose to split the data in
the most convenient way to have a distributed convergence speed similar to the
centralized convergence, which [17] does not discuss. Another difference is
that the tools they use require them to communicate for every nonlinear
computation, such as multiplications and comparison operations, which require
a number of communication rounds dependent on the number of bits in the
messages. In our case, the servers only send two messages per iteration and
the method we employ also allows us to batch vectors and perform operations in
parallel for all elements of a vector.

In [19, 20], the authors propose distributed/federated training and evaluation
with multi-party fully homomorphic encryption for linear, logistic and neural
network models, using stochastic gradient descent. While we inspired our
solution from the multi-party fully homomorphic encryption tool, their setup
is different than ours: the data is either distributed in cleartext locally at
the servers or other data providers perform the preprocessing; and the
computations are different, leading to different optimizations: e.g., [19]
uses a combination of distributed and centralized bootstrapping operations,
[20] has only one recurring variable (the model) to bootstrap.

While differential privacy for the solution of the optimization problem [21,
22] is out of the scope of this paper, it is an avenue for future research in
privately determining parameters of the optimization problem.

Encrypted control, surveyed in [23], has been recently gaining momentum, due
to the strong privacy guarantees it offers even when the controller is located
on an untrusted platform. In the case of controlling a system with known model
and linear controller parameters, the line of work [24, 25] has shown how to
perform the computations at subsequent time steps without the need of
bootstrapping or ciphertext reset. In contrast, we deal with both unknown
model matrices and nonlinearities, which prevent the application of their
methods.

Our previous work [26, 27] provides confidentiality for a different
formulation of a data predictive control problem. The Lasso formulation in our
current work does not have a closed-form solution as before, which complicates
the computations on encrypted data. We also use a different architecture at
the cloud, in order to completely remove the client involvement during the
computation of the control input.

An interesting remark is related to [28], where the authors propose to use
multiple controllers in parallel that perform asynchronous local bootstrapping
to ensure that at least one has a control input ready at each time step. In
our case, we prefer multiple servers to perform a distributed bootstrapping in
order to reduce the time it takes to refresh the ciphertexts.

##  II Problem formulation

For a covariate matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$, a vector of
outcomes $\mathbf{b}\in\mathbb{R}^{m}$, the variable
$\mathbf{x}\in\mathbb{R}^{n}$ and a penalty parameter $\lambda>0$, the Lasso
problem in its Lagrangian (unconstrained) form is given by:

$\min_{\mathbf{x}}~{}\frac{1}{2}||\mathbf{A}\mathbf{x}-\mathbf{b}||_{2}^{2}+\lambda||\mathbf{x}||_{1}.$
(1)

For dependent covariates, there is no closed-form solution to (1) and many
iterative optimization algorithms have been proposed in the literature [2, Ch.
5]. For example, Lasso problems can be solved using proximal methods or
augmented Lagrangian methods, such as the Alternating Direction Method of
Multipliers (ADMM) [29], [2, Ch. 5]. Splitting the objective function in the
ADMM way, we get:

$\displaystyle\begin{split}\min_{\mathbf{x},\mathbf{z}}~{}&\frac{1}{2}||\mathbf{A}\mathbf{x}-\mathbf{b}||_{2}^{2}+\lambda||\mathbf{z}||_{1}\\\
\text{s.t.}~{}&~{}\mathbf{x}-\mathbf{z}=\mathbf{0}.\end{split}$ (2)

Let
$S_{\alpha}(\mathbf{x})=(\mathbf{x}-\alpha\mathbf{1})_{+}-(-\mathbf{x}-\alpha\mathbf{1})_{+}$
denote the soft thresholding operator. The ADMM algorithm for (2) is:

$\displaystyle\mathbf{x}^{k+1}$
$\displaystyle=\arg\,\min\limits_{\mathbf{x}}\left(\frac{1}{2}||\mathbf{A}\mathbf{x}-\mathbf{b}||_{2}^{2}+\frac{\rho}{2}||\mathbf{x}-\mathbf{z}^{k}+\mathbf{w}^{k}||_{2}^{2}\right)$
$\displaystyle=\left(\mathbf{A}^{\intercal}\mathbf{A}+\rho\mathbf{I}\right)^{-1}\left(\mathbf{A}^{\intercal}\mathbf{b}+\rho(\mathbf{z}^{k}-\mathbf{w}^{k})\right)$
$\displaystyle\mathbf{z}^{k+1}$
$\displaystyle=\arg\,\min\limits_{\mathbf{z}}\left(\lambda||\mathbf{z}||_{1}+\frac{\rho}{2}||\mathbf{x}^{k+1}-\mathbf{z}+\mathbf{w}^{k}||_{2}^{2}\right)$
(3) $\displaystyle=S_{\lambda/\rho}(\mathbf{x}^{k+1}+\mathbf{w}^{k})$
$\displaystyle\mathbf{w}^{k+1}$
$\displaystyle=\mathbf{w}^{k}+\mathbf{x}^{k+1}-\mathbf{z}^{k+1}.$

While for general optimization problems, the ADMM might converge slowly, for
the Lasso problem it is known to have a fast convergence of a few (tens of)
iterations for a large range of parameter $\rho>0$ [29]. ADMM is a fast choice
in the cases where a very high precision of the optimal result is not
required, which is the case in noisy control problems, like the one we
investigate in Section V.

Goals and privacy requirements.  A client outsources problem (1) to a cloud
service that has to compute the optimal solution based on the data from the
client.

We consider the cloud service to be semi-honest, meaning it does not deviate
from the client’s specifications, but can process the data it receives to
extract private information for its own profit. The cloud service can be a
conglomerate of $K$ servers (see Figure 1), possibly belonging to different
organizations, offering the guarantee that not all $K$ servers collude.

Under this adversarial model, we require client data confidentiality, i.e., an
adversary corrupting at most $K-1$ of the servers should not be able to infer
anything about the client’s inputs and outputs, which consist of the values of
the matrix $\mathbf{A}$ and the vector $\mathbf{b}$, any intermediate values,
and solution $\mathbf{x}^{\ast}$. The penalty $\lambda$ and parameter $\rho$
can be chosen by the cloud service or chosen by the client, but are public.

Figure 1: Schematic diagram of the problem, with a client outsourcing its
encrypted data to a cloud service, that is authorized to compute on the data,
but not to decrypt it.

##  III Preliminaries

###  III-A Homomorphic Encryption (HE)

A HE scheme is called partially homomorphic if it supports the encrypted
evaluation of either a linear polynomial or a monomial, somewhat or leveled
homomorphic if it supports the encrypted evaluation of a polynomial with a
finite degree and fully homomorphic if it supports the encrypted evaluation of
arbitrary polynomials.

A HE scheme has a computational security parameter $\kappa$ if all known
attacks against the scheme, including breaking the encryption or
distinguishing between the encryptions of different messages, take
$2^{\kappa}$ bit operations. In practice, at least 128 bits of security are
preferred [30].

The newer generations of HE schemes are based on lattices and rely on the Ring
Learning with Errors [31, 32] hardness problem. Each operation on ciphertexts
introduces some noise, with additions being cheap in terms of noise amount,
but multiplications being expensive. If the noise amount exceeds a threshold,
decryption yields an incorrect result, so parameters of HE schemes need to be
designed accordingly. The level in leveled HE schemes concerns the supported
number of sequential multiplications. We say that a fresh ciphertext is on
level $L$ and a multiplication consumes a level. A ciphertext on level 0 does
not accept any more multiplications. Leveled HE schemes can be turned into
fully HE schemes if a bootstrapping operation is enabled, which can refresh
the ciphertext after the levels were consumed, such that further operations
are allowed while still guaranteeing correct decryption.

If done locally at a server with no access to the private key, bootstrapping
is a very expensive operation, consuming around 10 levels and introducing more
noise (from encrypted approximation of non-polynomial functions); see [33, 34]
for details. Apart from the computationally intensive bootstrapping procedure,
all the prior and posterior operations are impacted, since ciphertexts are
required to have an extra 10 levels, leading to very large scheme parameters
and ciphertexts size, making centralized bootstrapping undesirable.

Instead of performing bootstrapping locally, a server can ask the client to
refresh a ciphertext on level 0. However, this implies more computation,
communication and availability from the client, which is often prohibitive. A
preferable solution is to use two or more servers for the computation and the
refreshing. However, corrupting only two servers can be attainable by an
adversary. Increasing the number of servers decreases the probability of an
adversary corrupting them all. _Distributed bootstrapping trades substantial
computation power to one round of communication and does not introduce as much
noise as the centralized bootstrapping_ , as described below.

In multi-party HE schemes, the private key is additively secret shared between
a number of servers, meaning that no proper subset of the servers can decrypt.
An important assumption is that an adversary cannot corrupt all servers at
once, hence the private key is never recovered. In [35], a multi-party HE
scheme is described, where servers can carry out the homomorphic computations
locally and only need to interact for decryption and bootstrapping. In our
scenario, the decryption will be performed at the client so we are only
interested in distributed bootstrapping. Intuitively, distributed
bootstrapping requires each server to use its local secret share of the
private key to perform a partial decryption, mask this result and send it to
the other servers. Summing up all the partial decryptions results in a
refreshed ciphertext with the desired number of levels that can be correctly
decrypted to the original message. The masking needs to provide statistical
privacy of the message, so we require the mask to be $>80$ bits larger than
the size of messages and that no overflow occurs (the statistical security
parameter is generally smaller than the computational security parameter).
This means that distributed bootstrapping consumes around 3 levels; see more
details in [35, 19].

In this paper, we work with the multi-party version [19] of the leveled HE
scheme CKKS [36]. Each plaintext/ciphertext is a polynomial in the ring of
integers of a cyclotomic field. This enables the encoding of multiple scalars
in a plaintext/ciphertext and performing single instruction multiple data
($\mathrm{SIMD}$) operations, which can bring major computation and storage
improvements. Abstracting the details away, the $\mathrm{SIMD}$ operations
that can be supported are addition, element-wise multiplication by a plaintext
or ciphertext and data slot rotations in ciphertexts. We will use $+$ and
$\odot$ for $\mathrm{SIMD}$ addition and multiplication and
$\rho(\mathbf{x},i)$ to denote the row vector $\mathbf{x}$ rotated to the left
by $i$ positions ($i<0$ means rotation to the right).

We will denote by $\mathrm{E_{v0}}(\mathbf{x})$ the encryption of the vector
$\mathbf{x}$ followed by trailing zeros and by
$\mathrm{E_{v\ast}}(\mathbf{x})$ the encryption of the vector $\mathbf{x}$
followed by junk elements (elements whose value we do not care about).

###  III-B Polynomial approximation and Chebyshev series

As described above, homomorphic encryption can evaluate polynomials on
encrypted values. However, other operations such as trigonometric functions,
divisions or comparisons are not supported. Therefore, we prefer to evaluate a
polynomial approximation of the non-polynomial functions. We choose to work
with the Chebyshev polynomial series rather than the more common Taylor power
series due to better precision and smaller approximation error. Specifically,
the Chebyshev series polynomial interpolation is a near-minimax approximation
of a continuous function on the interval $[-1,1]$ [37].

However, polynomial approximation is not a panacea: for non-smooth functions,
it gives a reasonable error only on relatively small intervals and using high
degree polynomials. We choose to use this method, rather than other encrypted
computation tools that can exactly evaluate non-smooth functions at the cost
of more communication, knowing that we are dealing with noisy systems, where
the small approximation errors are absorbed by noise.

##  IV Encrypted distributed Lasso

The setting we consider is the following: the client encrypts its data
$\mathbf{A},\mathbf{b}$ and secret shares its private key to a number of
servers. The servers are responsible to compute and return the solution of
problem (2) to the the client.

###  IV-A Challenge: Evaluating non-polynomial functions

In the steps (3) of the ADMM algorithm for problem (2), the soft thresholding
function is non-polynomial, yet we need to evaluate it on encrypted data when
computing $\mathbf{z}^{k+1}$. We deal with this challenge by approximating the
soft thresholding function using a polynomial on a fixed interval via the
Chebyshev series (we hardcode the coefficients for this function). If the
interval is not $[-1,1]$, we first apply a linear transformation to bring the
inputs to this interval.

In the context of encrypted evaluation, a high polynomial degree increases the
number of levels necessary for the computations, as well as the number of
homomorphic multiplications between ciphertexts, which are expensive
operations (compared to plaintext-ciphertext multiplications or additions).
While we cannot consume fewer than $\left\lceil\log n\right\rceil$ levels to
evaluate a polynomial of degree $n$, we can reduce the $O(n)$ homomorphic
multiplications from the naive evaluation. Specifically, we implement the
Paterson-Stockmayer algorithm [38], which reduces the number of homomorphic
multiplications to $\left\lceil\sqrt{2n}+\log n\right\rceil+\mathcal{O}(1)$ by
recursively evaluating polynomials of smaller degree. We modify the Paterson-
Stockmayer algorithm that works with power series to work with Chebyshev
series. The benefit of this algorithm compared to the naive evaluation is
visible after degree 5 and grows with the degree. As an example, to evaluate a
non-monic polynomial of degree 25, we require 5 levels and 11 homomorphic
multiplications between ciphertexts.

###### Remark 1

The “stability” of the ADMM iterations allows the value
$\mathbf{x}^{k+1}+\mathbf{w}^{k}$ to stay within a fixed interval, given in
Lemma 1. In practice, we choose this interval from prior simulation.

###### Lemma 1

Define $\mathbf{M}:=\mathbf{A}^{\intercal}\mathbf{A}+\rho\mathbf{I}$,
$\mathbf{n}:=\mathbf{M}^{-1}\mathbf{A}^{\intercal}\mathbf{b}$ and
$c:=\sqrt{n}\lambda/\rho||2\rho\mathbf{M}^{-1}-\mathbf{I}||_{2}+||\mathbf{n}||_{2}$.
Since $\sigma:=||\rho\mathbf{M}^{-1}||_{2}$ is in $(0,1]$, we have the
following bounds for the quantity
$||\mathbf{x}^{k+1}+\mathbf{w}^{k}||_{\infty}$ in (3), for all
$k=1,\ldots,K_{\mathrm{iter}}$:

$\displaystyle||\mathbf{x}^{k+1}+\mathbf{w}^{k}||_{\infty}$
$\displaystyle\leq\sigma^{k}||\mathbf{n}||_{2}+\frac{1-\sigma^{k}}{1-\sigma}c,$
$\displaystyle\text{if }\sigma<1$ (4)
$\displaystyle||\mathbf{x}^{k+1}+\mathbf{w}^{k}||_{\infty}$
$\displaystyle\leq||\mathbf{n}||_{2}+kc,$ $\displaystyle\text{if }\sigma=1.$

###### Proof:

Let $\mathbf{y}^{k}:=\mathbf{x}^{k}+\mathbf{w}^{k-1}$ Manipulating (3), we
get:

$\displaystyle\mathbf{y}^{k+1}$
$\displaystyle=(\mathbf{I}-\rho\mathbf{M}^{-1})\mathbf{y}^{k}+(2\rho\mathbf{M}^{-1}-\mathbf{I})\mathbf{z}^{k}+\mathbf{n}$
$\displaystyle=\rho\mathbf{M}^{-1}\mathbf{y}^{k}+(2\rho\mathbf{M}^{-1}-\mathbf{I})(\mathbf{z}^{k}-\mathbf{y}^{k})+\mathbf{n}.$

The expression of the thresholding operator gives the following bound:
$-\lambda/\rho\leq\mathbf{z}_{i}^{k}-\mathbf{y}_{i}^{k}\leq\lambda/\rho$, for
$i=1,\ldots,n$. Then, using the triangle inequality and submultiplicative
property:

$\displaystyle||\mathbf{y}^{k+1}||_{2}\leq||\rho\mathbf{M}^{-1}\mathbf{y}^{k}||_{2}+||(2\rho\mathbf{M}^{-1}-\mathbf{I})(\mathbf{z}^{k}-\mathbf{y}^{k})||_{2}+||\mathbf{n}||_{2}$
$\displaystyle\leq||\rho\mathbf{M}^{-1}||_{2}||\mathbf{y}^{k}||_{2}+\sqrt{n}\lambda/\rho||2\rho\mathbf{M}^{-1}-\mathbf{I}||_{2}+||\mathbf{n}||_{2}.$

We select $\mathbf{z}^{0}=\mathbf{w}^{0}=\mathbf{0}$, compress the geometric
progression and use
$||\mathbf{y}^{k+1}||_{\infty}\leq||\mathbf{y}^{k+1}||_{2}$, getting the
bounds in (4). ∎

###  IV-B Challenge: Evaluating iterations

Depending on the precision we choose, the polynomial approximation can have a
high degree, implying the need of bootstrapping in order to continue
operations in the subsequent iterations. We resolve this challenge by making
use of multiple servers in order to realize a cheaper bootstrapping compared
to a centralized bootstrapping and a less burdensome solution than requesting
action from the client.

To this end, we turn to the distributed version of ADMM [29], such that we use
the servers both to ease the computation of the optimal solution and to ensure
privacy through encrypted computations. We split the matrix $\mathbf{A}$ and
vector $\mathbf{b}$ into $K$ parts, each to be held by a server, and rewrite
(2) as:

$\displaystyle\begin{split}\min_{\mathbf{x}_{1},\ldots,\mathbf{x}_{K},\mathbf{z}}~{}&\frac{1}{2}\sum_{i=1}^{K}||\mathbf{A}_{i}\mathbf{x}_{i}-\mathbf{b}_{i}||_{2}^{2}+\lambda||\mathbf{z}||_{1}\\\
\text{s.t.}~{}&~{}\mathbf{x}_{i}-\mathbf{z}=\mathbf{0},\quad
i=1,2,\ldots,K.\end{split}$ (5)

The ADMM algorithm for problem (5) is, for $i=1,\ldots,K$:

$\displaystyle\begin{split}\mathbf{x}_{i}^{k+1}&=\left(\mathbf{A}_{i}^{\intercal}\mathbf{A}_{i}+\rho\mathbf{I}\right)^{-1}\left(\mathbf{A}_{i}^{\intercal}\mathbf{b}_{i}+\rho(\mathbf{z}^{k}-\mathbf{w}_{i}^{k})\right)\\\
\mathbf{z}^{k+1}&=\frac{1}{K}S_{\lambda/\rho}\left(\sum_{i=1}^{K}\mathbf{x}_{i}^{k+1}+\sum_{i=1}^{K}\mathbf{w}_{i}^{k}\right)\\\
\mathbf{w}_{i}^{k+1}&=\mathbf{w}_{i}^{k}+\mathbf{x}_{i}^{k+1}-\mathbf{z}^{k+1}.\end{split}$
(6)

Each server is given ciphertexts corresponding to
$\mathbf{A}_{i},\mathbf{b}_{i}$. We assume that there is a preprocessing step
where servers can precompute convenient ciphertexts that will be used often in
the online iterations, such as
$1/\rho\mathbf{A}_{i}^{\intercal}\mathbf{b}_{i}$ and
$\rho\left(\mathbf{A}_{i}^{\intercal}\mathbf{A}_{i}+\rho\mathbf{I}\right)^{-1}$.
As in the unencrypted case, the servers can use the matrix inversion lemma to
compute an inversion of a smaller matrix, which saves in offline computation.
Online, each server locally computes the encryptions of $\mathbf{x}_{i}$ and
$\mathbf{w}_{i}$, then communicates to the other servers the local sum
$\mathbf{x}_{i}^{k+1}+\mathbf{w}_{i}^{k}$, such that all servers are then able
to compute $\mathbf{z}^{k+1}$. So far, the only communication necessary is the
same as in the unencrypted ADMM.

###  IV-C Challenge: Realizing the fewest bootstrapping operations

Distributed bootstrapping requires all parties to start by holding the same
ciphertext and all parties to obtain that refreshed ciphertext. Bootstrapping
the ciphertext encrypting $\mathbf{z}^{k+1}$ seems attractive, because it is
global and its evaluation involves the most sequential multiplications.
However, this is not enough: $\mathbf{w}_{i}^{k+1}$ loses levels through
$\mathbf{x}_{i}^{k+1}$, which is the result of a multiplication; so we would
need to bootstrap also before computing $\mathbf{z}^{k+1}$, not just after.

Instead, we do the following trick. Each server already has to compute and
send a ciphertext encrypting $\mathbf{x}_{i}^{k+1}+\mathbf{w}_{i}^{k}$ to the
other servers in order to compute the global iterate $\mathbf{z}^{k+1}$. This
means that each server can then construct a packed ciphertext $c^{k+1}$
encrypting
$\big{[}(\mathbf{x}_{1}^{k+1}+\mathbf{w}_{1}^{k})^{\intercal}(\mathbf{x}_{2}^{k+1}+\mathbf{w}_{2}^{k})^{\intercal}\,\ldots\,(\mathbf{x}_{K}^{k+1}+\mathbf{w}_{K}^{k})^{\intercal}\big{]}$
and distributedly bootstrap it. Afterwards, each server can extract the
refreshed ciphertext containing its local value
$\mathbf{x}_{i}^{k+1}+\mathbf{w}_{i}^{k}$, as well as a refreshed ciphertext
containing $\sum_{i=1}^{K}\mathbf{x}_{i}^{k+1}+\mathbf{w}_{i}^{k}$ by
repeatedly rotating and summing the refreshed ciphertext $c^{k+1}$ (this takes
$O(K)$ operations). From this value, each server can locally compute its
encrypted iterates $\mathbf{w}_{i}^{k+1}$ and $\mathbf{x}_{i}^{k+1}$, while
doing only one bootstrapping operation per iteration rather than two.

Apart from $\mathbf{x}_{i}^{k+1}+\mathbf{w}_{i}^{k}$, the servers send one
more message to complete the bootstrapping operation, so there are two rounds
of communication per iteration, one at the smallest admissible level (dictated
by bootstrapping) and the other at the full number of levels required,
computed below.

Assume that offline quantities are refreshed. Define $l_{B}$ to be the number
of levels for a statistically secure distributed bootstrapping and $l_{P}$ to
be the number of levels necessary for the evaluation of
$S_{\lambda/\rho}(\cdot)$ at a desired precision: this is the degree of the
approximation polynomial plus one, coming from the linear transformation to
the interval $[-1,1]$ (we merge the multiplication by $1/K$ in the Chebyshev
coefficients). Hence, the fresh ciphertexts need to have $L=l_{B}+l_{P}+1$
levels, if we bootstrap once per iteration. Because $l_{P}$ is usually higher
than 5, bootstrapping once every few iterations would lead to larger
parameters and ciphertexts.

###  IV-D Encrypted protocol

Protocol IV-D describes the steps for privately solving the distributed Lasso
problem.

We use an optimized diagonal method [39] for encrypted matrix-vector
multiplication. Consider a matrix $\mathbf{S}\in\mathbb{R}^{n\times n}$ and a
vector $\mathbf{p}\in\mathbb{R}^{n}$. Denote the diagonals of $\mathbf{S}$ by
$\mathbf{d}_{i}$, for $i=0,\ldots,n-1$. Let $n_{1}:=\lceil\sqrt{n/2}\rceil$
and $n_{2}:=n/n_{1}$. The corresponding result
$\mathbf{q}=\mathbf{S}\mathbf{p}$ is:

$\displaystyle\mathbf{q}=\sum_{i=0}^{n-1}\mathbf{d}_{i}\odot\rho(\mathbf{p},i)$
(7)
$\displaystyle=\sum_{j=0}^{n_{2}-1}\rho\left(\sum_{k=0}^{n_{1}-1}\rho(\mathbf{d}_{j\cdot
n_{1}+k},-j\cdot n_{1})\odot\rho(\mathbf{p},k);j\cdot n_{1}\right).$ (8)

With (8), we need $n_{1}+n_{2}=O(\sqrt{n})$ homomorphic rotations, given
$\rho(\mathbf{d}_{j\cdot n_{1}+k},-j\cdot n_{1})$, instead of $n$ if we use
(7). In both cases we require $n$ homomorphic multiplications.

For a rectangular matrix, we need extended diagonals but the method is the
same. The header of the function that achieves this is
$\mathrm{MultDiag}(\mathbf{S},\mathbf{p})$, and we pass the matrix
$\mathbf{S}$ as separate ciphertexts encoding diagonals rotated accordingly
and the vector $\mathbf{p}$ encoded in a ciphertext with trailing zeros. (In
line 9 in Protocol IV-D, a masking is performed in order to satisfy the latter
requirement.) We implement $\mathrm{MultDiag}$ such that it returns a
ciphertext that encodes the result $\mathbf{q}$ with trailing zeros. This
method can be parallelized. $\mathrm{MultDiag}$ is performed locally at the
servers (lines 3 and 11).

The header $\mathrm{ApproxSoftT}(\mathbf{p})$ represents the implementation of
the Chebyshev interpolation element-wise for $\mathbf{p}$, for a given set of
coefficients that specify the degree of the approximation and an interval for
which the coefficients are valid. Internally, the input is normalized to the
interval $[-1,1]$ (this also masks the junk elements), such that the output is
a ciphertext encoding the result with trailing zeros. $\mathrm{ApproxSoftT}$
is performed locally (line 8 in Protocol IV-D).

$\mathrm{DBoot}(\mathbf{p})$ is a distributed protocol, where all servers
start with the ciphertext of the vector $\mathbf{p}$ and all servers obtain a
ciphertext that contains the refreshed vector $\mathbf{p}$ having a
predetermined number of levels (line 6 in Protocol IV-D). It implies one round
of communication between all servers, as described in Section III-A.

There is also an offline protocol that computes the input as listed in
Protocol IV-D. We mention that the client distributes the rows of $\mathbf{A}$
and $\mathbf{b}$ to the servers and the shares of the private key. The servers
compute $\mathbf{m}_{i}$ by multiplication and $\mathbf{M}_{i}$ by the matrix
inversion lemma and a high degree polynomial approximation for the inversion
function, and use rotation and masking in order to obtain the diagonal
representation needed. The servers collaborate to bootstrap the ciphertexts,
in order to start Protocol 1 on the desired level.

Protocol 1: Distributed encrypted protocol for (5) with _equal servers_ and
_equal data split_

1:Public parameters: public key $\mathrm{pk}$, number of servers $K$, number
of maximum iterations $K_{\mathrm{iter}}$. $C$: secret key $\mathrm{sk}$.
$S_{1},\ldots,S_{K}$: encryption of
$\mathbf{M}_{i}=\rho(\mathbf{A}_{i}^{\intercal}\mathbf{A}_{i}+\rho\mathbf{I})^{-1}$,
encryption of
$\mathbf{m}_{i}=\frac{1}{\rho}(\mathbf{A}_{i}^{\intercal}\mathbf{b}_{i})$,
share of the secret key $\mathrm{sk}_{i}$, the Chebyshev coefficients for
evaluating $S_{\lambda/\rho}(\cdot)$ on a given interval.

2:$C$: $\mathbf{x}^{\ast}$.

3:$S_{i=1,\ldots,K}$: set initial values
$\mathrm{E_{v0}}(\mathbf{x}^{0}_{i})$, $\mathrm{E_{v0}}(\mathbf{w}^{0}_{i})$,
$\mathrm{E_{v0}}(\mathbf{z}^{0})$ (the value of $\mathbf{z}^{k}$ is previously
agreed upon);

4:for $k=0,\ldots,K_{\mathrm{iter}}-1$ do

5: $S_{i=1,\ldots,K}$:
$\mathrm{E_{v0}}(\mathbf{x}_{i}^{k})=\mathrm{MultDiag}(\mathbf{M}_{i},\mathbf{m}_{i}+\mathbf{z}^{k}-\mathbf{w}_{i}^{k})$;

6: $S_{i=1,\ldots,K}$: compute and send to other servers the rotation of the
sum
$\mathrm{E_{v0}}(\mathbf{v}_{i}):=\rho(\mathbf{x}^{k+1}_{i}+\mathbf{w}^{k}_{i},-(i-1)n)$;

7: $S_{i=1,\ldots,K}$: assemble
$\mathrm{E_{v\ast}}(\mathbf{v}):=\mathrm{E_{v\ast}}([\mathbf{v}_{1}\,\mathbf{v}_{2}\,\ldots\mathbf{v}_{K}])$
by summing own ciphertext and all received shifted ciphertexts;

8: $S_{i=1,\ldots,K}$: perform own part in the distributed bootstrapping and
get
$\mathrm{E_{v\ast}}(\mathbf{v}^{b}):=\mathrm{DBoot}(\mathrm{E_{v\ast}}(\mathbf{v}))$;

9: $S_{i=1,\ldots,K}$: extract its refreshed sum of local iterates
$\mathrm{E_{v*}}(\mathbf{v}^{b}_{i})=\rho(\mathrm{E_{v\ast}}(\mathbf{v}^{b}),(i-1)n)$;

10: $S_{i=1,\ldots,K}$: rotate and sum $\mathrm{E_{v\ast}}(\mathbf{v}^{b})$ to
get
$\mathrm{E_{v\ast}}(\sum_{i=1}^{K}\mathbf{x}^{k+1}_{i}+\mathbf{w}^{k}_{i})$,
then compute
$\mathrm{E_{v0}}(\mathbf{z}^{k})=\mathrm{ApproxSoftT}(\frac{1}{K}\sum_{i=1}^{K}\mathbf{x}^{k+1}_{i}+\mathbf{w}^{k}_{i},\lambda/\rho)$;

11: $S_{i=1,\ldots,K}$:
$\mathrm{E_{v0}}(\mathbf{w}^{k+1}_{i})=[\mathbf{1}_{S}^{\intercal}\,\mathbf{0}^{\intercal}]^{\intercal}\odot\mathrm{E_{v*}}(\mathbf{v}^{b}_{i})-\mathrm{E_{v0}}(\mathbf{z}^{k})$;

12:end for

13:$S_{1}$: compute
$\mathrm{E_{v0}}(\mathbf{x}_{1}^{K_{\mathrm{iter}}})=\mathrm{MultDiag}(\mathbf{M}_{1},\mathbf{m}_{1}+\mathbf{z}^{K_{\mathrm{iter}}-1}-\mathbf{w}_{1}^{K_{\mathrm{iter}}-1})$
and send it to the client $C$;

14:$C$: decrypt and extract $\mathbf{x}^{\ast}$.

###### Proposition 1

Protocol IV-D achieves client data confidentiality with respect to semi-honest
servers, assuming at least one of the servers is honest.

###### Proof:

We use two theoretical results in the proof. First, we need the underlying
homomorphic encryption scheme, CKKS, to be semantically secure, in order to
ensure that an adversary that does not have access to the private key of the
scheme cannot decrypt or even distinguish ciphertexts encrypting different
values. This has been proven in [36] and assumes that the Decisional Ring
Learning with Errors problem is computationally hard [31]. We select the
scheme parameters to ensure that this problem is hard in practice,
specifically that it achieves a security level of 128 bits, according to the
Learning with Errors estimator of [32]. Second, we need that the interactive
part of the multi-party CKKS scheme, in our case, the distributed
bootstrapping protocol, preserves the indistinguishability of the ciphertexts
and does not reveal the private key. This has been proven in [19] and assumes
that at least one servers is honest and that the masks used in
$\mathrm{DBoot}$ are statistically hiding. Our adversarial model indeed
considers that at most $K-1$ servers can be corrupted and we choose the masks
to be 80 bits larger than the messages, while ensuring enough levels such that
the result does not overflow.

In Protocol IV-D, the data of the client is sent encrypted to the servers.
Since the servers cannot all collude in order to reveal the private key, and
no information about the messages underlying is leaked by viewing or computing
on the respective ciphertexts, Protocol IV-D achieves client data
confidentiality. ∎

##  V Case study: data predictive control

Consider that the client wants to control a linear system with unknown model
parameters. The goal is to compute a reference-tracking LQR control, based
only on precollected input-output data. We can reformulate this control
problem inspired by the behavioral framework [11, 4, 5, 6].

We need the concept of a block-Hankel matrix. For the input signal
$\mathbf{u}=\left[\begin{matrix}\mathbf{u}_{0}^{\intercal}&\mathbf{u}_{1}^{\intercal}&\ldots&\mathbf{u}_{T-1}^{\intercal}\end{matrix}\right]^{\intercal}\in\mathbb{R}^{mT}$
and a positive integer $L$, this is given by the following:

$\mathbf{H}_{L}(\mathbf{u}):=\left[\begin{matrix}\mathbf{u}_{0}&\mathbf{u}_{1}&\ldots&\mathbf{u}_{T-L}\\\
\mathbf{u}_{1}&\mathbf{u}_{2}&\ldots&\mathbf{u}_{T-L+1}\\\
\vdots&&\ddots&\vdots\\\
\mathbf{u}_{L-1}&\mathbf{u}_{L}&\ldots&\mathbf{u}_{T-1}\end{matrix}\right].$

By definition, the signal $\mathbf{u}$ is persistently exciting of order $L$
if $\mathbf{H}_{L}(\mathbf{u})\in\mathbb{R}^{mL\times(T-L+1)}$ is full row
rank.

For input and output data $\mathbf{u}^{d}\in\mathbb{R}^{mT}$ and
$\mathbf{y}^{d}\in\mathbb{R}^{pT}$, we construct block-Hankel matrices for $M$
samples for the past data and $N$ samples for the future data, where
$S:=T-M-N+1$ and $\mathbf{U}_{p}\in\mathbb{R}^{mM\times
S},\mathbf{U}_{f}\in\mathbb{R}^{mN\times
S},\mathbf{Y}_{p}\in\mathbb{R}^{pM\times
S},\mathbf{Y}_{f}\in\mathbb{R}^{pN\times S}$:

$\mathbf{H}_{M+N}(\mathbf{u}^{d})=:\left[\begin{matrix}\mathbf{U}_{p}\\\
\mathbf{U}_{f}\end{matrix}\right],\quad\mathbf{H}_{M+N}(\mathbf{y}^{d})=:\left[\begin{matrix}\mathbf{Y}_{p}\\\
\mathbf{Y}_{f}\end{matrix}\right].\\\ $ (9)

Assume we have data richness, i.e., the precollected input is persistently
exciting [11]. This requires that the precollected input signal has length at
least $(m+1)(M+N+n)-1$.

Fix a time $t$ and let $\bar{\mathbf{u}}_{t}=\mathbf{u}_{t-M:t-1}$ be the
batch vector of the last $M$ inputs. The batch vector $\bar{\mathbf{y}}_{t}$
of the last $M$ outputs is defined similarly. If $M\geq n$, the standard LQR
problem can be reformulated as the data predictive control problem [5] in
(10), where the state representation is replaced with the precollected data.
According to the behavioral framework, an input-output trajectory of a linear
system is in the image of the block-Hankel matrices for the precollected data,
i.e., the constraint in (10), where $\mathbf{g}$ is a preimage of the
trajectory. $\mathbf{Q}$, $\mathbf{R}$ are LQR costs and $\mathbf{r}$ is the
desired reference. The first $m$ elements of $\mathbf{u}^{\ast,t}$ are input
into the system in a receding horizon fashion and $\mathbf{y}^{\ast,t}$ is the
predicted output.

$\displaystyle\begin{split}\min_{\mathbf{g},\mathbf{u},\mathbf{y}}~{}&~{}\frac{1}{2}\sum_{k=t}^{N+t-1}\left(||\mathbf{y}_{k}-\mathbf{r}_{k}||_{\mathbf{Q}}^{2}+||\mathbf{u}_{k}||_{\mathbf{R}}^{2}\right)\\\
s.t.~{}&~{}\left[\mathbf{U}_{p}^{\intercal}\,\,\,\mathbf{Y}_{p}^{\intercal}\,\,\,\mathbf{U}_{f}^{\intercal}\,\,\,\mathbf{Y}_{f}^{\intercal}\right]^{\intercal}\cdot\mathbf{g}=\Big{[}\bar{\mathbf{u}}_{t}^{\intercal}\,\,\,\bar{\mathbf{y}}_{t}^{\intercal}\,\,\,\mathbf{u}^{\intercal}\,\,\,\mathbf{y}^{\intercal}\Big{]}^{\intercal}.\end{split}$
(10)

In practice, there will be noise affecting the output measurement, as well as
precision errors induced by encryption, which might prevent an exact solution
to the equality constraints for the past data in (10). Hence, we prefer a
relaxation of the equality constraints via an $\ell_{2}$-least-squares
approach with penalty weights $\lambda_{y}$ and $\lambda_{u}$. Then, we
rewrite (10) as a minimization problem depending only on $\mathbf{g}$ by
enforcing $\mathbf{u}=\mathbf{U}_{f}\mathbf{g}$ and
$\mathbf{y}=\mathbf{Y}_{f}\mathbf{g}$. We also batch the objective function,
using the same notation $\mathbf{Q},\mathbf{R},\mathbf{r}$ for the batched
costs and reference.

Finally, to avoid overfitting due to noisy data, we penalize the magnitude of
$\mathbf{g}$ through an $\ell_{1}$-regularization with penalty parameter
$\lambda_{g}$ in (11). The intuition behind this choice comes from the fact
that in the noiseless data predictive control formulation, the block-Hankel
matrix of the trajectory data has an inherent low-rank structure. Choosing an
$\ell_{1}$-regularization acts like a convex relaxation of imposing a low-rank
constraint–see [10, Thm. 4.6] for more details.

$\displaystyle\begin{split}\min_{\mathbf{g}}~{}&\frac{1}{2}\left(||\mathbf{Y}_{f}\mathbf{g}-\mathbf{r}_{t}||_{\mathbf{Q}}^{2}+||\mathbf{U}_{f}\mathbf{g}||_{\mathbf{R}}^{2}\right)+\\\
&~{}\lambda_{y}||\mathbf{Y}_{p}\mathbf{g}-\bar{\mathbf{y}}_{t}||_{2}^{2}+\lambda_{u}||\mathbf{U}_{p}\mathbf{g}-\bar{\mathbf{u}}_{t}||_{2}^{2}+\lambda_{g}||\mathbf{g}||_{1}.\end{split}$
(11)

Notice that (11) is indeed a Lasso problem:

$\displaystyle\begin{split}\min_{\mathbf{g}}~{}&\frac{1}{2}||\mathbf{H}\mathbf{g}-\mathbf{J}\mathbf{f}_{t}||_{2}^{2}+\lambda_{g}||\mathbf{g}||_{1},\end{split}$
(12)

where we have
$\mathbf{J}:=\mathrm{blkdiag}\left(2\lambda_{y}\mathbf{I},\,\mathbf{Q},\,2\lambda_{u}\mathbf{I},\,\mathbf{R}\right)^{1/2},\mathbf{f}_{t}:=\Big{[}\bar{\mathbf{y}}_{t}^{\intercal}\,\,\,\mathbf{r}_{t}^{\intercal}\,\,\,\bar{\mathbf{u}}_{t}^{\intercal}\,\,\,\mathbf{0}^{\intercal}\Big{]}^{\intercal},\mathbf{H}:=\mathbf{J}\left[\mathbf{Y}_{p}^{\intercal}\,\,\,\mathbf{Y}_{f}^{\intercal}\,\,\,\mathbf{U}_{p}^{\intercal}\,\,\,\mathbf{U}_{f}^{\intercal}\right]^{\intercal}$.
In the case a hybrid regularization
$\lambda_{g}\|\mathbf{g}\|_{1}+\mu_{g}\|\mathbf{g}\|_{2}^{2}$ is preferred, we
can use the same formulation (12) and appropriately modify $\mathbf{H}$ and
$\mathbf{f}_{t}$.

Our goal is to provide a solution that outsources to a cloud service the
computation of the optimal solution $\mathbf{g}^{\ast,t}$ of (12) and of
$\mathbf{u}^{\ast,t}=\mathbf{U}_{f}\mathbf{g}^{\ast,t}$, while ensuring client
data confidentiality for all time steps, as described in Section II.

##  VI Encrypted data predictive control

Following the discussion in Section IV, we write the problem (12) as a
distributed problem with split variables:

$\displaystyle\begin{split}\min_{\mathbf{g}_{1},\ldots,\mathbf{g}_{K},\mathbf{z}}~{}&\frac{1}{2}\sum_{i=1}^{K}||\mathbf{H}_{i}\mathbf{g}_{i}-(\mathbf{J}\mathbf{f})_{i}||_{2}^{2}+\lambda_{g}||\mathbf{z}||_{1}\\\
\text{s.t.}~{}&~{}\mathbf{g}_{i}-\mathbf{z}=\mathbf{0},\quad
i=1,2,\ldots,K.\end{split}$ (13)

In the context of this data predictive problem, when we perform a homogenous
split of the data (a split of equal size), the distributed solution converges
very slowly to the global optimal solution. The reason for that is that the
homogeneous sub-problems have a different optimal solution than the global
solution. To gain intuition, consider splitting the component matrices
$\mathbf{U}_{p},\mathbf{U}_{f},\mathbf{Y}_{p},\mathbf{Y}_{f}$ of $\mathbf{H}$
equally between the $K$ servers. This means that each server solves a local
optimization problem for the same system that generated the values, but being
given fewer samples than necessary to characterize the behaviour of the
system, i.e., losing persistency of excitation. The problem remains when
allocating a random set of rows of the equal size to the servers.

To avoid this issue, we prefer to unequally split the problem. Specifically,
we designate Server 1 to have most of the rows and the rest of the servers to
hold fewer. Because the local solution of Server 1 is close to the central
solution, the empirical convergence to the optimal solution is much faster.
However, the more servers we add, the slower the convergence (if we do not
weight contributions differently). A valid option is to have only one server
do all the computation, i.e., central ADMM, and request help only for the
distributed bootstrapping from the rest of the servers (recall that we require
multiple servers both for an efficient distributed bootstrapping and for
security of the private key). But since the rest of the servers would be idle
while the central server performs the computation, we prefer to distribute
some of the computations to them as well.

Let the matrix $\mathbf{H}_{1}$ denote the first $(m+p)M+pN$ rows of matrix
$\mathbf{H}\in\mathbb{R}^{(m+p)(N+M)\times S}$. We split the remaining rows of
$\mathbf{H}$ into blocks of $mN/(K-1)$ rows, denoted $\mathbf{H}_{i}$ for
$i=2,\ldots,K$. We similarly split
$\mathbf{H}^{\intercal}\mathbf{J}\in\mathbb{R}^{S\times(m+p)(N+M)}$ into
$\bar{\mathbf{H}}_{1}$ and $\bar{\mathbf{H}}_{2}\ldots,\bar{\mathbf{H}}_{K}$
and $\mathbf{f}_{t}\in\mathbb{R}^{m}$ into $\mathbf{f}_{1,t}$ and
$\mathbf{f}_{2,t},\ldots,\mathbf{f}_{K,t}$. We prefer to use more of less
powerful devices in order to increase the security threshold (by splitting the
secret key into more values) and reduce the cost of operating the cloud
service. To this end, we shift some of the computations from the less powerful
servers to the more powerful Server 1 and remove online communication between
the client and the less powerful servers. Protocol 2 differs from Protocol
IV-D in this different allocation of computation, described below.

First, we split problem (13) such that $\mathbf{f}_{i,t}=\mathbf{0}$, for
$i=2,\ldots,K$, see (14). Second, Servers $2,\ldots,K$ have an easier offline
computation, since they have to invert substantially smaller matrices than
Server 1, using the matrix inversion lemma. The bootstrapping step is done the
same as in Protocol IV-D, after all parties broadcast their local sums.
However, we let only the more powerful Server 1 perform the summation
$\sum_{i=1}^{K}\mathbf{g}_{i}^{k+1}+\mathbf{w}_{i}^{k}$ and the evaluation of
the soft thresholding approximation, and then send the result
$\mathbf{z}^{k+1}$ to the other less powerful servers (the ciphertext
$\mathbf{z}^{k+1}$ will have only $l_{B}+1$ levels so communication is cheap).
Then, all servers continue with the computation of $\mathbf{w}_{i}^{k+1}$ and
finish the iteration.

$\displaystyle\begin{split}\mathbf{g}_{1}^{k+1}&=\left(\mathbf{H}_{1}^{\intercal}\mathbf{H}_{1}+\rho\mathbf{I}\right)^{-1}\left(\bar{\mathbf{H}}_{1}^{\intercal}\mathbf{f}_{1}+\rho(\mathbf{z}^{k}-\mathbf{w}_{1}^{k})\right)\\\
\mathbf{g}_{i}^{k+1}&=\rho\left(\mathbf{H}_{i}^{\intercal}\mathbf{H}_{i}+\rho\mathbf{I}\right)^{-1}(\mathbf{z}^{k}-\mathbf{w}_{i}^{k}),~{}i=2,\ldots,K\\\
\mathbf{z}^{k+1}&=\frac{1}{K}S_{\lambda_{g}/\rho}\left(\sum_{i=1}^{K}\mathbf{g}_{i}^{k+1}+\sum_{i=1}^{K}\mathbf{w}_{i}^{k}\right)\\\
\mathbf{w}_{i}^{k+1}&=\mathbf{w}_{i}^{k}+\mathbf{g}_{i}^{k+1}-\mathbf{z}^{k+1},~{}i=1,\ldots,K.\end{split}$
(14)

Moreover, because of the way we split the time-varying vector
$\mathbf{f}_{t}$, such that the elements corresponding to Servers $2,\ldots,K$
are 0, there is no need for them to update with the latest values of
$\mathbf{u}_{t}$ and $\mathbf{y}_{t}$. This way, only Server 1 needs to have a
connection with the client. The ciphertexts communicated to the client are on
level 0 (the predicted input $\mathbf{u}^{\ast,t}$), while the ciphertexts
communicated from the client ($\mathbf{u}_{t}$ and $\mathbf{y}_{t}$ for
assembling $\mathbf{f}_{1,t}$) are on level $l_{B}+2$.

Nevertheless, if the servers have different capacity, the more powerful server
will likely have to wait on the other servers for the bootstrapping
synchronization (requiring completion of the computation for
$\mathbf{g}_{i}^{k+1}$). In the idle time, the more powerful Server 1 can
perform multiple local updates, which heuristically helps with convergence in
our problem.

###### Proposition 2

Protocol 2 achieves client data confidentiality with respect to semi-honest
servers, assuming at least one of the servers is honest.

The proof follows from the proof of Proposition 1, regardless of having the
servers perform different tasks, since all tasks involve computations only on
encrypted data.

Protocol 2: Distributed encrypted protocol for (13) with _unequal servers_ and
_unequal data split_ for one time step $t$

1:Public parameters: public key $\mathrm{pk}$, parameters of the system and
offline trajectory $m,p,N,M,S$, the number of servers $K$, number of maximum
iterations $K_{\mathrm{iter}}$. $C$:
$(\mathbf{u}_{\tau},\mathbf{y}_{\tau})_{\tau=0,\ldots,t}$. $S_{1}$: encryption
of
$\mathbf{M}_{1}=\rho(\mathbf{H}_{1}^{\intercal}\mathbf{H}_{1}+\rho\mathbf{I})^{-1}$,
encryption of
$\mathbf{F}_{1}=\frac{1}{\rho}(\mathbf{M}_{1}\mathbf{H}_{1}^{\intercal}\mathbf{J})$,
encryption of $\mathbf{U}_{f}$, encryption of
$\bar{\mathbf{y}}_{t},\mathbf{r}_{t},{\mathbf{u}}_{t}$, share of the secret
key $\mathrm{sk}_{1}$, the Chebyshev coefficients for evaluating the soft
threshold function for a given interval and bias $\lambda_{g}/\rho$.
$S_{2},\ldots,S_{K}$: encryption of
$\mathbf{M}_{i}=\rho(\mathbf{H}_{i}^{\intercal}\mathbf{H}_{i}+\rho\mathbf{I})^{-1}$,
share of the secret key $\mathrm{sk}_{i}$, for $i=2,\ldots,K$.

2:$C$: $\mathbf{u}_{t+1}$

3:$C$: send to $S_{1}$ the ciphertexts $\mathrm{E_{v0}}(\mathbf{u}_{t})$,
$\mathrm{E_{v0}}(\mathbf{y}_{t})$, $\mathrm{E_{v0}}(\mathbf{r}_{t})$;

4:$S_{1}$: assemble the ciphertext
$\mathrm{E_{v0}}(\mathbf{f}_{1,t})=\mathrm{E_{v0}}(\left[\bar{\mathbf{y}}_{t}^{\intercal}\,\mathbf{r}_{t}^{\intercal}\,\bar{\mathbf{u}}_{t}^{\intercal}\right]^{\intercal})$;

5:$S_{i=1,\ldots,K}$: set initial values
$\mathrm{E_{v0}}(\mathbf{g}^{0}_{i})$, $\mathrm{E_{v0}}(\mathbf{w}^{0}_{i})$,
$\mathrm{E_{v0}}(\mathbf{z}^{0})$ (the value of $\mathbf{z}^{k}$ is previously
agreed upon);

6:for $k=0,\ldots,K_{\mathrm{iter}}-1$ do

7: $S_{1}$: compute
$\mathrm{E_{v0}}(\mathbf{g}_{1}^{k})=\mathrm{MultDiag}(\mathbf{F}_{1},\mathbf{f}_{1,t})+\mathrm{MultDiag}(\mathbf{M}_{1},\mathbf{z}^{k}-\mathbf{w}_{1}^{k})$;

8: $S_{i=2,\ldots,K}$: compute
$\mathrm{E_{v\ast}}(\mathbf{g}_{i}^{k})=\mathrm{MultDiag}(\mathbf{M}_{i},\mathbf{z}^{k}-\mathbf{w}_{i}^{k})$;

9: $S_{i=1,\ldots,K}$: compute and send to the other servers the rotation of
the sum
$\mathrm{E_{v0}}(\mathbf{v}_{i}):=\rho(\mathbf{g}^{k+1}_{i}+\mathbf{w}^{k}_{i},-(i-1)S)$;

10: $S_{i=1,\ldots,K}$: assemble
$\mathrm{E_{v\ast}}(\mathbf{v}):=\mathrm{E_{v\ast}}([\mathbf{v}_{1}\,\mathbf{v}_{2}\,\ldots\mathbf{v}_{K}])$
by summing own ciphertext and all received shifted ciphertexts;

11: $S_{i=1,\ldots,K}$: perform own part in the distributed boostratpping and
get
$\mathrm{E_{v\ast}}(\mathbf{v}^{b}):=\mathrm{DBoot}(\mathrm{E_{v\ast}}(\mathbf{v}))$;

12: $S_{i=1,\ldots,K}$: extract its refreshed sum of local iterates
$\mathrm{E_{v*}}(\mathbf{v}_{i})=\rho(\mathrm{E_{v\ast}}(\mathbf{v}),(i-1)S)$;

13: $S_{1}$: rotate and sum the refreshed $\mathrm{E_{v\ast}}(\mathbf{v}^{b})$
to obtain
$\mathrm{E_{v\ast}}(\sum_{i=1}^{K}\mathbf{g}^{k+1}_{i}+\mathbf{w}^{k}_{i})$,
then compute
$\mathrm{E_{v0}}(\mathbf{z}^{k})=\mathrm{EvalApproxSoftT}(\frac{1}{K}\sum_{i=1}^{K}\mathbf{g}^{k+1}_{i}+\mathbf{w}^{k}_{i},\lambda_{g}/\rho)$;

14: $S_{1}$: send to all the other servers $\mathrm{E_{v0}}(\mathbf{z}^{k})$;

15: $S_{i=1,\ldots,K}$: compute
$\mathrm{E_{v0}}(\mathbf{w}^{k+1}_{i})=[\mathbf{1}_{S}^{\intercal}\,\mathbf{0}^{\intercal}]^{\intercal}\odot\mathrm{E_{v*}}(\mathbf{g}_{i})-\mathrm{E_{v0}}(\mathbf{z}^{k})$;

16:end for

17:$S_{1}$: compute
$\mathrm{E_{v0}}(\mathbf{u}^{\ast})=\mathrm{MultDiag}(\mathbf{U}_{f},\mathbf{g}^{K}_{1})$
and send the result to the client $C$; $\triangleright$ or directly obtain
only the first $m$ components by multiplying by
$[\mathbf{I}_{m}\,\,\mathbf{0}]\mathbf{U}_{f}$

18:$C$: decrypt $\mathbf{u}_{t+1}$, plug it in the system to measure
$\mathbf{y}_{t+1}$.

##  VII Numerical results

We consider a data-driven temperature control of a 4x4 stable system
representing a building with four rooms, with sampling time $300$ seconds, and
$M=4$, $N=8$, $T=84$. We add process noise and measurement noise, both zero
mean Gaussian with covariance $0.01\mathbf{I}$. We choose the cost matrices
and regularization terms
$\mathbf{Q}=300\mathbf{I},\mathbf{R}=\mathbf{I},\lambda_{g}=300,\lambda_{y}=\lambda_{u}=3000$.
The data was distributed among 3 servers: Server 1 holds 64 rows and Servers 2
and 3 hold 16 rows each. Convergence for the Lasso problem associated to one
time step of the above problem occurred after 20 ADMM iterations, for
$\rho=1200$. Figure 2 reflects the tracking performance of the data predictive
control problem with these parameters.

Figure 2: Comparison between the tracking performance of the data predictive
controller solving (12) exactly via the CVX solver, and solving (12) via
distributed ADMM with 3 servers and approximating the soft thresholding
function with a degree-11 polynomial. The vertical dashed line marks the first
$M$ time steps, corresponding to the initial offline data. The curves
represent the temperature measurements in the four rooms of the system.

We evaluate Protocol 2 on Ubuntu 18.04 on a commodity laptop with 8 GB of RAM
and Intel Core i7, implemented using the PALISADE library [40], using 8
threads. We set the parameters such that we get a security level of 128 bits,
i.e., we use a ciphertext modulus of 436 bits and a ring dimension of
$2^{14}$. We obtain 6 decimal places precision for the results. The average
time for the first iteration is 2.75 seconds and for any iteration afterwards
is 1.98 seconds. The time for Server 1 to assemble the vector
$\mathbf{f}_{1,t}$ from the client and to compute the prediction is 0.61
seconds. The client needs 0.07 seconds to decrypt the control input and to
encrypt the new measurement and the input. This gives a total computation time
of 38.44 seconds per solving the optimization problem, not taking
communication into account. The setup takes 4.5 seconds, and is performed once
for all subsequent iterations.

If we artificially add a 150 ms delay of communication
(serialization/deserialization and transport) and assume the machines send the
messages sequentially to the other machines, then the total computation time
increases by 6.45 seconds (450 ms delay per iteration, and 450 ms delay for
communication between Server 1 and Client).

Overall, the maximum amount of memory Server 1 needs to have is 1.22 GB, while
Server 2 and 3 need 0.52 GB.

To simulate less powerful devices, we run Servers 2 and 3 on 2 threads instead
of 8. The total time necessary for the 20 iterations increases to 49.96
seconds. The majority of the difference comes from the final operation of
bootstrapping (which can be done asynchronously, i.e., it occurs after the
communication so servers do not have to wait for each other): the total
bootstrapping time increases from 0.91 seconds to 1.48 seconds. The rest of
the difference comes from the fact that computing $\mathbf{g}_{2}^{k+1}$ and
$\mathbf{g}_{3}^{k+1}$ takes 1.06 seconds compared to the 0.83 seconds that
Server 1 needs to compute $\mathbf{g}_{1}^{k+1}$.

Figure 3: Timing for steps in one online iteration for one server for solving
Lasso problems of various dimensions via encrypted distributed ADMM with three
servers (Protocol IV-D). The legend shows the operation that takes the most
time in the step.

Because the servers only need to synchronize in order to bootstrap (Server 2
and 3 also wait for $\mathbf{z}^{k+1}$ from Server 1, but in our instance this
does not create idle time since Server 1 is more powerful), Server 1 can
either wait for the other servers to finish the computation or can perform two
local updates of $\mathbf{g}_{1}$ . Heuristically, since the local solution of
Server 1 is closer to the global optimal solution, performing more local
iterations helps convergence (in our example, by needing 18 iterations instead
of 20). Nevertheless, the computation of the control input is still ready in
one sixth of the sampling time.

In Figure 3, we show how the time for one ADMM iteration varies with the
dimension of the problem, i.e., number of columns of matrix $\mathbf{A}$ in
(2). The scheme parameters are the same as described above. The blue bar shows
the time for lines 3 and 4 in Protocol IV-D, effectively consisting of the
matrix-vector multiplication. The yellow bar shows the time for lines 5 and 6,
representing the preparation for bootstrapping and the bootstrapping itself.
Finally, the red bar represents lines 7–9 of Protocol IV-D, consisting mostly
of the polynomial evaluation. We want to stress that the bootstrapping and
polynomial evaluation are made independent from the dimension of the problem
through packing, which represents a great advantage when increasing the
dimension. On the other hand, for large dimensions, the encrypted matrix
multiplication takes most of the computational and memory effort, and other
methods that decrease storage and number of operations at the cost of more
levels might be preferable.

##  VIII Future work

In this work, as well as in most of encrypted protocols for control and
optimization, we assumed knowledge of some bounds on the variables, which can
themselves leak information. An important avenue of research is to design
encrypted and differentially private prior experiments in order to compute
differentially private bounds on parameters such as costs, penalties and
number of iterations. We can involve the client in these experiments, by
having it control a small device with one share of the secret key. In this
way, the client has to first agree on the computation that is being
effectuated to let the experiment continue and decrypt the result.

We will also perform several optimizations to reduce the memory consumption at
the server machines, to make it even more amenable to small devices.

## References

  * [1] E. J. Candes, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” _Communications on Pure and Applied Mathematics_ , vol. 59, no. 8, pp. 1207–1223, 2006.
  * [2] T. Hastie, R. Tibshirani, and M. Wainwright, _Statistical learning with sparsity: the lasso and generalizations_. CRC press, 2015.
  * [3] Z. Zhang, Y. Xu, J. Yang, X. Li, and D. Zhang, “A survey of sparse representation: algorithms and applications,” _IEEE access_ , vol. 3, pp. 490–530, 2015.
  * [4] C. De Persis and P. Tesi, “On persistency of excitation and formulas for data-driven control,” in _58th Conference on Decision and Control_. IEEE, 2019, pp. 873–878.
  * [5] J. Coulson, J. Lygeros, and F. Dörfler, “Data-enabled predictive control: In the shallows of the DeePC,” in _18th European Control Conference_. IEEE, 2019, pp. 307–312.
  * [6] J. Berberich, J. Köhler, M. A. Muller, and F. Allgower, “Data-driven model predictive control with stability and robustness guarantees,” _IEEE Trans. on Automatic Control_ , 2020.
  * [7] J. Coulson, J. Lygeros, and F. Dörfler, “Distributionally robust chance constrained data-enabled predictive control,” _arXiv preprint arXiv:2006.01702_ , 2020.
  * [8] H. J. Van Waarde, J. Eising, H. L. Trentelman, and M. K. Camlibel, “Data informativity: a new perspective on data-driven analysis and control,” _IEEE Trans. on Automatic Control_ , 2020.
  * [9] J. Berberich, J. Köhler, M. A. Müller, and F. Allgöwer, “Robust constraint satisfaction in data-driven mpc,” in _59th Conference on Decision and Control_. IEEE, 2020, pp. 1260–1267.
  * [10] F. Dörfler, J. Coulson, and I. Markovsky, “Bridging direct & indirect data-driven control formulations via regularizations and relaxations,” _arXiv preprint arXiv:2101.01273_ , 2021.
  * [11] J. C. Willems, P. Rapisarda, I. Markovsky, and B. L. De Moor, “A note on persistency of excitation,” _Systems & Control Letters_, vol. 54, no. 4, pp. 325–329, 2005.
  * [12] I. Markovsky and P. Rapisarda, “Data-driven simulation and control,” _Intl. Journal of Control_ , vol. 81, no. 12, pp. 1946–1959, 2008.
  * [13] C. Zhang, M. Ahmad, and Y. Wang, “ADMM based privacy-preserving decentralized optimization,” _IEEE Trans. on Information Forensics and Security_ , vol. 14, no. 3, pp. 565–580, 2018.
  * [14] Y. Ye, H. Chen, M. Xiao, M. Skoglund, and H. V. Poor, “Privacy-preserving incremental ADMM for decentralized consensus optimization,” _IEEE Trans. on Signal Processing_ , vol. 68, pp. 5842–5854, 2020.
  * [15] Y. Lu and M. Zhu, “Privacy preserving distributed optimization using homomorphic encryption,” _Automatica_ , vol. 96, pp. 314–325, 2018.
  * [16] A. B. Alexandru, K. Gatsis, Y. Shoukry, S. A. Seshia, P. Tabuada, and G. J. Pappas, “Cloud-based quadratic optimization with partially homomorphic encryption,” _IEEE Trans. on Automatic Control_ , 2020.
  * [17] W. Zheng, R. Popa, J. E. Gonzalez, and I. Stoica, “Helen: Maliciously secure coopetitive learning for linear models,” in _IEEE Symposium on Security and Privacy_ , vol. 1, 2019.
  * [18] I. Damgård, V. Pastro, N. Smart, and S. Zakarias, “Multiparty computation from somewhat homomorphic encryption,” in _Annual Cryptology Conference_. Springer, 2012, pp. 643–662.
  * [19] D. Froelicher, J. R. Troncoso-Pastoriza, A. Pyrgelis, S. Sav, J. S. Sousa, J.-P. Bossuat, and J.-P. Hubaux, “Scalable privacy-preserving distributed learning,” _Proceedings on Privacy Enhancing Technologies_ , vol. 2021, no. 2, pp. 323–347, 2021.
  * [20] S. Sav, A. Pyrgelis, J. R. Troncoso-Pastoriza, D. Froelicher, J.-P. Bossuat, J. S. Sousa, and J.-P. Hubaux, “POSEIDON: Privacy-preserving federated neural network learning,” _arXiv preprint arXiv:2009.00349_ , 2020.
  * [21] S. Han, U. Topcu, and G. J. Pappas, “Differentially private distributed constrained optimization,” _IEEE Trans. on Automatic Control_ , vol. 62, no. 1, pp. 50–64, 2017.
  * [22] X. Zhang, M. M. Khalili, and M. Liu, “Improving the privacy and accuracy of ADMM-based distributed algorithms,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 5796–5805.
  * [23] M. S. Darup, A. B. Alexandru, D. E. Quevedo, and G. J. Pappas, “Encrypted control for networked systems–an illustrative introduction and current challenges,” _arXiv preprint arXiv:2010.00268_ , 2020.
  * [24] J. H. Cheon, K. Han, H. Kim, J. Kim, and H. Shim, “Need for controllers having integer coefficients in homomorphically encrypted dynamic system,” in _57th Conference on Decision and Control_. IEEE, 2018, pp. 5020–5025.
  * [25] J. Kim, H. Shim, and K. Han, “Dynamic controller that operates over homomorphically encrypted data for infinite time horizon,” _arXiv preprint arXiv:1912.07362_ , 2019.
  * [26] A. B. Alexandru, A. Tsiamis, and G. J. Pappas, “Towards private data-driven control,” in _59th Conference on Decision and Control_. IEEE, 2020, pp. 5449–5456.
  * [27] ——, “Data-driven control on encrypted data,” _arXiv preprint arXiv:2008.12671_ , 2020.
  * [28] J. Kim, C. Lee, H. Shim, J. H. Cheon, A. Kim, M. Kim, and Y. Song, “Encrypting controller using fully homomorphic encryption for security of cyber-physical systems,” _IFAC-PapersOnLine_ , vol. 49, no. 22, pp. 175–180, 2016.
  * [29] S. Boyd, N. Parikh, and E. Chu, _Distributed optimization and statistical learning via the alternating direction method of multipliers_. Now Publishers Inc, 2011.
  * [30] M. R. Albrecht, M. Chase, H. Chen, J. Ding, S. Goldwasser, S. Gorbunov, S. Halevi _et al._ , “Homomorphic encryption standard,” _IACR Cryptol. ePrint Arch._ , vol. 2019, p. 939, 2019.
  * [31] V. Lyubashevsky, C. Peikert, and O. Regev, “On ideal lattices and learning with errors over rings,” in _Intl. Conference on the Theory and Applications of Cryptographic Techniques_. Springer, 2010, pp. 1–23.
  * [32] M. R. Albrecht, R. Player, and S. Scott, “On the concrete hardness of learning with errors,” _Journal of Mathematical Cryptology_ , vol. 9, pp. 169–203, 2015, https://lwe-estimator.readthedocs.io/en/latest.
  * [33] J. H. Cheon, K. Han, A. Kim, M. Kim, and Y. Song, “Bootstrapping for approximate homomorphic encryption,” in _Intl. Conference on the Theory and Applications of Cryptographic Techniques_. Springer, 2018, pp. 360–384.
  * [34] H. Chen, I. Chillotti, and Y. Song, “Improved bootstrapping for approximate homomorphic encryption,” in _Intl. Conference on the Theory and Applications of Cryptographic Techniques_. Springer, 2019, pp. 34–54.
  * [35] C. Mouchet, J. Troncoso-Pastoriza, J.-P. Bossuat, and J.-P. Hubaux, “Multiparty homomorphic encryption from ring-learning-with-errors,” _Cryptology ePrint Archive, Report 2020/304_ , 2020.
  * [36] J. H. Cheon, A. Kim, M. Kim, and Y. Song, “Homomorphic encryption for arithmetic of approximate numbers,” in _Intl. Conference on the Theory and Application of Cryptology and Information Security_. Springer, 2017, pp. 409–437.
  * [37] J. C. Mason and D. C. Handscomb, _Chebyshev polynomials_. CRC press, 2002.
  * [38] M. S. Paterson and L. J. Stockmeyer, “On the number of nonscalar multiplications necessary to evaluate polynomials,” _SIAM Journal on Computing_ , vol. 2, no. 1, pp. 60–66, 1973.
  * [39] S. Halevi and V. Shoup, “Algorithms in HElib,” in _Annual Cryptology Conference_. Springer, 2014, pp. 554–571.
  * [40] “PALISADE Lattice Cryptography Library (release 1.10.4),” https://palisade-crypto.org/, Sep. 2020.


# An application of discrete-time SEIR model to the COVID-19 spread

u.a.rozikov, s. k. shoyimardonov  U.Rozikova,b,c a V.I.Romanovskiy Institute
of Mathematics of Uzbek Academy of Sciences; b AKFA University, 1st Deadlock
10, Kukcha Darvoza, 100095, Tashkent, Uzbekistan; c Faculty of Mathematics,
National University of Uzbekistan<EMAIL_ADDRESS>Sobirjon Shoyimardonov.
V.I.Romanovskiy Institute of Mathematics of Uzbek Academy of Sciences,
Tashkent, Uzbekistan<EMAIL_ADDRESS>
###### Abstract.

The Susceptible-Exposed-Infectious-Recovered (SEIR) model is applied in
several countries to ascertain the spread of the coronavirus disease 2019
(COVID-19). We consider discrete-time SEIR epidemic model in a closed system
which does not account for births or deaths, total population size under
consideration is constant. This dynamical system generated by a non-linear
evolution operator depending on four parameters. Under some conditions on
parameters we reduce the evolution operator to a quadratic stochastic operator
(QSO) which maps 3-dimensional simplex to itself. We show that the QSO has
uncountable set of fixed points (all laying on the boundary of the simplex).
It is shown that all trajectories of the dynamical system (generated by the
QSO) of the SEIR model are convergent (i.e. the QSO is regular). Moreover, we
discuss the efficiency of the model for Uzbekistan.

###### Key words and phrases:

COVID-19, quadratic stochastic operator, fixed point, discrete-time, SEIR
model, epidemic

###### 2000 Mathematics Subject Classification:

34D20 (92D25).

##  1\. Introduction

Since December 2019, healthcare systems worldwide have been struggling with
management of the coronavirus disease 2019 (COVID-19) pandemic. About 58000
publications related to COVID-19111https://publons.com/publon/covid-19/..
Unfortunately, the disease does not stop spreading and nowadays, about 152 025
280 worldwide infections over 219 countries and territories (May 01,
2021)222https://www.worldometers.info/coronavirus/.

In the book [2] basic definitions in epidemiology and some mathematical models
such as SI, SIS, SEIS, SARS, SIRS are given. To study epidemics at the
population level, one of the most traditional mathematical model is the
Susceptible, Exposed, Infectious, Recovered (SEIR) compartmental model. In
[1], it was adapted a traditional SEIR epidemic model to the specific dynamic
compartments and epidemic parameters of COVID 19, as it spreads in an age-
heterogeneous community.

In classic SEIR models, four compartments are usually considered:

  * $S(t)$

\- the susceptible population at time $t$ (i.e., the class of individuals who
are healthy but can contract the disease);

  * $E(t)$

\- the exposed population (individuals who have contracted the virus but are
not yet symptomatic);

  * $I(t)$

\- the infected population (class of individuals exhibiting signs and
symptoms, have contracted the disease and are now sick with it), it is assumed
that infected population are also infectious;

  * $R(t)$

\- the recovered population (the number of individuals who have recovered and
cannot contract the disease again, no longer infect others).

In a closed system which does not account for births or deaths, the sum of
these compartments remains constant in time, i.e.

$S(t)+E(t)+I(t)+R(t)=N.$

Let us consider SEIR model [1]:

(1) $\begin{cases}\frac{dS}{dt}&=-\beta S(I+qE)/N\\\\[5.69054pt]
\frac{dE}{dt}&=\beta S(I+qE)/N-E/\delta\\\\[5.69054pt]
\frac{dI}{dt}&=E/\delta-I/\gamma\\\\[5.69054pt]
\frac{dR}{dt}&=I/\gamma\end{cases}$

where, the parameter $\beta$ is the average number of contacts per person per
time, multiplied by the probability of disease transmission via a contact
between a susceptible individual, and an individual carrying the virus. The
carrier can be either infected or exposed, with the fraction ${SI}/{N^{2}}$
representing the likelihood of an arbitrary contact to be between a
susceptible and an infectious individual, and the fraction ${SE}/{N^{2}}$
corresponding to the likelihood of a contact to be between a susceptible and
an exposed individual.

The model allows the possibility that a contact with an exposed individual may
have different probability of transmission than that made with an infected
individual, which is reflected in the scaling factor $q$). The transition rate
at which people are exposed then takes the form $-d(S/N)/dt=\beta
S(I+qE)/N^{2}$, leading to the first equation of (1). The rate of transfer
from the exposed to the infectious stage is a fraction $1/\delta$ the number
of exposed individuals, where $\delta$ is the average time for an exposed
individual to become infectious. The rate of recovery is a fraction $1/\gamma$
the infectious population, where $\gamma$ is the average time it takes a
person to die or recover once in the infectious stage.

This model has already been used in its original form for an early assessment
of the epidemic in Wuhan, China [1].

In the equations (1) we do the following replacements:

$x=\frac{S}{N},\ \ y=\frac{E}{N},\ \ u=\frac{I}{N},\ \ v=\frac{R}{N},\ \
a=\frac{1}{\delta},\ \ b=\frac{1}{\gamma}$

and obtain

(2) $\begin{cases}\frac{dx}{dt}&=-\beta x(u+qy)\\\\[5.69054pt]
\frac{dy}{dt}&=\beta x(u+qy)-ay\\\\[5.69054pt] \frac{du}{dt}&=ay-
bu\\\\[5.69054pt] \frac{dv}{dt}&=bu\end{cases}$

where $a>0,b>0,\beta\geq 0,q\geq 0$. Here we notice that $x+u+y+v=1$.

In this paper we study discrete time dynamical system, associated to the
system (2), which is generated by the operator $V$ defined as

(3) $V:\left\\{\begin{aligned} x^{(1)}&=x-\beta x(u+qy)\\\ y^{(1)}&=y-ay+\beta
x(u+qy)\\\ u^{(1)}&=u-bu+ay\\\ v^{(1)}&=v+bu\end{aligned}\right.$

##  2\. Reduction to Quadratic Stochastic Operators

_The quadratic stochastic operator_ (QSO) [4], [5], [11] is a mapping of the
simplex

(4) $S^{m-1}=\\{x=(x_{1},...,x_{m})\in\mathbb{R}^{m}:x_{i}\geq
0,\sum\limits_{i=1}^{m}x_{i}=1\\}$

into itself, of the form

(5)
$F:x^{\prime}_{k}=\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{m}P_{ij,k}x_{i}x_{j},\qquad
k=1,...,m,$

where the coefficients $P_{ij,k}$ satisfy the following conditions

(6) $P_{ij,k}\geq 0,\quad
P_{ij,k}=P_{ji,k},\quad\sum\limits_{k=1}^{m}P_{ij,k}=1,\qquad(i,j,k=1,...,m).$

Thus, each QSO $F$ can be uniquely defined by a cubic matrix
$\mathbb{P}=(P_{ij,k})_{i,j,k=1}^{m}$ with conditions (6).

Note that each element $x\in S^{m-1}$ is a probability distribution and each
such distribution can be interpreted as a state of the corresponding
biological system.

For a given $\lambda^{(0)}\in S^{m-1}$ the _trajectory_ (orbit)
$\\{\lambda^{(n)};n\geq 0\\}$ of $\lambda^{(0)}$ under the action of QSO (5)
is defined by

$\lambda^{(n+1)}=F(\lambda^{(n)}),\;n=0,1,2,...$

The _main problem_ in mathematical biology consists in the study of the
asymptotical behaviour of the trajectories.

######  Definition 1.

A QSO $F$ is called regular if for any initial point $\lambda^{(0)}\in
S^{m-1}$, the limit

$\lim_{n\to\infty}F^{n}(\lambda^{(0)})$

exists, where $F^{n}$ denotes $n$-fold composition of $F$ with itself (i.e.
$n$ time iterations of $F$).

######  Proposition 1.

For the operator (3) we have $V\left(S^{3}\right)\subset S^{3}$ if and only if

(7) $a,b,\beta\in[0,1],\ \ \mbox{and}\ \ \beta q\leq 1.$

Moreover, under this conditions the operator $V$, defined in (3), is a QSO.

###### Proof.

Using $x+y+u+v=1$ we rewrite the operator (3) (as in [7], [8], [9]):

(8) $V:\left\\{\begin{aligned} x^{(1)}&=x(x+y+u+v)-\beta x(u+qy)\\\
y^{(1)}&=y(x+y+u+v)(1-a)+\beta x(u+qy)\\\
u^{(1)}&=u(x+y+u+v)(1-b)+ay(x+y+u+v)\\\
v^{(1)}&=v(x+y+u+v)+bu(x+y+u+v)\end{aligned}\right.$

From this system and definition of QSO we get the following relations:

(9) $\begin{array}[]{cccc}\begin{aligned} {\scriptstyle
P_{11,1}}&={\scriptstyle 1},\\\ {\scriptstyle 2P_{14,1}}&={\scriptstyle 1},\\\
{\scriptstyle P_{22,2}}&={\scriptstyle 1-a},\\\ {\scriptstyle
2P_{12,3}}&={\scriptstyle a},\\\ {\scriptstyle 2P_{23,3}}&={\scriptstyle
1+a-b,}\ \ \\\ {\scriptstyle 2P_{34,3}}&={\scriptstyle 1-b},\\\ {\scriptstyle
2P_{23,4}}&={\scriptstyle b},\\\ {\scriptstyle 2P_{34,4}}&={\scriptstyle
1+b},\\\ \end{aligned}\par\begin{aligned} {\scriptstyle
2P_{12,1}}&={\scriptstyle 1-\beta q},\\\ {\scriptstyle
2P_{12,2}}&={\scriptstyle 1-a+\beta q},\\\ {\scriptstyle
2P_{23,2}}&={\scriptstyle 1-a},\\\ {\scriptstyle 2P_{13,3}}&={\scriptstyle
1-b},\\\ {\scriptstyle 2P_{24,3}}&={\scriptstyle a},\\\ {\scriptstyle
2P_{13,4}}&={\scriptstyle b},\\\ {\scriptstyle 2P_{24,4}}&={\scriptstyle
1},\\\ {\scriptstyle P_{44,4}}&={\scriptstyle 1},\\\
\end{aligned}\par&\begin{aligned} {\scriptstyle 2P_{13,1}}&={\scriptstyle
1-\beta},\\\ {\scriptstyle 2P_{13,2}}&={\scriptstyle\beta},\\\ {\scriptstyle
2P_{24,2}}&={\scriptstyle 1-a},\\\ {\scriptstyle P_{22,3}}&={\scriptstyle
a},\\\ {\scriptstyle P_{33,3}}&={\scriptstyle 1-b},\\\ {\scriptstyle
2P_{14,4}}&={\scriptstyle 1},\\\ {\scriptstyle P_{33,4}}&={\scriptstyle b},\\\
{\scriptstyle\text{other }}&{\scriptstyle
P_{ij,k}=0}\par\end{aligned}\end{array}$

Now it is easy to see that the conditions (6) on $P_{ij,k}$ are equivalent to
conditions (7). ∎

######  Remark 1.

Based on the work [1] we note that transmission rate $\beta$ changes between
$0.1-0.3,$ infection rate $\frac{1}{\delta}=a$ between $0.07-0.5$ (incubation
period 2-14 days), recovery rate $\frac{1}{\gamma}=b$ between $0.05-0.1$
(infectious period of 10-20 days). In addition, $\beta q$ is a probability of
transmission of the contacting with an exposed individual, so $\beta q\leq 1.$
Thus, the conditions (7) of Proposition 1 confirm the biological studies.
Therefore, below we consider the operator (3) with condition (7).

######  Proposition 2.

If parameters satisfy (7) then operator (3) is regular.

###### Proof.

For any initial point $\lambda^{(0)}\in S^{3}$ the trajectory
$\lambda^{(n)}=(x^{(n)},y^{(n)},u^{(n)},v^{(n)})=V^{n}(\lambda^{(0)})\in
S^{3}$, is given as

(10) $\left\\{\begin{aligned} x^{(n+1)}&=x^{(n)}-\beta
x^{(n)}(u^{(n)}+qy^{(n)})\\\ y^{(n+1)}&=y^{(n)}-ay^{(n)}+\beta
x^{(n)}(u^{(n)}+qy^{(n)})\\\ u^{(n+1)}&=u^{(n)}-bu^{(n)}+ay^{(n)}\\\
v^{(n+1)}&=v^{(n)}+bu^{(n)}\end{aligned}\right.,\ \ n\geq 0.$

Since $\lambda^{(n)}\in S^{3}$ we have that all its coordinates between $0$
and $1$. Therefore, by (10) we have $x^{(n+1)}\leq x^{(n)}$ and $v^{(n+1)}\geq
v^{(n)}$. Thus $x^{(n)}$ and $v^{(n)}$ are convergent as monotone and bounded
sequence. To show that $y^{(n)}$ and $u^{(n)}$ also have limits, consider
$a_{n}=x^{(n)}+y^{(n)}$, $b_{n}=u^{(n)}+v^{(n)}$ then by (10) we have

$a_{n+1}=a_{n}-ay^{(n)}\leq a_{n},\ \ b_{n+1}=b_{n}+ay^{(n)}\geq b_{n},$

i.e., $a_{n}$ is decreasing and $b_{n}$ is increasing. Both sequences are
positive and bounded (since $a_{n}+b_{n}=1$). Hence these sequences also have
limits. Consequently, from the following equalities it follows that $y^{(n)}$
and $u^{(n)}$ are convergent

$y^{(n)}=a_{n}-x^{(n)},\ \ u^{(n)}=b_{n}-v^{(n)}.$

This completes the proof (according to Definition 1). ∎

##  3\. Type of fixed points

Recall that a fixed point of the operator $V$ is a solution of $V(x)=x.$

One can see that the set of fixed points of the operator (3) is

(11) ${\rm Fix}(V)=\\{\Lambda(x)=(x,0,0,1-x):\forall x\in[0,1]\\}.$

Note that two vertices $e_{1}=(1,0,0,0),e_{2}=(0,0,0,1)$ of the simplex
$S^{3}$ belong in ${\rm Fix}(V).$

######  Remark 2.

Since $V$ is a continuous and regular operator the limit point of each
trajectory is a fixed point of the operator. Therefore, by Proposition 2 we
have that independently on the initial point the second and third coordinate
(i.e. $y^{(n)}$ and $u^{(n)}$ of the trajectory has limit

$\lim_{n\to\infty}y^{(n)}=\lim_{n\to\infty}u^{(n)}=0.$

######  Definition 2.

[3]. A fixed point $p$ for $F:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ is
called _hyperbolic_ if the Jacobian matrix $\textbf{J}=\textbf{J}_{F}$ of the
map $F$ at the point $p$ has no eigenvalues on the unit circle.

There are three types of hyperbolic fixed points:

  * (1)

$p$ is an attracting fixed point if all of the eigenvalues of $\textbf{J}(p)$
are less than one in absolute value.

  * (2)

$p$ is an repelling fixed point if all of the eigenvalues of $\textbf{J}(p)$
are greater than one in absolute value.

  * (3)

$p$ is a saddle point otherwise.

######  Proposition 3.

The fixed point $\Lambda(x^{*})$ is a nonhyperbolic for any $x^{*}\in[0,1].$

###### Proof.

In the system (3) we take first three equations:

(12) $W:\left\\{\begin{aligned} x^{(1)}&=x-\beta xu-\beta qxy\\\
y^{(1)}&=y(1-a)+\beta xu+\beta qxy\\\ u^{(1)}&=u(1-b)+ay\\\
\end{aligned}\right.$

where $x+y+u\leq 1.$
For the operator $W$ we calculate a Jacobian as follows:

$J=\left[\begin{array}[]{cccccc}1-\beta u-\beta qy&-\beta qx&-\beta x\\\ \beta
u+\beta qy&1-a+\beta qx&\beta x\\\ 0&a&1-b\\\ \end{array}\right]$

and at the fixed point $\Lambda(x^{*})$ it has the following form:

(13) $J(\Lambda(x^{*}))=\left[\begin{array}[]{cccccc}1&-\beta qx^{*}&-\beta
x^{*}\\\ 0&1-a+\beta qx^{*}&\beta x^{*}\\\ 0&a&1-b\\\ \end{array}\right]$

Obviously, for the matrix $J(\Lambda(x^{*}))$ at least one eigenvalue always
equal to 1. Similarly, using $x+y+u+v=1,$ for any triple of equations from the
system (3) it can be shown that Jacobian matrix always has an eigenvalue on
the unit circle. The proof is completed. ∎

Let $J(\Lambda(\alpha)),$ be the Jacobian matrix of the operator $W$, at fixed
point $\Lambda(\alpha)$, $\alpha\in[0,1]$. Then (13) has the following
characteristic equation:

(14) $(1-\mu)[(1-a+q\alpha\beta-\mu)(1-b-\mu)-a\alpha\beta]=0$

and the roots of this equation are:

(15) $\mu_{1}=1,\ \ \mu_{2}=1-\frac{a+b-q\alpha\beta-\sqrt{D}}{2},\ \
\mu_{3}=1-\frac{a+b-q\alpha\beta+\sqrt{D}}{2},$

where $D=(b-a+q\alpha\beta)^{2}+4a\alpha\beta\geq 0.$

By Proposition 3, the fixed point $\Lambda(\alpha)$ is a nonhyperbolic fixed
point.

######  Lemma 1.

Let $\Lambda(\alpha)$ be the fixed point of the operator $W$ and
$\mu_{2},\mu_{3}$ are eigenvalues defined in (15). Then

$|\mu_{2}|=\left\\{\begin{array}[]{lll}1\ \ \ \ \ \text{if}\ \
\alpha=\frac{ab}{\beta(a+bq)},\\\\[5.69054pt] <1\ \ \text{if}\ \
\alpha<\frac{ab}{\beta(a+bq)},\\\\[5.69054pt] >1\ \ \text{if}\ \
\alpha>\frac{ab}{\beta(a+bq)},\end{array}\right.$

and

$|\mu_{3}|<1\ \ \text{for any}\ \ \alpha\in[0,1].$

###### Proof.

Step-1. First we consider the condition $|\mu_{2}|<1.$

$|\mu_{2}|<1\Leftrightarrow 0<a+b-q\alpha\beta-\sqrt{D}<4$

since $a+b-q\alpha\beta-4<0$ one has $a+b-q\alpha\beta-4<\sqrt{D}.$

Let $a+b-q\alpha\beta>0,$ i.e., $\alpha<\frac{a+b}{q\beta}.$ Then

$\sqrt{D}<a+b-q\alpha\beta\Leftrightarrow
4a\alpha\beta+(b-a+q\alpha\beta)^{2}<(a+b-q\alpha\beta)^{2}\Leftrightarrow\alpha<\frac{ab}{\beta(a+bq)}.$

Thus, if $\alpha<\min\\{\frac{ab}{\beta(a+bq)},\frac{a+b}{q\beta}\\}$ then
$|\mu_{2}|<1.$ But always $\frac{ab}{\beta(a+bq)}<\frac{a+b}{q\beta}.$ Assume
the opposite:

$\frac{ab}{\beta(a+bq)}\geq\frac{a+b}{q\beta}\Leftrightarrow(a+b)(a+bq)\leq
abq\Leftrightarrow a^{2}+ab+b^{2}q\leq 0.$

Hence, $|\mu_{2}|<1$ if and only if $\alpha<\frac{ab}{\beta(a+bq)}.$ Moreover,
$|\mu_{2}|=1$ if $\alpha=\frac{ab}{\beta(a+bq)}$ and $|\mu_{2}|>1$ if
$\alpha>\frac{ab}{\beta(a+bq)}.$
Of course, there exists $\alpha$ satisfying the last inequality, because, by
Remark 1, $b\leq\beta$ and from this one has
$\frac{ab}{\beta(a+bq)}\leq\frac{ab}{\beta a}=\frac{b}{\beta}\leq 1.$
Step-2. Now, we consider the condition $|\mu_{3}|<1.$

$|\mu_{3}|<1\Leftrightarrow 0<a+b-q\alpha\beta+\sqrt{D}<4.$

If $q\alpha\beta-a-b<0$ i.e., $\alpha<\frac{a+b}{q\beta}$ then left inequality
holds true. If $\alpha\geq\frac{a+b}{q\beta}$ then from $q\alpha\beta-
a-b<\sqrt{D}$ one has $\alpha>\frac{ab}{\beta(a+bq)}.$ Thus, the inequality
$0<a+b-q\alpha\beta+\sqrt{D}$ is always true. Consequently, we have to check
$\sqrt{D}<4-(a+b-q\alpha\beta).$ Note that $4-(a+b-q\alpha\beta)>0.$

$\sqrt{D}<4-(a+b-q\alpha\beta)\Leftrightarrow
4a\alpha\beta+(b-a+q\alpha\beta)^{2}<16-8(a+b-q\alpha\beta)+(a+b-q\alpha\beta)^{2}$
$\Leftrightarrow
a\alpha\beta+2(a+b-q\alpha\beta)-b(a-q\alpha\beta)-4<0\Leftrightarrow$
$\Leftrightarrow
a<\frac{4+2q\alpha\beta-2b-bq\alpha\beta}{2-b+\alpha\beta}=\frac{(2-b)(2+q\alpha\beta)}{2-b+\alpha\beta}<2+q\alpha\beta.$

But last inequality and therefore $|\mu_{3}|<1$ always holds. This completes
the proof of the lemma. ∎

Let $\bar{x}$ be a fixed point. Then (see [12])

$E^{s}(\bar{x})=\emph{{\rm span}\\{eigenvectors associated with eigenvalues of
modulus$<1$\\}}$

is called a stable eigenspace

$E^{u}(\bar{x})=\emph{{\rm span}\\{eigenvectors associated with eigenvalues of
modulus$>1$\\}}$

is called an unstable eigenspace

$E^{c}(\bar{x})=\emph{{\rm span}\\{eigenvectors associated with eigenvalues of
modulus$=1$\\}}$

is called a center eigenspace.

Using Lemma 1 and [12] we define the followings about dimensions of
eigenspaces:

$\text{dim}E^{s}(\Lambda(\alpha))=1,\text{dim}E^{c}(\Lambda(\alpha))=2,\text{dim}E^{u}(\Lambda(\alpha))=0,\
\ \text{if}\ \ \alpha=\frac{ab}{\beta(a+bq)},$
$\text{dim}E^{s}(\Lambda(\alpha))=2,\text{dim}E^{c}(\Lambda(\alpha))=1,\text{dim}E^{u}(\Lambda(\alpha))=0,\
\ \text{if}\ \ \alpha<\frac{ab}{\beta(a+bq)},$
$\text{dim}E^{s}(\Lambda(\alpha))=1,\text{dim}E^{c}(\Lambda(\alpha))=1,\text{dim}E^{u}(\Lambda(\alpha))=1,\
\ \text{if}\ \ \alpha>\frac{ab}{\beta(a+bq)},$

Moreover, if $\bar{x}$ is a fixed point, then from the center manifold theory
([13]), there are the following invariant manifolds:
locally stable manifold $W_{loc}^{s}$ tangent to $E^{s}$ at $\bar{x};$
locally center manifold $W_{loc}^{c}$ tangent to $E^{c}$ at $\bar{x};$
locally unstable manifold $W_{loc}^{u}$ tangent to $E^{u}$ at $\bar{x}.$
locally center-stable manifold $W_{loc}^{cs}$ tangent to $E^{c}\oplus E^{s}$
at $\bar{x};$
locally center-unstable manifold $W_{loc}^{cu}$ tangent to $E^{c}\oplus E^{u}$
at $\bar{x};$

The stable and unstable manifolds are unique, but center, center-stable, and
center-unstable manifolds may not be unique [12].

##  4\. The limit points of trajectories

In this section we study the limit behavior of trajectories of initial point
$\lambda^{0}\in S^{3}$ under operator (3), i.e. the sequence
$V^{n}(\lambda^{0})$, $n\geq 1$.

######  Proposition 4.

If $\beta=0$ then for any initial point
$\lambda^{0}=\left(x^{0},u^{0},y^{0},v^{0}\right)\in S^{3}$ (except fixed
point) the trajectory has the following limit

$\lim_{n\to\infty}V^{(n)}(\lambda^{0})=(x^{0};0;0;1-x^{0})$

###### Proof.

If $\beta=0$ then $x^{(n)}=x^{0},\ \ y^{(n)}=(1-a)^{n}y^{0}\rightarrow 0.$
Since $v^{(n+1)}\geq v^{(n)},$ for any $n\in N,$ one has the sequence
$v^{(n)}$ has limit as a bounded increasing sequence. Moreover,
$u^{(n)}=1-x^{(n)}-y^{(n)}-v^{(n)}$ so the sequence $u^{(n)}$ also has a limit
$\bar{u}.$ We get a limit from $u^{(n+1)}=u^{(n)}(1-b)+ay^{(n)}$ and from
positiveness of $b$ one has $\bar{u}=0.$ Thus, the proof is completed. ∎

######  Remark 3.

For further results, we assume that $\beta>0$ and the conditions of the
Proposition 1 are satisfied.

######  Lemma 2.

The sets

$X=\\{(x,y,u)\in[0,1]^{3}:x=0,y+u\leq 1\\},$ $Y=\\{(x,y,u)\in[0,1]^{3}:x\leq
1,y=u=0\\},$ $M=\\{(x,y,u)\in[0,1]^{3}:x+y+u\leq 1,ay-\beta x(u+qy)>0,bu-
ay>0\\}$

are invariant sets with respect to operator (12).

###### Proof.

The invariantness of $X$ and $Y$ are straightforward. To show invariantness of
$M$, for shortness denote

$A=ay-\beta x(u+qy),\ \ B=bu-ay.$

Then first we rewrite $bu^{(1)}-ay^{(1)}$ as

$bu^{(1)}-ay^{(1)}=b(u-bu+ay)-a(y-ay+\beta x(u+qy))$ (16)
$=b(u-B)-a(y-A)=(1-b)B+aA.$

Next rewrite $ay^{(1)}-\beta x^{(1)}(u^{(1)}+qy^{(1)})$ as

$ay^{(1)}-\beta x^{(1)}(u^{(1)}+qy^{(1)})=a(y-ay+\beta x(u+qy))$
$-\beta(x-\beta x(u+qy))(u-bu+ay+q(y-ay+\beta x(u+qy)))$
$=a(y-A)-\beta(x-\beta x(u+qy))(u-B+q(y-A))$ $=ay-\beta x(u+qy)-aA+\beta
x(B+qA)+\beta^{2}x(u+qy)(u^{(1)}+qy^{(1)})$ $=(1-a)A+\beta
x(B+qA)+\beta^{2}x(u+qy)(u^{(1)}+qy^{(1)}).$

Since $u^{(1)}\geq 0,y^{(1)}\geq 0$ and by the conditions for the parameters
(by Proposition 1) in the last expression, each term is non-negative, so this
completes the proof of the lemma.

∎

######  Remark 4.

Note that the invariant set $Y$ belongs to the intersection of two surfaces
$ay-\beta x(u+qy)=0,bu-ay=0$ and the extension of $Y$ to a simplex coincides
with the set of fixed points $\Lambda(x)=(x;0;0;1-x),x\in[0,1].$

######  Lemma 3.

For any initial point $(x^{0},y^{0},u^{0})\notin M,$ there exists $k\geq 2$
such that $(x^{(k)},y^{(k)},u^{(k)})\in M.$

###### Proof.

Suppose the opposite, for any initial point $(x^{0},y^{0},u^{0})\notin M,$ and
for any $n\in N,$ $(x^{(n)},y^{(n)},u^{(n)})\notin M.$ It means that

(17) $ay^{(n)}-\beta x^{(n)}(u^{(n)}+qy^{(n)})\leq 0\ \ \mbox{or}\ \
bu^{(n)}-ay^{(n)}\leq 0.$

in addition, by the operator (3) one has

(18) $\begin{cases}y^{(n+1)}=y^{(n)}-[ay^{(n)}-\beta
x^{(n)}(u^{(n)}+qy^{(n)})]\\\\[5.69054pt]
u^{(n+1)}=u^{(n)}-[bu^{(n)}-ay^{(n)}]\end{cases}$

from the systems (17) and (18) we have that $y^{(n)}$ or $u^{(n)}$ increasing
sequence, but by Remark 2 both sequences converges to zero. This contradiction
completes the proof. ∎

######  Theorem 1.

Assume (7) is satisfied. Then for any initial point
$\lambda^{0}=\left(x^{0},u^{0},y^{0},v^{0}\right)\in S^{3}$ (except fixed
point) the trajectory has the following limit

$\lim_{n\to\infty}V^{(n)}(\lambda^{0})=(\bar{x};0;0;1-\bar{x})$

where $\bar{x}$ depends on parameters and initial point $\lambda^{0}$ and
$\bar{x}<\frac{ab}{\beta(a+bq)}.$

###### Proof.

Existence of the limit is known by Proposition 2. First, we consider the
intersection of two surfaces $ay-\beta x(u+qy)=0,bu-ay=0$ which are
represented in the invariant set $M$ (Lemma 2). Solving the system of
equations $ay-\beta x(u+qy)=0,bu-ay=0$ we have the parametric equations of two
direct lines $x=x,y=0,u=0$ and $x=\frac{ab}{\beta(a+bq)},y=y,u=ay/b.$ (as in
Fig.1).

Figure 1. Red lines are intersection of surfaces $ay-\beta x(u+qy)=0,bu-ay=0$.

Based on Lemma 2 and Lemma 3 we have that $\bar{x}<\frac{ab}{\beta(a+bq)}.$
Biologically it means that total number of susceptible at the limit does not
exceed $\frac{ab}{\beta(a+bq)}.$ In Fig. 1 dotted red line represents the
peaks of disease for initial points not in $M.$
Dependence of $\bar{x}$ to initial state and parameters can be represented as
following: from the last equation of the system (3) we find

$u=\frac{v^{(1)}-v}{b}\Rightarrow u^{(1)}=\frac{v^{(2)}-v^{(1)}}{b},$

substituting this into third equation we get

$y=\frac{v^{(2)}+(b-2)v^{(1)}+(1-b)v}{ab},$

using them and second equation of the system one has

$x=\frac{v^{(3)}+(a+b-3)v^{(2)}+(3+ab-2a-2b)v^{(1)}+(a+b-ab-1)v}{\beta(qv^{(2)}+(a+bq-2q)v^{(1)}+(q-bq-a)v)}.$

Consequently, using $x+y+u+v=1$ we can find the following recurrence formula
based on $v^{(n+3)},v^{(n+2)},v^{(n+1)},v^{(n)}$:

$abv^{(n+3)}=-\beta
q\left(v^{(n+2)}\right)^{2}-\beta(2bq-4q+a+aq)v^{(n+2)}v^{(n+1)}-$
$-\beta[q(1-a)(1-b)+(q-a-
bq)]v^{(n+2)}v^{(n)}-\beta(a+bq-2q)(a+b-2)\left(v^{(n+1)}\right)^{2}-$
$-\beta[(1-a)(1-b)(a+bq-2q)+(a+b-2)(q-a-bq)]v^{(n+1)}v^{(n)}-$
$-\beta(1-a)(1-b)(q-a-bq)\left(v^{(n)}\right)^{2}-ab[a+b-\beta q-3]v^{(n+2)}-$
$-ab[3+ab-2a-2b-\beta(a+bq-2q)]v^{(n+1)}-$ $-ab[a+b-ab-1-\beta(q-a-
bq)]v^{(n)},$

where for a given initial point $(x^{0},y^{0},u^{0},v^{0})$ using the system
(3) one can define $v^{(1)},v^{(2)}.$ This recursion gives dependence of
$v^{(n)}$ and its limit on the initial point. But the problem of finding of an
explicate formula for $\bar{x}=\bar{x}(\lambda^{0})$ remains open. ∎

##  5\. Discussion

In this work, we proved that after some time the COVID-19 virus disappears
(without counting mutation). But, if the disease transmission is large enough,
then the big part of population will be infectious. Using computer analysis we
have studied the disease spreading for concrete values of parameters. For
Uzbekistan, we assumed that incubation period is 10 days ($a=0.1$), infectious
period is 15 days ($b=0.066$), transmission rate is $\beta=0.12$, scaling
factor is $q=1$ and we get the following results:

(1) peak of the disease is 140 days from the start of the epidemic (real case:
138 days, August 4, 2020
333https://www.worldometers.info/coronavirus/country/uzbekistan/ ) Fig.2.

(2) completion of the disease is 300 days from the start of the epidemic (real
case: 290 days, January 5, 2021
444https://www.worldometers.info/coronavirus/country/uzbekistan/ ) Fig. 3.

Note that each coordinates of initial state must be multiplied to total
population size $N$ of a considering country. Finally, the model is working
for COVID-19 disease without counting coronavirus mutation, (i.e., new
stamps). One aspect of our future work is focused to apply this model and more
realized mathematical models for the COVID-19 stamps.

Figure 2.
$a=0.1,b=0.66,\beta=0.12,q=1,x^{0}=0.99999,u^{0}=0.00001,y^{0}=0.0,n=140\,(days)$.
Figure 3.
$a=0.1,b=0.66,\beta=0.12,q=1,x^{0}=0.99999,u^{0}=0.00001,y^{0}=0.0,n=300\,(days)$.

For any initial point (except fixed point) the trajectory located between
surfaces $ay-\beta x(u+qy)=0$ and $bu-ay=0$ after the peak of the disease
(Fig.4).

Figure 4. $a=0.1,b=0.66,\beta=0.12,q=1$.

In [10] we gave a prediction of COVID-19 pandemic in Uzbekistan. There our
method is an observation based projection similar than the classic Mooreís Law
in microelectronics [6].

## References

  * [1] A.Radulescu, C.Williams, K.Cavanagh, Managment strategies in a SEIR-type model of COVID 19 community spread. Scientific Reports, 10:21256 (2020).
  * [2] M. Martcheva, An Introduction to Mathematical Epidemiology. Springer, 2015, 462 p.
  * [3] R.L. Devaney, An Introduction to Chaotic Dynamical System. Westview Press, 2003, 336 p.
  * [4] R.N. Ganikhodzhaev, F.M. Mukhamedov, U.A. Rozikov, Quadratic stochastic operators and processes: results and open problems, Inf. Dim. Anal. Quant. Prob. Rel. Fields. 14(2) (2011), 279-335.
  * [5] Y.I. Lyubich, Mathematical structures in population genetics, Springer-Verlag, Berlin, 1992.
  * [6] J. Lu, A new, simple projection model for COVID-19 pandemic. https://doi.org/10.1101/2020.03.21.20039867.
  * [7] S.K. Shoyimardonov. A non-linear discrete-time dynamical system related to epidemic SISI model. To appear in Commun. Math. 2021\. Availabe in arXiv:2008.12677v1.
  * [8] U.A. Rozikov, S.K. Shoyimardonov, Ocean ecosystem discrete time dynamics generated by l-Volterra operators. Inter. Jour. Biomath. 12(2) (2019) 24 pages.
  * [9] U.A. Rozikov, S.K. Shoyimardonov, R.Varro, Planktons discrete-time dynamical systems. To appear in Nonlinear Studies, 28(2) (2021). Availbe in arXiv:2001.01182.
  * [10] U.A. Rozikov, S.K. Shoyimardonov, A Model for COVID-19 pandemic and prediction for Uzbekistan. Bulletin Inst. Math.2 (2020), 19–21.
  * [11] U.A. Rozikov,  Population dynamics: algebraic and probabilistic approach. World Sci. Publ. Singapore. 2020, 460 pp.
  * [12] O.Galor, Discrete dynamical systems. Springer, 2007, 158 p.
  * [13] G.Osipenko, Center Manifolds. In: Meyers R. (eds) Encyclopedia of Complexity and Systems Science. Springer, New York, 2009, 17 pp.