{'fraction_non_alphanumeric': 0.03910503692273779, 'fraction_numerical': 0.010691061390940152, 'mean_word_length': 4.376970487139979, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Data-driven control on encrypted data

We provide an efficient and private solution to the problem of encryption-aware data-driven control. We investigate a Control as a Service scenario, where a client employs a specialized outsourced control solution from a service provider. The privacy-sensitive model parameters of the client's system are either not available or variable. Hence, we require the service provider to perform data-driven control in a privacy-preserving manner on the input-output data samples from the client. To this end, we co-design the control scheme with respect to both control performance and privacy specifications. First, we formulate our control algorithm based on recent results from the behavioral framework, and we prove closeness between the classical formulation and our formulation that accounts for noise and precision errors arising from encryption. Second, we use a state-of-the-art leveled homomorphic encryption scheme to enable the service provider to perform high complexity computations on the client's encrypted data, ensuring privacy. Finally, we streamline our solution by exploiting the rich structure of data, and meticulously employing ciphertext batching and rearranging operations to enable parallelization. This solution achieves more than twofold runtime and memory improvements compared to our prior work.


I. INTRODUCTION
As cloud services become pervasive, inspiring the term "Everything as a Service" [1], privacy of the underlying data and algorithms becomes indispensable. Distributed systems coordination and process control enterprises can benefit from such services by commissioning cloud servers to aggregate and manipulate the increasing amounts of data or implement reference governors. For example, in the context of smart building automation, Control as a Service (CaaS) businesses [2] offer high-performance algorithms to optimize a desired cost, while achieving the required control goals. Clients like hospitals, factories, commercial and residential buildings might prefer such specialized outsourced solutions rather than investing in local control solutions. However, the data that is outsourced, stored and processed at a CaaS provider's cloud server consists of privacy-sensitive measurements such as user patterns, energy and temperature measurements, occupancy information etc. Nowadays, data can be monetized and used to disrupt the normal functionality of a system, leading to frequent database leakages and cyberattacks. This compels researchers and practitioners to design decision algorithms that are optimized for both efficiency and privacy. In this work, we provide a proof of concept for such a private cloud-based control algorithm.
In a client-server context, homomorphic encryption provides an appealing tool that allows one server to evaluate a functionality over the encrypted data of the client without the need to decrypt the data locally, at independently specified precision and security levels. In contrast, secure multi-party computation tools require collaboration between multiple cloud servers to execute the encrypted computations, and differential privacy introduces a trade-off between the accuracy and privacy levels. As discussed in [3], private single client computing can be enabled using homomorphic encryption schemes. Moreover, this technology has matured enough for it to be deployed in practice, especially in the healthcare [4] and financial [5] sectors. Since its genesis in [6], fully homomorphic encryption (FHE) has been developed substantially, in terms of more theoretically efficient leveled constructions [7]- [9], bootstrapping methods [10], [11], computational and hardware optimizations. There are many libraries that implement various FHE schemes and capabilities [12]- [15] and a tremendous number of papers that build on them.

A. Related work
Recently, private control algorithms operating on encrypted data have received a lot of interest [16]- [21]. These works can be classified with respect to the type of control algorithms they implement and the cryptographic tools they use. Linear control algorithms with public model and gains are considered in most works and are implemented using partially homomorphic encryption (PHE) schemes. Nonlinear control algorithms with public model are considered in [18], [19], where PHE is combined either with secure multiparty computation schemes or the client partakes in the computation. In [16], [20], somewhat or fully homomorphic encryption schemes are preferred, either to guarantee more privacy for the system parameters (that are still known at the client) or to achieve more complex control algorithms. Like our current work, [21] is based on data-driven techniques and uses leveled homomorphic encryption. There, the cloud computes the value function in a reinforcement learning task over multiple time steps in a one-shot way, while we collect encrypted data from the client at every time step and compute on it iteratively. Finally, our current work enhances and improves [22] in the following aspects: we propose a new encrypted solution that reduces more than twofold the time and the memory requirements of the previous solution, allowing us to simulate more realistic and more secure scenarios; we prove the close relationship between the standard data-driven LQR and our reformulation; and we present many more extensions.
As mentioned before, most private cloud-based control schemes assume that the model of the system is known to the client or public. However, the model might not be always available and has to be computed from data. Identification and data-driven control of unknown systems have been extensively studied in classical and recent literature [23]- [32]. These methods are usually designed with the objective of sample-efficiency and control performance, without considering a private implementation. For example, the standard system identification-certainty equivalence control architecture might require solving non-convex problems [33] or might involve Singular Value Decomposition (SVD) [24], which are prohibitive from an encrypted evaluation perspective. On the other hand, the behavioral framework [26]- [28], [30], [32], [34], which has received renewed interest recently, is more compatible with modern encryption tools and allows us to use less costly encrypted operations. In the behavioral framework, an alternative data-driven representation of linear systems is considered: if the collected input-output data are rich enough (input is persistently exciting), then the model is implicitly characterized by the collected data [26].

B. Contributions
Our goal is to perform online data-driven control on encrypted data, while maintaining the privacy of the client's uploaded input-output data, desired setpoint and control actions. In the context of private implementation, we need to depart from the typical techniques used in the non-private versions and make compromises between tracking performance and privacy. Our controller is based on the behavioral framework and the control performance is captured by the LQR cost. We consider two versions and adapt them to account for an encrypted implementation: i) an offline version, where precollected data from the system is encrypted and handed to the cloud to compute an offline feedback control law, and ii) an online version, where the cloud server computes the control at every time step based on both the offline precollected data and on the encrypted measurements received from the client. We employ a new leveled homomorphic encryption scheme optimized for efficiency and precision as the powerful base tool for encrypting and evaluating the private data on the untrusted cloud machine. At the same time, we rewrite the computation of the control input as a low-depth arithmetic circuit by using mathematical tools such as the Schur complement to facilitate matrix inversions, and carefully batching the private values into ciphertexts. Specifically, our contributions are the following: • Formulate a setpoint tracking problem that learns the control action from the available data at every time step and is amenable to encrypted implementation. • Prove the closeness between this formulation and a data-driven LQR problem. • Extend this approximation to allow the online collection and incorporation of samples to improve robustness and adaptation to new data and to errors arising from the use of encryption. • Design storage and computation efficient encrypted versions of the offline and online data-driven problems, by carefully encoding the values to facilitate parallelization and re-engineering the way the operations are performed. • Implement the algorithms for a standard security parameter and realistic client and server machines. Present the results for a temperature control problem and showcase the storage and complexity improvements compared to the previous algorithm in [22]. We emphasize that the techniques for performing the encrypted computations efficiently are one of main contributions of the paper. These tools can easily be generalized and used in implementing other related problems in an encryption-friendly way, such as learning and adaptive control algorithms, recursive least squares, sequential update of inverses etc.

C. Notation
For a vector x ∈ R n , we denote by x [i] the i'th element of the vector. We denote the L 2 norm of a vector x by x 2 and the spectral norm of a matrix M by M 2 . For a matrix P ∈ S n ++ the P-quadratic norm is x 2 P := x Px. The vector e i is the vector of all zeros, except for position i − 1 where it has 1. For a set of consecutive integers, we define [n] := {0, . . . , n − 1}. The notation E(·) denotes the encryption of a quantity, with more parameters to be specified in the text.

II. PROBLEM FORMULATION
Let x ∈ R n , u ∈ R m , y ∈ R p be the state, control input and measurement of a linear system in (1). A client owns system (1) and contracts a cloud service to provide the control for this system: as depicted in Figure 1. A private CaaS should solve the control problem that minimizes the cost associated to setpoint tracking while satisfying the unknown system's dynamics (1). Fig. 1: Client-server diagram for data-driven control.
Control requirements. Assume we are given a batch of offline input-output data for system (1). Then, at every time t, given u 0:t−1 and y 0:t−1 we want to compute in a private and receding horizon fashion u * ,t ∈ R N m in order to track the reference r t , for some costsQ,R, which is the solution of the LQR optimization problem: with no prior knowledge about the system model A, B, C. We also want (2) to perform well under unknown small process and measurement noise and bias.
In this work, we focus on investigating an approximation of problem (2) without hard constraints on the inputs or outputs of the system. In Section III, we describe a data-driven reformulation of (2) that is encryption-friendly. Privacy requirements. In the scenario we consider, the cloud service should not be able to infer anything about the client's private data, which consists of the input signals u, the output signals y, the model A, B, C and state x (the last four quantities being unknown in a data-driven control problem), and any intermediate values. The costsQ,R can be chosen by the cloud, as part of the CaaS service or chosen by the client. The cloud service is considered to be semi-honest, which means that it does not deviate from the client's specifications. We assume this to be the case because the cloud server is under contract.
The formal privacy definition is provided in Section VIII and can be written in the standard simulation paradigm [35,Ch. 7], i.e., all information obtained by the server after the execution of the protocol (while also keeping a record of the intermediate computations) can be obtained solely from the inputs and outputs available to the server. Efficiency requirements. The computation time of the control action at the current time step should not exceed the sampling time, which is the time until the next measurement is fetched. In the private case, we expect a large overhead for complex controllers and consider applications with sampling time of the order of minutes. Moreover, we require the client to be exempt from heavy computation and communication, and the bulk of the computation (in the allowable time) should be performed at the server side.

III. ENCRYPTION-AWARE LEARNING FORMULATION
We can formulate this control problem inspired by the behavioral and subspace identification frameworks explored in [26]- [28], [30], [32]. First, we introduce some preliminary concepts. A block-Hankel matrix for the input signal u = u 0 u 2 . . . u T −1 ∈ R mT is given by the following, for a positive integer L: Furthermore, the signal u is persistently exciting of order L if H L U ∈ R mL×(T −L+1) is full row rank. Let us construct block-Hankel matrices for the "past" and "future" input and output data, u ∈ R mT , respectively, y ∈ R pT , for M samples for the past data and N samples for the future data, and S := T − M − N + 1: A. Offline feedback data-driven control problem Assume that we are given precollected input and output data, with U p , Y p , U f , Y f the respective past and future Hankel matrices, for some past and future horizons M, N . Assume also that the precollected input is persistently exciting.
Assumption 1 (Data richness). We assume the offline precollected input trajectory u is persistently exciting of order M + N + n, where n is the order of the system [26].
Fix a time t and letū t = u t−M :t−1 be the batch vector of the last M inputs. The batch vector of the last M outputsȳ t is defined similarly. If M ≥ n, then the standard LQR problem (2) can be re-formulated as the following data-driven control problem [30]: where the state representation has been replaced with the precollected data and g represents the preimage of the system's trajectory with respect to the precollected block-Hankel matrices. Note that u * ,t [1:m] , the first m elements of u * ,t , will be input into the system in a receding horizon fashion, and y * ,t is the predicted output.
We now depart from the behavioral control problem (5) considered in the existing literature and explore a more encryption-friendly form.
First, we rewrite (5) as a minimization problem depending only on g by enforcing u = U f g and y = Y f g. Second, 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 constraint of (5). Hence, we prefer a least-squares penalty approach to the equality constraint in (5) with regularization weights λ y and λ u . Finally, to reduce overfitting and precision errors due to encryption, we also penalize the magnitude of g through two-norm regularization. We opt for two-norm regularizations to obtain better efficiency for the encrypted algorithm, as well as more robustness with respect to noise, and uniqueness of the solution g * ,t . This yields the following formulation, where Q = blockdiag(Q, . . . ,Q) and R = blockdiag(R, . . . ,R) and we reuse the notation r t for the batch reference signal: Note that the resulting problem (6) is an approximation of (5). Finally, (6) can be written as a quadratic program: where M ∈ R S×S is: Since (7) is a strongly convex optimization problem (ensured by the regularization term λ g I), we can find the optimal value for g by zeroing the gradient of the objective function: Going back to the control input for the current time step we want to compute, we obtain from (9): from which we select the first m elements and input them to the unknown system. As seen from (10), the controller has the form of a dynamic output-feedback law, where the feedback terms are computed using only the offline precollected data. More details about the choices we made in assembling this data-driven control problem formulation are given in Section IV and Appendix A.

B. Closeness between approximate and original problems
We prove that if we select small enough regularization coefficient λ g and large enough penalty coefficients λ y , λ u , then the solution of the approximate problem (6) is very close to the minimum norm solution of the behavioral control problem (5). The main challenge in the proof is that the past and future Hankel matrices are rank-deficient. Hence the behavioral control problem (5) involves singular matrices in both the objective and the constraints, and has multiple optimal solutions. This requires an involved analysis, where we deal with pseudo-inverses and subspaces, shown in Appendix A. To formally state the result we need two definitions. The set of optimal solutions of (5) is denoted by G opt : The minimum norm element of G opt is defined as Theorem 1. Consider the original behavioral control problem (5) and its approximation (6), with penalty coefficient λ y = λ u = λ > 0 and regularization coefficient λ g > 0. Let g min be the minimum norm solution of the behavioral problem (5) as defined in (12). Let g * be the optimal solution of the approximate problem (6). Then: as (λ g , λ) → (0, ∞) restricted on the set λ g > 0, λ > 0.
Although the behavioral problem has infinite solutions, due to the regularization term the approximate problem (6) will only return a solution close to the minimum norm one g min .

C. Online feedback data-driven control problem
To satisfy Assumption 1, it is necessary that the precollected input signal has length at least (m + 1)(M + N + n)−1, cf. [30]. In practice, Assumption 1 might be violated if less precollected data is available. Another issue with the precollected data is that it can be affected by perturbations, e.g., measurement noise. To alleviate these issues, we prefer an online algorithm, where the Hankel matrices H M +N U and H M +N Y are updated at each time step for anotherT steps with the used control input and the corresponding output measured (while keeping the block Hankel form). This approach empirically ensures richness of the data and robustness to perturbations. However, it is important thatT is not too large, in order to prevent overfitting.
Note that the precollected and the online data in this online algorithm belong to different trajectories. For this reason, we must compute the Hankel matrix for each data set separately, then append them in a single matrix [36]. This means that the online adaptation of the matrices U p , U f , Y p , Y f can only start at time t = M + N − 1, when we obtain enough data u 0 , y 0 , . . . , u M +N −1 , y M +N −1 to fill the first Hankel matrix column for the online data set. We call this phase between t = M − 1 and t = M + N − 1 trajectory concatenation.
The data-driven LQR algorithm is given in Algorithm 1.
Algorithm 1 Online data-driven control algorithm  Compute u * ,t = U f g * ,t and obtain (10). 8: Input to the system u t = u t, * [1:m] and measure the output y t .

9:
Updateū t andȳ t to be the last M components of ū t u t and ȳ t y t , respectively. The decryption primitive of a homomorphic encryption scheme is a homomorphism from the space of encrypted messages, or ciphertexts, to the space of unencrypted messages, or plaintexts. Homomorphic encryption allows one server to evaluate multivariate polynomial functionalities over the encrypted data of the client. An encryption 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. Leveled homomorphic schemes can be turned into fully homomorphic schemes by a bootstrapping operation. We direct the interested reader to the following surveys on homomorphic encryption [37], [38].
The common term for such a multivariate polynomial functionality is arithmetic circuit and the logarithm of the degree of the polynomial is the multiplicative depth of the circuit. Each operation evaluated on ciphertexts introduces some noise, which, if it overflows, can prevent correct decryption. Multiplications introduce the most noise, so we focus on the number of sequential multiplications allowed in an instance of a leveled homomorphic scheme, which we call multiplicative budget.
In this paper, we work with a leveled homomorphic encryption scheme. Specifically, we use the version of the CKKS scheme [9], optimized to run on machine word size of 64-bit integer arithmetic [39], [40] instead of multiprecision integer arithmetic. We give more details with respect to the primitives of this scheme in Appendix B. We chose this scheme because it can perform operations on encrypted real numbers with a smaller error than other leveled homomorphic schemes. Each real number is multiplied by a positive integer scaling factor and truncated, as commonly done, but the real advantage of this scheme is that one can remove the extra scaling factor occurring in the result after a multiplication, through a rescaling procedure, at very little error. This manages the magnitude of the underlying plaintexts, which could otherwise cause overflow in a large depth circuit.
A ciphertext's size grows with the number of sequential multiplications it supports. Each ciphertext, respectively plaintext, is characterized by a level, a scaling depth and a number of moduli. A new ciphertext (freshly encrypted, rather than obtained as the result of operations on other ciphertexts) is on level L, scaling depth 1 and has as many moduli as the multiplicative budget minus one. The number L minus the current level corresponds to the number of rescaling operations previously performed on the ciphertext and scaling depth corresponds to the number of multiplications without rescaling that have been performed. After a multiplication followed by a rescaling procedure, the number of levels decreases by one, the scaling depth remains the same and one modulus is dropped. To avoid some technicalities and parallel the notion of circuit depth, in the rest of the manuscript, when we refer to the depth of a ciphertext, we refer to the multiplicative depth of the arithmetic circuit used in order to obtain that ciphertext. If rescaling is done after every multiplication, then the depth will be the scaling depth minus one.
In the rest of the paper, we will analyze the multiplicative depth of a quantity x, obtained after evaluating an arithmetic circuit. Out of all circuits that output x on the same inputs, we choose to evaluate the one with the minimum depth d(x). In turn, this means that the ciphertext encrypting the resulting quantity E(x) will have d(x) + 1 consumed levels. Note that the fewer the number of levels necessary in the computation, the cheaper the operations on ciphertexts are.
In the CKKS scheme, each plaintext is a polynomial in the ring of integers of a cyclotomic field with dimension ringDim, which enables us to encode multiple scalar values in a plaintext and a ciphertext. Hence, we obtain ciphertexts that contain the encryption of a vector. This is a standard practice in lattice-based homomorphic encryption that allows performing single instruction multiple data (SIMD) operations, and can bring major computation and storage improvements when evaluating an arithmetic circuit. We can pack up to ringDim/2 values in one plaintext/ciphertext using a Discrete Fourier Transform [9]. Packing can be thought of as the ciphertext having ringDim/2 independent data slots. Abstracting the details away, the SIMD operations that can be supported are addition, element-wise multiplication by a plaintext or ciphertext and data slot permutations that can achieve ciphertext rotations (e.g., used for summing up the values in every slot). In what follows, we will use + and for SIMD addition and multiplication and ρ(x, i) to denote the row vector x rotated to the left by i positions (i < 0 means rotation to the right). However, once the values are packed, we should try to avoid extractions of individual values from the ciphertexts to improve computation time and circuit multiplicative depth.
In what follows, we denote by e v0 (x) the encoding of the vector x followed by trailing zeros in one plaintext and by e v * (x) the encoding of the vector x followed by junk elements (elements whose value we do not care about, usually obtained through partial operations). We denote by e vv (x) the encoding of the vector x repeated [x x . . .] in one plaintext. Finally, we also define the encoding of a vector repeated element-wise: for x = x [1] x [2] . . . x [n] , this encoding is e vr0 (x) = x [1] . . . x [1] x [2] . . . x [2] . . . x [n] . . . x [n] 0 0 . . . , where each element is repeated for a specific number of times. Similarly, we use the notation e vr * (·) when the repeated elements are followed by junk elements. When constructing ciphertexts through encryption, we use similar notations for the encryption of the corresponding encodings of a vector: A scheme has a security parameter κ if all known attacks against the scheme take 2 κ bit operations. For CKKS, the security parameter is determined according to the Decisional Ring Learning with Errors problem hardness [41]. We are also interested in the semantic security of this encryption scheme, which, at a high level, states that any two ciphertexts are computationally indistinguishable to an adversary that does not have the secret key.

V. CO-DESIGN OF ENCRYPTED CONTROLLER
According to the requirements in Section II, the challenges for the encrypted data-driven control can be summarized as: • The computations are iterative and not readily formulated as (low-depth) arithmetic circuits.
• The problem is computationally intensive: it requires large storage, large matrix inversions and many consecutive matrix multiplications. • The precision loss due to the private computations might affect the control performance. To deal with these challenges, we design an encrypted version of the closed-form solution (9) of the control problem stated in Section II and manipulate the computations involving matrix inverses to reduce the multiplicative depth needed. We employ the CKKS homomorphic scheme described in Section IV to address the precision challenge.
First, we approximate problem (5) into optimization problem (7), as shown in Section III, to simplify the encrypted computations and to avoid infeasibility due to noise and encryption errors. Second, we aim to write the computation of the solution of (7) as a low-depth arithmetic circuit. Despite being preferable to an iterative algorithm for solving the optimization problem (7) (that would increase the depth at every iteration), the closed-form solution (9) involves the encrypted inversion of a matrix, which cannot be generally written as a low-depth arithmetic circuit. This is not a problem in the offline feedback version of the problem (7), because this matrix inversion is required only once and such complex computations can be all performed offline, leaving only three encrypted matrix-vector multiplications to be performed at each time step, as shown in Section VI. However, in the online feedback version of the problem, this matrix inversion and many consecutive matrix-vector multiplications are required at every time step. Nevertheless, we leverage the special structure of the matrix to be inverted by using Schur's complement [42, Ch. 0] to reformulate the inverse computation as some lower multiplicative depth matrix-vector multiplications and one scalar division. Furthermore, to avoid performing the division, which is costly on encrypted data, the server sends to the client the denominator such that the client can return its inverse to the server, turning division into multiplication.
We also perform further valuable optimizations in terms of ciphertext packing, adding redundancy in the encoded values, and manipulating the computations to be performed in a SIMD mode, which are crucial to the tractability of the solution.
The multiplicative depth of the arithmetic circuit computing u * ,t [0:m−1] can be reduced through reordering of the intermediate operations. Given a circuit that computes x and a circuit that computes y, with multiplicative depth d(x) and d(y), the multiplicative depth of the circuit that evaluates them in parallel and then computes their product is: d(xy) = max(d(x), d(y)) + 1. The multiplicative depth should not be confounded with the number of multiplications. Judiciously choosing the order in which to perform the multiplications can reduce the multiplicative depth of the result. For example, consider we want to compute We will use this trick in Section VII.
Finally, as explained in Section IV, each ciphertext is created with a number of moduli corresponding to the depth of multiplications it can support before the underlying plaintext is corrupted by noise. After each multiplication and rescaling, a modulus is removed from the ciphertext, leading to more efficient computations. We exploit this fact in our solution, by making sure we compute each operation on ciphertexts that have the minimum required number of moduli. This happens by directly encrypting the message at a lower number of moduli and by removing the number of moduli on the fly.

VI. OFFLINE FEEDBACK ENCRYPTED SOLUTION
Recall the closed-form solution of the optimization problem (7). We now compute the multiplicative depth of the quantities of interest for consecutive time steps. First, the quantitiesū 0 ,ȳ 0 , the reference signal and measurements are freshly encrypted at every time step t, hence they are inputs to the circuit that computes u t and have a multiplicative depth of 1. This means, Second, we assume that all the quantities obtained offline will have multiplicative depth 0. We do not address the offline computations in this paper; since these computations depend only on the offline data and are one-time, expensive secure solutions can be used, e.g., the cloud could perform the encrypted computations, then perform bootstrapping to refresh the ciphertexts [10], [11]. Hence, the cloud has fresh encryptions of the following products: and the underlying messages are inputs to the circuit and they have multiplicative depth 0. Then: For t = 0, we obtain d(u 0 ) = max(d(r 0 ) + 1, d(ȳ 0 ) + 1, d(ū 0 ) + 1) = 1. Generalizing: At time t + 1,ū t+1 will be updated by u t . We next present three encrypted methods of computing (14), that exhibit trade-offs in depth, memory and communication.
If we individually encrypt each element of the quantities in (14) in a separate ciphertext (the inputs to the circuit computing u t will be all elements, rather than three vectors and three matrices),ū t+1 will have the multiplicative depth of u t and d(u t ) = t + 1, but memory-wise, a step will require (m + 1)(pN + pM + mM ) ciphertexts.
It is more efficient from both storage and computation points of view to encode a vector instead of a scalar in a ciphertext and perform the matrix-vector multiplications by a diagonal method. This vector encoding requires a different analysis.
Diagonal method [43], [44]. We explore a method for efficient matrix-vector multiplication where the diagonals of the matrix are encrypted in separate ciphertexts (the notion of diagonal is extended for rectangular matrices).
1) Tall matrix: Consider we want to multiply a tall matrix S ∈ R u×v , u ≥ v, by a vector p ∈ R v . To this end, we extract the extended diagonals d i of S such that d i will have as many elements as the number of rows u, for i = 0, . . . , v − 1. The corresponding result q = Sp is computed as follows, wherep is p concatenated with itself as many times such thatp ∈ R u : This way of computing the encrypted matrix-vector multiplication is highly efficient when using ciphertext batching. Specifically, we encrypt each extended diagonal of S in one ciphertext, which gives us a storage of only v ≤ u ciphertexts and we encrypt p repeatedly in one ciphertext. The encrypted result is then computed as: 2) Wide matrix: In the case of a wide matrix S ∈ R u×v with u ≤ v, instead of using extended diagonals, we use reduced diagonals. One reduced diagonald i , for i = 0, . . . , v − 1 will have as many elements as the number of rows u. The storage is not as efficient since we have to encrypt v ≥ u ciphertexts for the diagonals. Apart from this, the computation is identical to (16) and produces E v0 (q).
We can improve the storage (at the cost of some extra rotations) when the number of columns divides the number of rows. In that case, we again extract the extended diagonals d i ∈ R v for i = 0, . . . , u − 1. The computation of the matrix-vector multiplication takes the following form: When u << v, it is useful to append rows of zeros such that u divides v, and then use the above method. Let us return to (14). In general, A r , A y , A u are wide matrices, but there can exist particular cases when m >> p or when we want to compute all u * ,t , leading to tall matrices.
We assume that at the onset of time step t, the cloud server has i.e., each extended diagonal (whose exact definition varies with the shape of the matrix: tall or wide) of the matrices A r , A y , A u is encrypted in a separate ciphertext. In order to be able to use the diagonal methods, the cloud server first obtains E vv (ū t ), E vv (ȳ t ), E vv (r t ) by rotating and adding E v0 (ū t ), E v0 (ȳ t ), E v0 (r t ). This vector packing substantially reduces memory: 3 + max(m, pN ) + max(m, pM ) + mM ciphertexts in the worst case and 3 + min(m, pN ) + min(m, pM ) + m in the best case, depending on the shape of the matrix and the type of diagonal chosen, but will require an extra level. From these quantities, the cloud server computes and sends back to the client one ciphertext containing E v0 (u t ) or E v * (u t ), depending on the type of diagonals chosen. Assume the cloud server obtained E v0 (u t ). After that, the cloud server has to create E vv (ū t+1 ) from E v0 (ū t ) and E v0 (u t ), so it: , the last part appearing because of the rotation to the left; • masks the result in order to truly obtain E v0 (ū t+1 ) (otherwise, the summing and rotation would not produce the correct result); • repeatedly rotates and adds to obtain E vv (ū t+1 ). Because of the CKKS encoding through the Discrete Fourier Transform, (which enables SIMD multiplications, unlike other possible encodings), the masking operation needs to consume a level in order to preserve precision [9]. Obtaining E v * (u t ) instead of E v0 (u t ) does not change the analysis, since the same masking applied in the third step also removes the junk elements in E v * . Furthermore, the analysis (not shown here) of other less efficient matrix-vector multiplication methods, such as row and column methods, also require an extra masking.
For t ≥ 1, equation (15) can be updated the same way E vv (ū t+1 ) is updated. However, it is the same cost for the client to encrypt E v0 (y t ) and E vv (ȳ t+1 ), so it can encrypt and send the latter. Whenever the allocated multiplicative budget is exhausted, the server can ask the client to send a fresh encryption ofū t , at little extra cost.
Nevertheless, a reasonable and inexpensive option is to ask the client to send along with the encryption ofȳ t a fresh encryption ofū t at every time step. This implies that d(u t ) = 1, i.e., a multiplicative budget of only 1 is required for computing the control input for no matter how many time steps. In constructing the algorithm in Section VII, we prefer this option for the trajectory concatenation phase.

A. Computing the solution using arithmetic circuits
Being capable of evaluating polynomials via homomorphic encryption theoretically gives us the possibility of evaluating any function arbitrarily close (using Taylor series for example). In practice, the multiplicative depth and loss of precision of these approximations limit the types of functions we can evaluate. Consequently, division is still prohibitive. Furthermore, large matrix inversion should be judiciously performed, in order to avoid computing products of as many factors as the number of rows.
In this section, we will use the following shorter notation: where u p,f , y p,f are vectors added at the end of the Hankel matrices U p,f , Y p,f . We will drop the time subscripts t for conciseness, using prime to denote the next time step. Let us investigate the inversion of matrix M. For step t = 0, all the values involved in computing u * in (10), are precollected and can be computed offline. This includes the inverse M −1 and other matrix products. However, at the next time step, cf. line 13 in Algorithm 1: HU := HU hu , HY := HY hy (17) The last m elements on the last column of HU , respectively of hu, are the values of u t at the previous time step. Then, the matrix M ∈ R (S+1)×(S+1) : is a rank-1 update of matrix M. Let Specifically, M will have the following form: Schur's complement [42,Ch. 0] gives an efficient way of computing M −1 from M −1 (assuming M −1 exists), by inverting a scalar s and computing a few matrix-vector multiplications: To avoid performing the division 1/s on encrypted data, the server will ask the client to perform it on plaintext and send back an encryption of result.

B. Reducing memory and depth of the arithmetic circuit
In this paper, for simplicity, we consider diagonal cost matrices Q and R. Otherwise, the multiplication by these matrices requires more complicated encrypted operations.
At a given time, the cloud server has (in an encrypted form) HU, HY,ū,ȳ and M −1 , along with the unencrypted costs Q, R, penalties λ g , λ u , λ y and reference signal r. The cloud has to compute the equivalent formulation of (10): There is a vast number of parameters to tune in the encrypted implementation, such as the packing method from the messages into the plaintexts, the storage redundancy, the order of performing the operations and choice of refreshing some ciphertexts. The trade-offs that these different versions bring are in terms of ciphertext storage, key storage (especially permutation key storage), type and number of operations, precision at the end of the operations and total multiplicative depth of the resulting circuit; these goals are intricately intertwined. For example, designing the circuit to have a lower multiplicative depth reduces the size of the ciphertexts and reduces the encryption load at the client, but might involve storing more ciphertexts and performing more computations at the server, compared to a version with a higher multiplicative depth. Thus, the main difficulty in making the computations tractable is to astutely batch-encode the vectors and matrices into ciphertexts in order to reduce the memory, and manipulate the operations in order to reduce depth, number of operations and storage.
We describe here how to change the flow of operations in (20)- (23) in order to minimize the number of ciphertexts representing the relevant quantities: HU, HY, M −1 , m, µ,ū,ȳ, u, and the number of operations, while keeping the depth to a minimum. We make use of the feature of lattice-based homomorphic encryption schemes of encoding a vector of values into a single plaintext, which is then encrypted in a single ciphertext. Since we are also dealing with matrices, we explored several options of how to encode the matrices in ciphertexts: columns, rows, (hybrid) diagonals, vectorized matrix (see [44]). Note that the same linear algebraic operation, implemented for different encodings, leads to a different number of stored ciphertexts, multiplicative depth of circuit, and number of SIMD operations. For the problem we tackle in this paper, we found that encoding each column of a matrix into a separate ciphertext minimizes the number of stored ciphertexts and number of operations, while keeping the same depth as when each element is encoded in a different ciphertext and the usual (unencrypted) method of performing the operations is used, as stated in Proposition 1. However, not all the matrices' columns are encoded in the same fashion, as we will see next. In summary, having some redundancy inside the ciphertexts, i.e., having some values encoded multiple times, helps optimize the depth and complexity.
The following proposition summarizes the complexity of the online feedback algorithm, when carrying out the steps in this section. The proof of Proposition 1 is constructive and is given in the remaining of this subsection and Appendix C.
Proposition 1. We evaluate the arithmetic circuit for the online feedback algorithm on encrypted data corresponding to Algorithm 1 for one time step with O((S + t) 2 ) operations, O(S + t) ciphertexts and O(S + t) rotation keys, at depth 2t + 4, where S is the number of columns of the offline generated block-Hankel matrix and t is the number of online samples accumulated so far. Figure 2a shows how to obtain an inner product between two encoded vectors using SIMD multiplications, rotations and additions. Notice that due to rotations, the resulting vector will have the relevant scalar in its first slot, and junk (partial sums) in the following slots. Figure 2b shows the most efficient way (in terms of required number of operations) to obtain the product between a matrix and a vector when the matrix is encoded as separate columns. Specifically, we need the elements of the vector to be repeated and separately encoded.  Let us look first at M −1 m . In order to use the efficient method outlined in Figure 2b, we need to have a separate ciphertext that encodes each element of m, repeated as many times as the number of columns in M −1 . In turn, this suggests that we should use the method outlined in Figure 2c when computing m as in (20), rather than the method outlined in Figure 2a which would require extra masking and rotations afterwards. Specifically, we have to repeat each element of hy for S times, encode the resulting vector in a ciphertext and do the same for each column of HY and the elements in the diagonal matrix Q and λ y . Then, we perform an element-wise multiplication between the ciphertext encoding hy and the ciphertext encoding the elements of Q and λ y , and then perform the inner product with the encoded columns of HY. The same steps are taken for hu, HU, R and λ u .

Remark 1.
To account for all the time steps where new samples will be collected, whenever we encode values repeatedly, we need to encode from the beginning S +T copies of the elements, where S is the initial number of columns in M andT is the total number of new samples intended to be collected. DefineM := (mM −1 ) (mM −1 ). Then: ∀i = 0, . . . , S − 1: . . .
The next piece is to compute s as follows: The server asks the client to invert s. After the server receives from the client 1/s, it can compute 1 sM , as required in (23), by replacing e i in (25)  For that, we also encodeū andȳ with repeated elements, so we can compute the inner product as in Figure 2c, for the elements of Z in (24) stored repeatedly. We also prefer some redundancy in storage to avoid online processing, by having a ciphertext for the first m rows of U f , encoded without repetition. This allows us to compute, for i = 0, . . . , m − 1: To send only one ciphertext for the result u (instead of m ciphertexts), while avoiding an extra masking, the server will pack the vectors υ i (which have an encoding with trailing zeros, whereas the encoding of u [i] does not) into one ciphertext and ask the client to perform the corresponding summation. The cloud service performs the encrypted version of (20)-(23), taking the steps outlined in this subsection, which we describe in detail in Appendix C in equation (55)-(57), along with the specific encoding of each quantity.

C. Precision discussion
Scaling: As described in Section IV and Appendix B, every encrypted operation introduces some noise that corrupts the least significant bits of the encrypted values, which accumulates with the depth of the circuit. This implies that values with a smaller magnitude are affected more by this noise. In the problem we address, the smallest values are found in M −1 t , hence scaling it by a positive integer factor α and accordingly changing the computations such that at every time step we obtain αM −1 t and αu t (such that the client can remove the scaling) allows us to increase the precision of the result. In order to incorporate this scaling without increasing the depth, the server has to perform more operations than before, e.g., recomputation from scratch of some quantities. The scaling factor should be selected by the client, and if its magnitude is sensitive, it can be sent as a ciphertext to the cloud server. We can modify the circuit such that this scaling does not affect the depth.
Precision loss when computing the Schur complement: Given the particularity of the matrix that has to be inverted, at each time step, some precision bits of the Schur complement are lost (depending on the regularization parameter and the magnitude of the measurements), which incurs a loss in precision in the subsequent computations.
Without the regularization term λ g I added to the objective function of (6), the resulting matrix M does not have full rank in the noiseless case (noise helps, but the inverse is not numerically stable). This suggests that the Schur complement (22) will have the same order of magnitude as λ g , regardless of the values in M, which we show in the following. Let H := HY HU and h := hy hu be a new column obtained by adding new samples, as described in Section VII-A, and P := blkdiag(λ y I, Q, λ u I, R). Then, We now prove that, in the noiseless case, the Schur complement has the same order of magnitude as λ g . Formally: Lemma 1. The following statements characterize the relation between s and λ g : The proof is given in Appendix C-B. The term in (b) h P 1/2 (H P 1/2 ) † (P 1/2 H) † P 1/2 h is small for slowly varying systems. Under small random noise, the results in Lemma 1 are not exact, but empirically follow closely.
This analysis is important from an encrypted implementation perspective. Lemma 1 shows that, for fixed cost values in P but regardless of the values in H and h, i.e., the values of the y measurements and of the u control actions, s will be close to the regularization parameter λ g . In the encrypted computations, we will compute µ = λ g + h Ph and mM −1 m each with a fixed precision of x bits. For large values of the measurements, we obtain that s ≈ λ g << h Ph, which means that there will be a cancellation in the most significant bits (e.g., the first y MSBs), leading to a loss of precision in s, which will now have only x − y bits of precision. This is particularly crucial since such a cancellation of the MSBs will happen at every time step we accumulate new data samples, leading to a cascading loss of precision. Furthermore, due to using the Residue Number System implementation of the encryption scheme, which improves the efficiency by orders of magnitude by using native sizes of 64 bits to store the "residue ciphertexts" rather than multi-precision arithmetic of arbitrary size, the precision cannot be increased indefinitely (see [39], [40] for technical details).
As a side note, decreasing or increasing the values in the cost matrix P independently from λ g does not solve this precision issue, because it implies changing the problem (it is equivalent to increasing or decreasing λ g ).
The conclusion of this discussion is that, apart from the role it has for regularization and noise robustness, the parameter λ g also affects the precision of the solution (meaning the difference between the encrypted solution and the unencrypted solution), which suggests we need to pick λ g to not be too small. This introduces a precision versus accurate convergence trade-off. We will comment more on this trade-off in Section X-B.

D. Considerations for continuous running
Accumulating new samples serves to robustify the algorithm and adapt to new disturbances. However, an important caveat is that adding too many new samples damages the performance of the algorithm, both because of overfitting and because of the problem becoming intractable as the number of variables in (7) grows. In future work, we will investigate this issue more, along with having a sliding window of sample collection.
Nevertheless, there are cases where the number of new samples to be collected leads to a circuit of a depth higher than the preselected multiplication budget. Our options for continuing the computations are: (i) Restore the initial precollected Hankel matrices which bypasses the refreshing step altogether. Advantages: no extra computations needed. Disadvantages: this causes oscillations in the control actions.
(ii) Stop adding new information to the matrices after the multiplication budget is exhausted. Advantages: no extra computations needed. Disadvantages: the multiplicative budget has to be large enough such that enough samples are collected.
(iii) Pack the matrix into a single ciphertext as described in (28) and ask the client to refresh it.
Advantages: the server can continue collecting values for any desired time, without extra multiplication depth. Disadvantages: the client has to decrypt, encrypt and send another ciphertext; the rotation keys necessary for packing can occupy a lot of storage (but compared to [22], here we only need S + t such rotations, not (S + t) 2 /2).
(iv) Bootstrap the ciphertext of M −1 t . The computation advancements [11] regarding the bootstrapping procedure suggest that it is likely to locally resolve the refreshing step. Advantages: the server can continue collecting values for any desired time without the client's intervention. Disadvantages: the initial multiplication budget has to be larger to also allow for the bootstrapping circuit.
In the solution we implemented, we chose option (iii). This gives us flexibility on the maximum multiplicative depth of the circuit, which will now be 2t ref resh + 4, at little extra cost for the client.
We can pack as many values as half the ring dimension ringDim in one ciphertext. Then, the number of ciphertexts the server can pack M −1 t into is (S +t)(S +t)/(ringDim/2) . This is viable since the ring dimension is in general large (e.g. ≥ 4092), in order to accommodate a reasonable multiplicative budget, plaintext precision and standard security parameter. The client only has to decrypt, re-encrypt and send back this number of ciphertexts. The server then uses one extra level to perform the reverse of (28) . These multiplications can be absorbed in the same initial multiplicative depth.
VIII. PRIVACY OF THE ENCRYPTED ALGORITHMS Definition 1. (Client privacy w.r.t. semi-honest behavior of server) Let x C be the private input of a client and x S be the input of a server. The client wants the server to evaluate a functionality f and return the result f (x C , x S ). We will denote by f S (x C , x S ) the corresponding result at the server. Let Π be a two-party protocol for computing f . The view of the server during an execution of Π on the inputs ( is (x S , r, msg), where r represents the outcome of the server's internal coin tosses, and msg represents all the messages it has received from the client. For a deterministic functionality f , we say that Π privately computes f w.r.t. to the server if there exists a probabilistic polynomial-time algorithm, called simulator and denoted S S , such that, for any inputs x C , x S : It is common for a server to be only employed for computing a result that is then sent to the client, and the server has no output formally, i.e., f S (x C , x S ) = ∅.
In our case, the functionality f in Definition 1 represents the data-driven algorithm, and f S (x C , x S ) represents the output of the server after evaluating the functionality on the given inputs. Fix a number of time steps T . We want to prove a secure evaluation of the data-driven control functionality for T steps, after which we assume the protocol ends. In reality, because the security of the encryption scheme used is based on computational problems, after a long time (years), it is recommended that the secret key is changed; so T represents a natural stopping point and is not restrictive.
Under the aforementioned condition, the output of the protocol at the server, f S (x C , x S ), is actually the empty set. Furthermore, all the intermediate messages in the view of the server will be ciphertexts encrypted with the client's public key using the CKKS scheme. Finally, the ring dimension of the CKKS scheme is selected such that, for the multiplicative budget of the desired functionality, an adversary cannot brute force the encryption. These are the main observations that will be used in proving the subsequent theorems.
There is a preprocessing protocol, where the server gets and computes the encrypted information it will have to use for the online control of the system, i.e., the encrypted data trajectories, costs and parameters. This preprocessing protocol is proven secure under the same Definition 1, hence, there is an ideal functionality F pre , which is essentially indistinguishable from the real-world functionality of the preprocessing protocol. We work in a F pre -hybrid model, which means that we can securely compose this functionality with the functionality of the online protocols (see [45] for more details about simulation proofs in the hybrid model and composition theorems). Theorem 2. The encrypted offline data-driven control algorithm in Section VI achieves privacy with respect to the server.
Proof. At time t ≤ T , the client's input is HU, HY, u 0:t−1 , y 0:t−1 , r 0:t−1 and pk, where pk is the public key of the leveled homomorphic scheme used. We assume that, after the preprocessing protocol, the server has input Q, R, λ g , λ u , λ y , pk and receives as messages E v0 (diag i A r ), E v0 (diag i A y ), E v0 (diag i A u ). At time t, the server receives from the client either E vv (ū t ), E vv (ȳ t ), E vv (r t ) or E vv (ȳ t ), E vv (r t ). After applying the functionality as described in Section VI, the server obtains E v0 (u t ) or E v * (u t ) and sends it to the client.
The proof immediately follows from the fact that the server only receives fresh ciphertexts from the client. Since the homomorphic encryption scheme is semantically secure, these ciphertexts are computationally indistinguishable from random encryptions, so we can construct a simulator for the server that replaces the true messages from the client by random encryptions of the same size. The input-output distribution of such a simulator will be indistinguishable from the input-output distribution of the true protocol.
Theorem 3. The encrypted online data-driven control algorithm in Section VII achieves privacy with respect to the server.
Proof. The input of the client is the same as in the proof of Theorem 2. The input of the server is Q, R, λ g , λ u , λ y , pk and after the preprocessing protocol, it has the corresponding ciphertexts for HU 0 , HY 0 , M −1 0 , U f 0 . At time t < T , while collecting new samples, the server has ciphertexts corresponding to HU t−1 , HY t−1 , M −1 t−1 ,ū t−1 ,ȳ t−1 , receives from the client E vr0 (u t ), E vr0 (y t ), E vr0 (r t ), then computes E vr * (s t ) and receives from the client E vr0 (1/s t ). The analysis in Section VII-C does not reveal private information: although the cloud server knows λ g , it never gets s t in cleartext. Finally, after applying the functionality described in Section VII and Section C, the server obtains the corresponding ciphertexts of HU t , HY t , M −1 t , U f t ,ū t ,ȳ t and E v * (υ t ), and sends the latter to the client, which decrypts it and computes u t . During the times designated for refreshing M −1 t , the server packs it into a single ciphertext, sends it to the client and receives E v0 (M −1 t ), which uses new uncorrelated randomness. After the collection of the new samples ends, i.e.,T ≤ t ≤ T , the server only updates the ciphertexts corresponding to the quantitiesū t andȳ t , but otherwise computes the same output.
Since the server only receives fresh encryptions of the private data, we can construct a simulator for the server that replaces the true messages from the client by random encryption of corresponding size and apply the server's functionality. Privacy follows because of the semantic security of the underlying CKKS scheme.

IX. EXTENSIONS
We now briefly describe extensions and possible modifications to the online algorithm.

A. Inputting blocks of control inputs
For simplicity of exposition, we illustrated the encrypted computations outline for when the server sends only the first component of u * ,t to the client. In practice, it is common that more consecutive components from the N computed control inputs are applied. This is also beneficial for the encrypted solution, since the most expensive computations, the inversion of M t and its subsequent update, need to be computed only once for multiple time steps. Similarly, the communication rounds between the client and the server would be reduced. We recommend this for larger systems.

B. Batch collection of new samples
We can still use the Schur complement to compute the inverse of a rank-x update of the matrix M t , where x is how many new data samples we want to accumulate at once. The server then packs the matrix representing the Schur complement into one ciphertext and asks the client to invert it and send it back. The circuit depth will increase accordingly.

C. Sliding window of sample collection
In order to manage the growth in the number of variables, i.e., number of columns of HU t and HY t , after collecting a sample and updating M −1 t , we would like to remove the first sample collected. We can again achieve this by using Schur's complement and the Woodbury matrix inversion lemma. We revert to the notation used in Section V for simplicity of the exposition and we show how to remove the first sample, i.e., the first row and first column of M (for removing rows and columns inside the matrix, we will have to multiply by permutation matrices). This gives the expression in (23) and the following equal expression for (M ) −1 : By the matrix inversion lemma, Removing a sample after collecting one sample incurs an increase by two in the total depth of the computation and asking the client to perform the division of l 1 .

D. Reducing accumulated noise
As discussed in Section VII-C, noise accumulates in the stored ciphertexts of M −1 t at the server, because of the fixed precision imposed by the native dimension of 64 bits (in the Residue Number System implementation of the encryption scheme). However, before the noise grows too much, the server could try to keep track of it and reduce it by using the fact it can also easily compute M t as in (18). Denote the error on M −1 t by X. Then: It is safe to assume ||X|| ∞ ≤ ||M t || ∞ (other norms can be chosen), with < 1. Furthermore, if satisfies the relation ||M t || 2 ∞ ≤ 1, then ||XM t X|| ∞ ≤ ||X|| ∞ . This would mean that the server can compute a better approximation of M −1 t than M −1 t + X as M −1 t + XM t X. However, this second approximation would imply a larger depth for M −1 t , so the best option is for the server to compute this approximation right before it asks the client to refresh the ciphertext.

E. Function privacy
The algorithm run at the CaaS service provider might be proprietary. In that case, the client should not get extra information about that algorithm, apart from what it agreed to receive. This notion is captured by the function privacy property of a homomorphic encryption scheme, which states that a ciphertext should not expose information about the function that was evaluated in order to create that ciphertext. Inherently, the CKKS scheme does not preserve function privacy because of the partial information stored in the junk elements. The common ways to achieve function privacy is for the server to mask out the irrelevant slots before handing a ciphertext to the client to decrypt, either by a multiplicative mask (zeroing out the junk elements) or by adding large enough randomness to them. The multiplicative mask increases the depth by one, whereas the randomness masking can lead to an increase in the plaintext modulus since the randomness has to be large enough to drown the relevant information.
In our scenario, the client is allowed to obtain the numerator and denominator of the control input at every time step and also the inverse matrix M −1 , which is only determined from the inputs and outputs that the client knows. The server can zero out the extra slots in E vr * (s t ) and E v * (υ t ). This masking would add one extra level to the respective ciphertext, but the level increase would not build up, since these masks are not used in posterior computations. Further investigation on keeping the service provider's algorithm and costs private from the client will be subject to future work.

X. IMPLEMENTATION AND EVALUATION
We focused on testing out our proposed algorithms on a temperature control problem, because first, smart buildings can support flexible encryption methods for the collected data, and second, the sampling time is of the order of minutes, which allows for complex computations and communication between the client and the server to take place.
A. Comparison with algorithm in [22] We first want to underline the improvements in both runtime and memory that the algorithm presented in the current work brings compared to the algorithm in [22]. For the same parameters considered in the zone temperature problem exemplified in [22], which had one input and one output, used a ring dimension of 2 12 and 21 levels in the ciphertext, we used the same commodity laptop with Intel Core i7, 8 GB of RAM and 8 virtual cores at 1.88 GHz frequency. We obtained an improvement of 2.4x in terms of memory: 2.25 GB compared to 5.48 GB. We obtained an improvement of 4x in terms of total running time, and specifically, the maximum runtime per step dropped from 308.8 s (previous algorithm) to 75 s (current algorithm). Some of the changes (e.g., in terms of better thread parallelization) we applied to our current algorithm can also be applied to the algorithm in [22]. These improvements allowed us to simulate the encrypted control for larger systems with better security parameters, as described in the next parts. B. 2x2 system with 73 bits of security 1) Setup: We considered a temperature control problem for the ground floor of a building with two zones. The sampling time is of T s = 420 s. The system parameters are: n = 2, m = 2, p = 2, M = 4, N = 4, T = 32 and the system is stable. To mimic a realistic example, we considered two types of disturbances: unknown, characterized by a zero mean Gaussian process noise with covariance 0.001I and a zero mean Gaussian measurement noise with covariance 0.01I, and known, slowly varying disturbances caused by the exterior temperature and the ground temperature, which we simulate as varying uniformly at random between [24.5 • , 25.5 • ] and [9.5 • , 10.5 • ], respectively. We choose the cost matrices and regularization terms Q = I, R = 10 −3 I, λ g = 5, λ y = λ u = 1.
For the offline data collection, we assume a different initial point than for the online computation, and randomly sample the offline input trajectory so that the corresponding output trajectory lies in the interval [10 • , 40 • ]. The M, N and T parameters mean that we start from 25 columns in the Hankel matrices. We note that, due to noise, the offline feedback control algorithm has suboptimal performance for these parameters-see Figure 3. Thus, we collect data online for 15 more time steps, which means adding another 15 columns to the Hankel matrices. Afterwards, we run the system with fixed gains. Figure 3 shows the performance of the setpoint tracking with these parameters and under noise and disturbances.
2) Implementation details: For better benchmarking capabilities (less variability between runs), for these experiments we used an AWS EC2 c5.2xlarge machine, with 8 virtual cores, 3.4GHz frequency and 16 GB of RAM. We implemented 1 the encrypted solutions proposed in this paper using the PALISADE library [14].
We use a ring dimension of 2 14 and 15 levels (for a refresh step after 5 collected samples). The resulting ciphertext modulus of a fresh ciphertext is of 760 bits (the first plaintext modulus has 60 bits and the following plaintext moduli have 50 bits). This gives a security size of 73 bits, according to the LWE estimator [41]. Choosing a smaller refreshing time increases the communication rounds between the client and the server (one extra exchange at every 35 minutes in this case), but at the same time decreases the multiplicative depth of the circuit, reducing the computational complexity and the size of the parameters. 3) Precision: Figure 3 shows the tracking performance for the encrypted and unencrypted versions of the online and offline feedback algorithms. Due to the precision loss incurred by working with a deep circuit on encrypted data, we see a small offset in the measured output compared to the unecrypted measured output for the online feedback algorithm. Since the offline feedback algorithm is a very shallow circuit (depth 1) when the client returns a fresh encryption of the control input to the server, it incurs no precision loss.
We obtain a maximum magnitude of 0.129 • C and an average of 0.045 • C for the difference between the unencrypted measured output and the encrypted output during the online simulation of a noisy process. For the control input, the maximum difference is 0.248 kWatts and an average difference of 0.053 kWatts. A typical such sample is depicted in Figure 3. Notice that the maximum magnitude of the difference is smaller than the variation in the disturbances and does not affect the performance and stability of the control algorithm. One can increase the plaintext modulus in order to reduce the error introduced by encryption (e.g. 53 bits instead of 50 bits in the moduli will lead to an almost tenfold decrease in the differences). However, if we keep the ring dimension constant, this reduces the security level to 69 bits. In Section X-C, we show results for better precision at a standard security level. As described in Section VII-C, the loss of precision can be alleviated by a larger regularization parameter λ g . For example, keeping all other parameters the same but increasing λ g = 10 will lead to a maximum magnitude of the difference in the measurements of 0.037 • and the maximum magnitude of the difference in the inputs of 0.043 kWatts, but the tracking performance is slightly reduced (from 0.06 • C to 0.12 • C deviation from the reference in the noiseless case).

4) Offline:
The encrypted offline feedback algorithm is very fast, since it has only depth 1 when the client sends to the server encryptions ofū t andȳ t . A ring dimension of 2 12 gives a security parameter of 126 bits, with a ciphertext having a size of 110 KB. The peak online RAM is 42 MB. The average runtime for a time step is 57 ms for the server and 12 ms for the client. 5) Online: We now examine the encrypted online control algorithm simulation adding details that complement the information in Section VII.
Memory: The peak RAM for both the client and the server incurred during the three phases for 45 time steps was 6.26 GB, out of which, 3.9 MB is the public key, 2 MB is the secret key, 15.7 MB is the relinearization key and 4.6 GB are the evaluation rotation keys (corresponding to a key for each of the 336 rotation indices), generated offline. The client discards the evaluation keys after it sends them to the server. At every iteration, the client receives two ciphertexts from the server (the Schur complement and control action) and sends three ciphertexts (inverse of the Schur complement, measurement, control action-we assumed here that the setpoint is constant). During the refresh time steps, one extra ciphertext is sent and received (the inverse matrix). The size of these ciphertexts depends on the level where the ciphertext is at. Specifically, for our selected parameters described in the beginning of the subsection, the ciphertexts are 2.77 MB at maximum size and 240 KB at the minimum size.
Time: The runtimes for the encrypted computation are given in Figure 4, where three different phases are depicted: trajectory concatenation, online sample collection and static update. The total offline initialization time for key generation and ciphertext encryption is 60 s for both the client and the server. In all online phases, the client only performs cheap and fast computations that take less that 0.3 s. The first online phase is the trajectory concatenation phase, as described in Section III-C, from t = M − 1 = 3 to t = M + N − 1 = 7, where the server computes the encrypted control action only with the precollected data. This is very efficient, despite computing on the maximum ciphertext sizes, and results in a total computation time for one time step of 2 s. Before t = M − 1, the client applies random inputs.
The second phase is the online collection of new input-output samples, which implies modifying the Hankel matrices and computing the inverse matrix via the Schur complement at every time step M −1 t from step 8 to step 22. This phase is split in as many parts as the refreshing time dictates. In the simulation depicted in Figure 4, this corresponds to three parts. At the established refresh times (t = 12, 17), the server packs the matrix M −1 t into one ciphertext, sends it to the client to re-encrypt it with the maximum number of moduli, and unpacks it back into column-wise ciphertexts. First, the increase in runtime between the subsequent time steps 13 and 12, respectively 18 and 17, is given by the fact the ciphertexts are returned to the maximum number of moduli. Second, the increase in computation time from the first refresh at time 12 and 13 (47.23 s and 104.6 s) to the second refresh at time 17 and 18 (63.9 s and 144.2 s) is given by the fact that the server has to deal with more collected samples ciphertexts than in the beginning of phase two. The intermediate decrease in the computation time is given by the decrease in the ciphertext size (we make sure to compute with only the minimal number of moduli required).
The third phase corresponds to the computations after stopping the collection of new samples, which starts from time step 23 and can go for the rest of the desired simulation time (accounting for Remark 2 in Appendix C). The third phase is only slightly more computationally intensive than the first. The time for the client halves compared to the previous two phases, because the client has to decrypt and encrypt ciphertexts with two moduli. The running time for the server substantially decreases compared to phase two but doubles compared to phase one, because we use a less efficient matrix-vector multiplications in order to minimize the multiplicative depth. Nevertheless, the running time required for computing the control input at one time step is around 3 s.

C. More security
In practice, the desired security parameter for a cryptographic application is at least 100 bits. The refresh time after 5 time steps (35 minutes) gives a ciphertext modulus of 802 bits and the refresh time after 8 time steps (56 minutes), which means a deeper circuit, gives a ciphertext modulus of 1120 bits. In Table I, we present simulation results where we choose the ring dimension 2 15 such that we satisfy this security requirement for the chosen ciphertext modulus. We simulate for the above 2x2 system and for a 4x4 system (where M = 4, N = 4, T = 64 and we again collect 15 more online samples, which gives a maximum number of 72 unknowns). We use 53 bits of precision for the plaintext moduli, such that we obtain 0.034 • C and 0.008 • C maximum, respectively average difference between the encrypted and unencrypted measurements, and 0.05 kWatts and 0.008 kWatts maximum, respectively average difference for the input.
We use different AWS cloud machines in order to satisfy the RAM and sampling time requirements. The c5.2xlarge machine (16 GB RAM) used 8 threads, the c5.4xlarge machine (32 GB RAM ) used 13 threads and the c5.9xlarge machine (72 GB RAM) used 28 threads.
From Table I, we see the trade-off between a longer refresh time (the server asks the client to refresh a ciphertext after more time, which implies a deeper circuit and less security for the same ring dimension) and the memory consumption and computation time. The conclusion is that the refresh time is a tuning knob between security level, communication and computational complexity which should be determined according to the application's specifications.
Note that as the performance of the underlying homomorphic encryption library improves, these improvements will reflect in our computation times as well. XI. CONCLUSIONS AND FUTURE WORK We proposed an online data-driven control algorithm designed to be encryption-friendly, i.e., allow both good setpoint tracking and efficient encrypted computations for privacy-preserving control. We achieve this by choosing a regularized relaxed convex formulation of the control optimization problem, such that we obtain a closed form solution for the control input, which we prove to be close to the solution of the noiseless original problem, under conditions on the regularization parameters. Moreover, this solution can be computed using a low-depth arithmetic circuit using the Schur's complement formula for matrix inversion and manipulating the sequential matrix-vector operations into different forms. These computations are amenable for an efficient encrypted solution using leveled homomorphic encryption and carefully encoding vectors into ciphertexts. The client encrypts its measurements with such an encryption scheme and sends them to a cloud server, which returns the encrypted control input. We discussed how to achieve an optimized implementation of the encrypted computations and how to make sure that the computations can run for a desired number of time steps. Finally, we discussed extensions that can be implemented to enhance the functionality of this encrypted control algorithm. These tools can be used to tackle other related problems, such as adaptive control or that involve sequential matrix updates.
The LQR data-driven approximation algorithm we presented in this paper is the first step towards private datadriven control. We plan to include constraints on the inputs and the outputs. This will substantially increase the complexity of the encrypted computations, since comparisons are not inherently polynomial operations. We plan to explore other encryption-friendly methods for private data-driven control, that are less data-hungry than the behavioral framework and more amenable for larger systems.

APPENDIX A PROOF OF CLOSENESS
In this section, we present a detailed proof of Theorem 1. First, we characterize the solution set of the original behavioral problem (5). Then, we characterize the solution of the approximate problem (6). To prove closeness, we make suitable use of the inversion lemma, the pseudo-inverse (denoted by †) limit definition, and Singular Value Decomposition (SVD).

A. Notation and preliminaries
Before we proceed with the proof, we consider the following definitions: whereM is matrix M defined in (8), without the penalty terms. Similarly,M 0 is matrix M without any regularization or penalty terms. For compactness, we denote the matrix of past precollected data as D p and the vector of the last M outputs and inputs as d t : Both matricesM 0 , D p are singular. Let the SVD of the past precollected data D p be: where Σ ∈ R q×q is an invertible diagonal matrix, containing the non-zero singular values of D p with q = rank(D p ), while matrices E E ⊥ , F F ⊥ are orthonormal. Recall that the set of optimal solutions of the original behavioral problem (5) is denoted by G opt -see (11). The minimum norm element of G opt is denoted by g minsee (12). For completeness we present some inversion formulae.
Lemma 2 (Inversion formulae). Let K, L, V be matrices of appropriate dimensions. Then: a) Provided the respective inverses exist: b) Let V V ⊥ be an orthonormal matrix. Then, if L is invertible: Proof. Part a) is standard, also known as the Woodbury matrix identity. To prove (34), we use the identity V ⊥ V ⊥ + VV = I and verify that the properties of the inverse hold. To prove (35), we use (34).

B. Original behavioral problem
In the following lemma, we characterize the solution set G opt of the original constrained behavioral problem (5). In general, there are infinite optimal solutions sinceM 0 and D p have low rank.
Lemma 3. Consider optimization problem (5). Recall the SVD decomposition of D p in (33). Then a) g ∈ G opt is an optimal solution of (5) if and only if g lies in the following affine subspace: for some µ ∈ R q , where q is the rank of the past precollected data D p . b) Vector Y f Qr t lies in the column space ofM 0 . c) The above affine subspace can be equivalently described by the closed form expression: where g i , i = 1, 2, 3 are orthogonal with each other and where s is arbitrary. d) As a result, the minimum norm element of G opt is Proof. Proof of a) Since d t lies in the range space of D p , we can replace the singular equality constraint where now F has full rank. If we also replace u, y with U f g and Y f respectively, and we ignore the constant terms, then we obtain the optimization problem: which is equivalent to (5). The Lagrangian of (42) is where µ ∈ R q are the Lagrange multipliers. The result follows by applying standard first order optimality conditions.

Proof of b) Observe thatM
It is clear that Y f Qr t lies in the column space of Y f Q 1/2 U f R 1/2 . Then, using the fact that for a matrix A, Range(A) = Range(AA ) (the proof easily follows from using SVD on A), we also show that Y f Qr t lies in the column space ofM 0 . Proof of c) First, we solve subequation: Since matrix F F ⊥ is an orthonormal basis, we can express all solutions g as: where F ⊥ ξ captures the components of g in the kernel of F , and ξ is to be determined. Second, we solve subequation:M Multiplying from the left by the orthonormal matrix F F ⊥ results in the equivalent system of subequations: The former subequation (44) just determines the Lagrange multiplier µ. Replacing (43) into (45) gives: By the standard property of the pseudo-inverse [42, Ch. 2], we can write all possible solutions ξ as: where s is arbitrary. Moreover, the two components ξ 1 , ξ 2 are orthogonal. Combining (43) and (46) gives: where all components are orthogonal to each other. Proof of d) Follows from b), orthogonality of the three summands and setting the only free parameter s to zero.

C. Approximate problem
Due to the regularization term, the solution of the approximate problem (6) is unique-see (9). We can show that if we chose large enough λ, then we can show that the solution of (6) converges to an intermediate one.
Lemma 4. Consider optimization problem (6) with λ u = λ y = λ > 0 and λ g > 0 and let g * be the optimal solution. Define the intermediate solution: Then, as λ → ∞ and λ g → 0 + , where we used big-O notation to hide quantities which do not depend on λ, λ g . Proof.
Step a). Recall that g * is given by (9). We will use the inversion lemma to bring g * in a form closer tō g. Using the notation of this section, (9) can be written as: where the second equality follows from D p D p = FΣ 2 F and the last equality follows from F F = I. By Lemma 2 a): Step b). Next, define the difference: where the equality follows from The error g * −ḡ can be written as: To complete the proof, it is sufficient to show that the matricesM Step c) Boundedness ofM −1 F F M −1 F −1 follows by Lemma 5.
Step d) The following matrix is always bounded: It follows fromM (λ g + M 0 2 )I, which in turn implies where , denote comparison with respect to the positive semi-definite cone.

D. Proof of Theorem 1
Consider the intermediate solutionḡ defined in (47). By the triangle inequality we have: In Lemma 4, we show that the second error term decays to zero as fast as λ −1 . To complete the proof, we invoke the following Lemma, which states that the first error term decays to zero as well.
Lemma 5. Consider the minimum norm solution g min defined in (12) and the intermediate solutionḡ defined in (47). Let λ g > 0. Then: a) The following limit holds b) The following matrix is bounded c) The intermediate solution converges to the minimum norm solution as λ g goes to zero Proof. Proof of a). Let T :=M 1/2 0 F ⊥ . Since F ⊥ F ⊥ = I: The result follows from the pseudo-inverse limit definition and T † = (T T) † T [42, Ch 1].
Proof of b). Note that the norm remains unchanged if we multiply from the right by F .

By (35) and the triangle inequality
where the last equality follows from the fact that multiplying from the right by F F ⊥ leaves the norm unchanged.
For the other submatrix: (49), we immediately obtain boundedness since: Hence, both submatrices have bounded norm, which implies boundedness for the original matrix. Proof of c). Using (35), we can rewriteḡ in a form that resembles g min : By (52), we get that the second term converges to For the third term, recall from Lemma 3 b), that Y f Qr t lies in the column space ofM 0 , which implies that Y f Qr t =M 0M † 0 Y f Qr t . Using this identity and (49), we also get that the third term converges to The result follows by (41), (53) and the convergence of the aforementioned terms.

APPENDIX B CKKS SCHEME
We present a schematic view of the CKKS leveled homomorphic cryptosystem, introduced in [9]. The plaintext space is the ring of polynomials R := Z[X]/ Φ M (X) , where Φ M (X) is the M -th cyclotomic polynomial of degree N := φ(M )/2, i.e., X N + 1. N is also called the ring dimension for the plaintext space. The ciphertext space for level l is the ring (Z Ql [X]/ Φ M (X) ) 2 , where Q l is the ciphertext modulus corresponding to level l. Define also the residue ring R l := R/Q l R.
The CKKS scheme supports an efficient packing of up to N/2 real values in a single plaintext, respectively a single ciphertext. A real-valued message in C N is encoded in a plaintext in Z[X]/ Φ M (X) via the inverse canonical embedding map and rounding-which corresponds to evaluation of the inverse Discrete Fourier Transform (with the primitive M -th roots of unity exp(−2πi/M )). The decoding procedure is given by the forward Discrete Fourier Transform. This batching technique allows performing Single Instruction Multiple Data (SIMD) operations on the ciphertexts.
Each operation introduces some noise in the ciphertext, see [9] for the magnitude of the different noises. The multiplication operation (along with rescaling) introduces the largest noise. Hence, the multiplicative budget is set from the beginning as a parameter of the scheme, and correct decryption is only guaranteed for the evaluation of circuits of smaller multiplicative depth. To deal with precision, in the CKKS scheme, a scaling factor is multiplied to each value. After a multiplication (including multiplication by plaintexts), one scaling factor needs to be removed in order to keep the correctness. This is done through a rescaling procedure which decreases the ciphertext modulus.
The semantic security of the CKKS scheme assumes that the Decisional Ring Learning with Errors problem is computationally hard [46].

A. Ciphertext packing
Assume without loss of generality that we shift the time axis to the left by M + N , such that the trajectory concatenation mentioned in Section III-C is performed before t = 0. This will simplify the circuit depth expression. In the following, we will address some of the encoding methods we investigated, with the goals to minimize the depth of the resulting circuit and to reduce the memory consumption and time complexity.
1) Scalar encoding: each ciphertext only encodes and encrypts a scalar, i.e., element of a vector or matrix. This version has the smallest multiplicative depth compared to the other packing methods, no rotation key storage, but the largest ciphertext storage and the largest number of operations. It returns m ciphertexts instead of one to the client, unless we add more processing to mask and rotate the elements of u t to obtain a single ciphertext.
2) "Hybrid" encoding: some ciphertexts encode and encrypt individual elements while other ciphertexts encode and encrypt vectors. Specifically, we pack the input and output vectors and columns of input and output matrices, but hold the entries of the matrix M −1 t and intermediate vectors m t as individual ciphertexts. This version returns only one vector for the client to decrypt and involves the circuit depth given in (54).
The "hybrid" encoding was implemented in [22]. The conclusions there were that for the best compromise of depth versus client computation and communication, we ask the client to send back an encryption of 1/s t , along with an encryption of the input u t and output y t . This gives the depth expression of the denominator and numerator of u t , for t > 1): However, this "hybrid" solution still incurred a large storage: namely O((S + t) 2 ) ciphertexts and O((S + t) 2 ) rotation keys (only O(S + t) rotation keys when there is no refreshing and packing of the matrix M −1 t required), and O((S + t) 2 ) computations at the cloud, with large hidden constants. This suggests the need for a more tractable solution, involving a more efficient encoding.
3) Vector encoding: each ciphertext encodes and encrypts multiple values, i.e. a vector or a column matrix. Via this type of encoding, we aim to keep the minimum possible depth as in the previous two versions, but minimize the memory requirements for ciphertexts, from O((S + t) 2 ) to O(S + t). At the same time, we want to decrease the rotation key storage from O((S + t) 2 ) to O(S + t) and keep the number of computations at number of computations at O((S + t) 2 ).
With the previous "hybrid" encoding on the inputs in Appendix C-A2, used in [22], we found it preferable to encode the entries of M −1 t in individual ciphertexts in order to avoid masking when computing the desired result, which would have increased the depth. However, we found that using redundancy in the stored variables, i.e., repeating the elements in vector encoding of the input columns, along with completely rewriting the way the operations are performed, can alleviate this issue and keep the same depth as before.
We will use the encryption of the repeated elements encoding, E vr0 (·), for the columns of the Hankel matrices HU t and HY t , the vectorsū t ,ȳ t , the inputs and outputs u t , y t , the reference r t , and to encode the costs. In the encoding, each element is repeated for S +T times, whereT is the maximum number of online collected samples.
The precollected input and output measurements have a multiplicative depth of 0. To reduce the total circuit depth, at time t + 1, although the cloud service has computed the encryption of u t , the client sends a fresh encryption of it. This allows us to have, ∀t ≥ 0: d(colsHU t ) = d(colsHY t ) = d(ū t ) = d(ȳ t ) = d(r t ) = 0.
The update of the row representation of U f needs a masking, so d(rowsU f ) = 1. From the offline computations, we assume that we have fresh encryptions for the columns of M −1 0 , hence d(colsM −1 0 ) = 0. The encrypted operations, their flow and their motivations are described in Section VII-B. Compared to the previous approach in [22], we observed that we can avoid extra computations and a possible extra level in the computations if the server asks the client to compute 1/s t immediately after it computes s t . Hence, the server can compute directly the encryptions of M −1 t and u t instead of the numerators N M −1 t and N u t . The encrypted operations are described in (55)-(57), which correspond to lines 6 and 7 from Algorithm 1. EvalSum S+T refers to rotating and summing batches of S +T slots in the ciphertexts and v t := (m t M −1 t−1 ) . E vr * (µ t ) = EvalSum S+T (E vr0 (col S +1 HY t−1 ) e vr0 (λQ) E vr0 (col S +1 HY t−1 )+ + E vr0 (col S +1 HU t−1 ) e vr0 (λR) E vr0 (col S +1 HU t−1 ) + e vr0 (λ g ) E vr * (m t[i] ) = EvalSum S+T (E vr0 (col S +1 HY t−1 ) e vr0 (λQ) E vr0 (col i HY t−1 )+ + E vr0 (col S +1 HU t−1 ) e vr0 (λR) E vr0 (col i HU t−1 ), i ∈ [S + 1] E vr * (s t ) = E vr * (µ t ) − S i=0 E vr * (m t[i] ) ρ (E v0 (v t ), i) .
After receiving a fresh encryption of 1/s t from the client: (1/s t )), S + 1 , i ∈ [S + 2] E vr * (Z t[j] ) = EvalSum S+T (E vr0 (col j HY t ) e vr0 (λQ) E vr0 (ȳ t , r t ) + E vr0 (col j HU t ) e vr0 (λR) E vr0 (ū t )) After decryption: Following (55)-(57), we get the depth expressions for t ≥ 1: Remark 2. We make an abuse of notation when writing that we obtain E v0 (·) and E vr0 (·) after the rotation of a ciphertext with trailing zeros. Actually, the initial non-zero values will be shifted to the end of the encoded vector and we need to maintain a counter on how many times we can perform these rotations. In the vector packing case, we need the parameters to satisfy the following rules, before a masking is necessary: where the first line is for the online collection of samples (T is the total number of samples collected online), the second line is for the correct packing of the matrix into one ciphertext and the last line is the rule for after finalizing the collection of new samples t ≥T (instead of masking, here we can ask the client to prune the junk elements inū,ȳ).
If the rules in Remark 2 are not satisfied for the desired parameters of an application, the values can be split into the necessary number of ciphertexts (each packing up to ringDim/2 values) and the computations can be readily extended to deal with this.
(d) The proof is laborious so, for conciseness, we do not present the complete derivations. We show that f (λ g ) is decreasing by showing that ∂f (λg) ∂λg := − g(λg) λ 2 g takes only negative values, which we show by proving lim λg→0 + g(λ g ) ≥ 0 and ∂g(λg) ∂λg ≥ 0. To this end, we used again the pseudoinverse limit definitions and that (P 1/2 H)(P 1/2 H) † is an orthogonal projection onto the range of (P 1/2 H) and hence, its eigenvalues are in {0, 1}.