# Homomorphically encrypted gradient descent algorithms
for quadratic programming

André Bertolace, Konstantinos Gatsis, Kostas Margellos The authors are with
the Department of Engineering Science, University of Oxford, Oxford OX1 3PJ,
U.K. E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>(corresponding author) ©2023 IEEE. Personal use
of this material is permitted. Permission from IEEE must be obtained for all
other uses, in any current or future media, including reprinting/republishing
this material for advertising or promotional purposes, creating new collective
works, for resale or redistribution to servers or lists, or reuse of any
copyrighted component of this work in other works.

###### Abstract

In this paper, we evaluate the different fully homomorphic encryption schemes,
propose an implementation, and numerically analyze the applicability of
gradient descent algorithms to solve quadratic programming in a homomorphic
encryption setup. The limit on the multiplication depth of homomorphic
encryption circuits is a major challenge for iterative procedures such as
gradient descent algorithms. Our analysis not only quantifies these
limitations on prototype examples, thus serving as a benchmark for future
investigations, but also highlights additional trade-offs like the ones
pertaining the choice of gradient descent or accelerated gradient descent
methods, opening the road for the use of homomorphic encryption techniques in
iterative procedures widely used in optimization based control. In addition,
we argue that, among the available homomorphic encryption schemes, the one
adopted in this work, namely CKKS, is the only suitable scheme for
implementing gradient descent algorithms. The choice of the appropriate step
size is crucial to the convergence of the procedure. The paper shows firsthand
the feasibility of homomorphically encrypted gradient descent algorithms.

##  I Introduction

Homomorphic encryption (HE) is a ground-breaking mathematical method that
enables the analysis or manipulation of encrypted data without revealing its
content [14]. In doing so, HE permits the secure delegation of data processing
to third-party cloud providers. Several encryption schemes, such as Paillier
[26] or El Gamal [9] are partially HE schemes111Partially HE schemes enable
the implementation of either addition or multiplication on encrypted data, but
not both, whereas fully HE schemes enable the implementation of both addition
and multiplication operations.. In 2009, Gentry [13] proposed the first fully
HE scheme${}^{\text{1}}$. The computational overhead of the scheme was
significant, but it showed that such schemes are indeed possible. Since then
these approaches have been further developed. Currently, the state of the art
schemes are BFV [3, 10], CKKS [7], BGV [4], and GSW [15]. The computational
overhead remains large but it has been brought down to a level where these
schemes can be implemented in practice.

In applications of control and decision-making, the benefits of delegating
data processing without giving away access to the data are tremendous. The use
of HE schemes applied in control theory is at its infancy, however, encrypted
linear controllers have been implemented. Most results use partially HE
schemes [1, 2, 28, 6, 11, 12], but approaches using fully HE schemes also
exist [21, 29, 23]. In addition, [27] provides a detailed overview of the
current status of research in the encrypted control for networked systems and
[20] a comparison of different encrypted control approaches.

Yet, the implementation of algorithms in a HE setup is far from trivial. For
instance, many HE schemes use random noise to guarantee the security of the
encrypted data. This noise compounds at every arithmetic operation, resulting
in a limited number of sequential arithmetic operations performed by an
encryption circuit, in special, multiplication operations.

In this paper, we would like to understand the limits imposed by HE
computation on challenging computation tasks beyond the controller
implementation problems studied in the literature. Specifically, we consider
the problem of solving quadratic programming (QP) problems. QP is commonly
used in several control problems like those arising in state estimation under
minimum square error and model predictive control. Numerically solving such a
task often requires iterative methods (gradient descent) and the limit on the
multiplication depth of HE circuits is a major challenge for iterative
procedures. As a result, given the HE multiplication depth limits, we would
like to determine the most appropriate iterative methods for QP. In our case,
we adapt and implement gradient descent (GD) and accelerated gradient descent
(AGD) algorithms to solve a QP in a HE manner. Our contributions are
threefold:

  1. 1.

We argue that among the available HE schemes CKKS is the only scheme suitable
to handle GD and AGD iterations in a HE setup as it allows handling real-
valued operations, an important feature, especially in the selection of an
appropriate step-size that ensures the convergence of the underlying
algorithm.

  2. 2.

We implement our own HE matrix multiplication algorithm that is more
efficient, in terms of multiplication depth, than other algorithms proposed in
the literature [18].

  3. 3.

We demonstrate that in the HE setup the condition number of the matrix of the
quadratic term of the objective function plays an important role in
determining which algorithm is preferred. The encrypted version of AGD is the
preferred algorithm only for matrices with higher condition numbers. This is
in contrast to plain-text optimization, where AGD is typically preferable due
to its superior convergence rate. The reason being once an encryption circuit
is defined, the number of sequential arithmetic operations is fixed but at
each step, AGD performs one extra multiplication when compared to GD. These
extra iterations pays off, specially for matrices with lower conditional
number.

Other works in the literature, [1], also proposed to solve QPs in a
secure/distributed manner. We differ from this work by using fully homomorphic
encryption instead of partial homomorphic encryption schemes. Furthermore,
[29], presented an encrypted model predictive control scheme for linear
constrained systems, also using partial homomorphic encryption schemes and not
solving the QP but instead using the corresponding piece-wise affine control
law if explicitly given.

It should be mentioned that our paper focuses on unconstrained quadratic
problems due to inherent limitations in operations that are allowed to be
performed by available HE schemes. However, our analysis not only quantifies
these limitations on prototype examples thus serving as a benchmark for future
investigations, but also highlights additional trade-offs like the ones
pertaining the choice of GD or AGD methods, opening the road for the use of HE
methods in iterative methods widely used in optimization based control.

Finally, to emphasize that this line of research is still at an early stage,
we note that our implementation using state of the art HE tools (Microsoft
SEAL) permits just a modest amount of gradient steps, beyond what would be
required in practical applications.

The rest of the paper is organized as follows: Section II provides some
background information on gradient methods, while Section III discusses our
choice for an HE scheme. Section IV discusses the implementation of the
suggested scheme and provides an extension to the matrix multiplication
operation. Section V provides a detailed numerical study on the application of
GD and AGD algorithms for QPs, while Section VI concludes the paper.

##  II Descent algorithms for unconstrained quadratic programming

QP has been a very successful tool for modeling many real-life problems. It is
extensively used in applications that involve the variance minimization, such
as in the formulation of portfolio optimization problems or in solving the
ordinary least square (OLS) problem. In fact, many problems in physics on
engineering can be formulated as some form of energy minimization problem, in
which the energy can simply be formulated as a quadratic form,

$\underset{{x\in\mathbb{R}^{n}}}{\mbox{{min}}}\;\Big{\\{}{\frac{1}{2}x^{T}Qx+p^{T}x}\
\Big{|}\ {}\Big{\\}},$

where $Q\in\mathbb{R}^{n\times n}$ and $p\in\mathbb{R}^{n}$. Note that the QP
is convex if $Q\succcurlyeq 0$.

This unconstrained QP has a closed form solution, it requires however the
inversion of a matrix, a procedure that involves other operations than
additions and multiplications, posing hence a challenge for its implementation
in a HE setup. An alternative solution is to use gradient descent methods to
solve the QP problem. This is a class of iterative algorithms that provide a
simple way [5] to minimize a differentiable function $f$,

$\underset{{x\in\mathbb{R}^{n}}}{\mbox{{min}}}\;\Big{\\{}{f(x)}\ \Big{|}\
{}\Big{\\}}.$

Starting at an initial estimate, it iteratively updates

$x_{t+1}=x_{t}-\eta\nabla f(x_{t}),$

where $\nabla f(x_{t})$ denotes the gradient of $f$ calculated at $x_{t}$ and
$\eta$ the step-size, until reaching a desired tolerance in the solution.
Particularly to the QP case, the gradient takes a linear form, $\nabla
f(x)=Qx+p$.

Methods of this type have a convergence rate which is independent of the
dimension $n$ of the solution space. This feature makes them particularly
attractive for optimization in very high dimensions [5]. The convergence is
however deeply linked to the step-size $\eta$, be it too small, the algorithm
may take too long to converge, be it too high, it may diverge.

Properties such as smoothness or strong convexity of the objective function
$f$ do play a relevant role in choosing $\eta$ and a variant of the algorithm
with faster convergence.

###### Definition II.1

A continuous differentiable function $f$ is $\beta$-smooth if the gradient
$\nabla f$ is $\beta$-Lipschitz, i.e.,

$\left\lVert\nabla f(x)-\nabla f(y)\right\rVert\leq\beta\left\lVert
x-y\right\rVert,~{}\forall x,y\in\mathbb{R}^{n}.$

###### Definition II.2

A function $f$ is $\alpha$-strongly convex, with $\alpha>0$, if for any $x$,
$y$ it satisfies the following sub-gradient inequality, i.e.,

$f(x)-f(y)\leq\nabla f(x)^{T}(x-y)-\frac{\alpha}{2}\left\lVert
x-y\right\rVert^{2},~{}\forall x,y\in\mathbb{R}^{n}.$

Given these definitions, an immediate consequence is that if $f$ is twice
differentiable, then $f$ is $\alpha$-strongly convex if the eigenvalues of the
Hessian of $f$ are larger than or equal to $\alpha$. A quadratic $f$, as in
our QP, is $\lambda_{\max}$-smooth and $\lambda_{\min}$-strongly convex, where
$\lambda_{\max},\lambda_{\min}>0$ are, respectively, the maximum and minimum
eigenvalues of the matrix $Q$. The ratio
$\kappa=\frac{\lambda_{\max}}{\lambda_{\min}}$ is the condition number of the
matrix $Q$ and it plays an important role in the convergence of descent
algorithms. For reference, both Nesterov’s accelerated gradient descent (AGD)
and Gradient Descent (GD) methods [5] for the smooth and strongly convex
quadratic function converge exponentially fast (see Table I):

$f$ |  Algorithm |  Convergence rate |  Iterations needed
---|---|---|---
$\alpha$-conv., $\beta$-smooth |  AGD |  $R^{2}\exp\biggl{(}-\frac{t}{\sqrt{\kappa}}\biggr{)}$ |  $\sqrt{\kappa}\log\biggl{(}\frac{R^{2}}{\epsilon}\biggr{)}$
$\alpha$-conv., $\beta$-smooth |  GD |  $R^{2}\exp\biggl{(}-\frac{t}{\kappa}\biggr{)}$ |  $\kappa\log\biggl{(}\frac{R^{2}}{\epsilon}\biggr{)}$
TABLE I: Convergence rate and other parameters for different algorithms.
$R=\left\lVert x_{1}-x^{*}\right\rVert$ is the distance from the initial
estimate to the optimal value. The number of iterations is directly derived
from the convergence rate for a fixed tolerance (in terms of distance from
optimal value) $\epsilon$.

##  III Homomorphically encrypted arithmetic

###  III-A Homomorphic encryption schemes

HE schemes have been developed using different approaches. BFV and BGV perform
operations modulo integer whereas CKKS implements approximated fixed point
arithmetics. The security of these schemes is based on the Ring Learning With
Errors (RLWE) problem, a variant of the Learning With Errors problem (LWE), in
which the goal is to distinguish random linear equations, which have been
perturbed by a small amount of noise from uniform ones [25]. The HE cipher is
defined by a pair $E,D$, of encryption-decryption algorithms respectively. $E$
takes a public key $pk$ along with a message $m$ as inputs and outputs a
cipher-text $\mathtt{c}$, as $\mathtt{c}=E(pk,m)$. The decryption algorithm,
$D$, takes a secret key $sk$ along with the cipher-text $\mathtt{c}$ as inputs
and outputs the message $m=D(sk,\mathtt{c})$. The algorithms are parameterized
by a security parameter $\lambda$ which plays a direct role in the derivation
of the $sk$. In addition, these schemes exploit the structure of polynomial
rings for its plain-text and cipher-text spaces, the cyclonomic polynomial,
$R[\mathbb{Z}_{q}]=\mathbb{Z}_{q}[X]/(X^{N}+1)$. All schemes make use of
random variables with values sampled from a discrete Gaussian distribution
with a pre-defined variance and random variables sampled from a ternary
distribution $\\{-1,0,1\\}$ [7].

BFV: In the BFV scheme, [3, 10], the plain-text and cipher-text spaces are
defined by two distinct rings, $R[\mathbb{Z}_{t}]$ and $R[\mathbb{Z}_{q}]$,
where $t$ and $q$ are parameters of the plain-text and cipher-text
coefficients, respectively.

BGV: The BGV scheme, [4], is similar to the BFV. The plain-text and cipher-
text spaces are defined by two distinct rings, $R[\mathbb{Z}_{t}]$ and
$R[\mathbb{Z}_{q}]$.

CKKS: The CKKS scheme, [7], is often quoted as being the most efficient method
to perform approximate HE computations over real and complex numbers [19]. It
can be considered as a noisy channel [24], a simple encryption/decryption
procedure adds noise to the original message. The scheme exploits the
structure of integer polynomial rings for its plain-text and cipher-text
spaces, $R[\mathbb{Z}_{q}]$ and $R[\mathbb{C}]$. This polynomial ring is
combined with a canonical embedding transformation
$\sigma:\mathcal{S}\rightarrow\mathbb{C}^{N}$ that encodes/decodes a
vector222The space size is actually $N/2$ because the roots of the cyclonomic
polynomial lie on the unit circle and are pairwise complex conjugate. in
$\mathbb{C}^{N}$ to/from the ring of cyclonomical complex polynomials
$R[\mathbb{C}]$. To encode a message $x\in R[\mathbb{C}]$ one applies the
inverse embedding transformation to get $\mu=\sigma^{-1}(x)\in R[\mathbb{C}]$,
then scale $\mu$ by a factor $\Delta=2^{p}$ and round to obtain the plain-text
$m=\lfloor\Delta\cdot\mu\rceil\in R[\mathbb{Z}_{q}]$.

###  III-B Scheme choice

We claim that the most suitable choice for HE versions of the GD algorithm
(similar considerations hold for the AGD one) is the CKKS scheme. The main
reason for such claim is related to the selection of the step-size $\eta$.
Note that $\eta$ should be sufficiently small for GD to converge. This is
summarized in the following proposition; it is a standard result but we
present a proof below for completeness.

###### Proposition III.1

Consider a QP with $Q\succcurlyeq 0$ with $Q=Q^{\top}$, and let
$\lambda_{\max}$ denote the maximum eigenvalue of $Q$. The GD method converges
for any $\eta<\frac{2}{\lambda_{\max}}$.

###### Proof:

Given $f(x)=\frac{1}{2}x^{T}Qx+p^{T}x$, we have that $\nabla f(x)=Qx+p$ and
the iterative GD procedure takes the form $x_{t+1}=(I-\eta Q)x_{t}-\eta p$.
Let $x^{*}$ be an unconstrained minimizer of $f$. As such, $\nabla
f(x^{*})=Qx^{*}+p=0$, which in turn implies that $x^{*}=(I-\eta Q)x^{*}-\eta
p$. We thus have that $x_{t+1}-x^{*}=(I-\eta Q)(x_{t}-x^{*})$ and consequently
$\left\lVert x_{t+1}-x^{*}\right\rVert\leq\left\lVert I-\eta
Q\right\rVert\left\lVert(x_{t}-x^{*})\right\rVert\leq\left\lVert I-\eta
Q\right\rVert^{t+1}\left\lVert(x_{0}-x^{*})\right\rVert$. The latter implies
that $\lim_{t\to\infty}\left\lVert x_{t+1}-x^{*}\right\rVert=0$ if the maximum
eigenvalue of $(I-\eta Q)$ is less than $1$, which can be achieved if
$\eta<\frac{2}{\lambda_{\max}}$. ∎

A direct consequence of this fact is that if using BGV or BFV that require
integer step-sizes, one can only ensure convergence for matrices with
$\lambda_{\max}<2$ that is the only choice that allows for an integer step-
size $\eta$. The minimum then value of such step-size would be $\eta=1$, which
in turn may lead to an erratic numerical behaviour. Additionally, to be able
to use BFV or BGV, one would need to limit the calculations to integer
matrices $Q\in\mathbb{Z}^{n\times n}$, or manipulate $Q\in\mathbb{R}^{n\times
n}$ to be made integer. Towards this direction, [20, 22], suggest the
following manipulations:

  * •

”Scaling-up” the real numbers by a factor, say $10^{8}$, replicating a fixed
point arithmetic, and proceed by calculating using the given integer numbers.

  * •

Converting the matrix $Q$ by finding an invertible matrix
$T\in\mathbb{R}^{n\times n}$ such that $TQT^{-1}\in\mathbb{Z}^{n\times n}$.

The former is not a practical solution as the result of multiple
multiplications will overflow and the output after the decryption will be
incorrect [20], whereas the latter implies limiting ourselves to matrices $Q$
in which every eigenvalue has an integer real and imaginary part [20, 22]. In
summary, BGV and BFV are only suitable schemes for integer matrices or
matrices that have integer eigenvalues, which for our setting would require
$\lambda_{\max}<2$. This would imply working only with identity matrices,
$Q=I$, if we working with integer matrices which are symmetric and positive
definite. As such, for the purpose of an iterative methodology like GD and
AGD, CKKS is preferable.

##  IV Homomorphically encrypted gradient descent algorithms

###  IV-A Algorithm description

For the HE version of gradient descent methods, let us start by defining the
following arithmetic operators:

  * •

$\mathord{\leavevmode\hbox to8.9pt{\vbox
to8.9pt{\pgfpicture\makeatletter\hbox{\hskip 0.75pt\lower-0.75pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\pgfsys@setlinewidth{1.5pt}\pgfsys@invoke{ } {}{{}}{}
{}{}{}\pgfsys@moveto{3.69899pt}{0.0pt}\pgfsys@lineto{3.69899pt}{7.398pt}\pgfsys@stroke\pgfsys@invoke{
} {}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{3.69899pt}\pgfsys@lineto{7.398pt}{3.69899pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}$/$\mathord{\leavevmode\hbox
to7.19pt{\vbox to1.5pt{\pgfpicture\makeatletter\hbox{\hskip 0.75pt\lower
1.52637pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\pgfsys@setlinewidth{1.5pt}\pgfsys@invoke{ } {}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{2.27637pt}\pgfsys@lineto{5.69089pt}{2.27637pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}$: the
addition/subtraction of two cipher-texts;

  * •

$\bullet$: the multiplication of two cipher-texts;

  * •

$\odot$: the multiplication of a plain-text and a cipher-text;

Let us further assume that the user calculates $\lambda_{\min}$,
$\lambda_{\max}$, and sends these as plain-text, i.e. not encrypted, values to
the solver. Together with these constants, the user also sends the encrypted
matrix and vector, $\mathtt{Q}=E(pk,Q)$ and $\mathtt{p}=E(pk,p)$ that
determine the QP. The encrypted version of the descent algorithms will still
proceed in an iterative fashion. The only difference is that one would be
iterating over cipher-texts $\mathtt{c_{t}}$ instead of plain-text $x_{t}$.
When iterating over cipher-texts, two steps deserve special attention, the
stopping criteria $\left|\mathtt{c_{t+1}}\mathord{\leavevmode\hbox
to7.19pt{\vbox to1.5pt{\pgfpicture\makeatletter\hbox{\hskip 0.75pt\lower
1.52637pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\pgfsys@setlinewidth{1.5pt}\pgfsys@invoke{ } {}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{2.27637pt}\pgfsys@lineto{5.69089pt}{2.27637pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\mathtt{c_{t}}\right|>\epsilon$
and the matrix multiplication procedure, referred to MMult (Algorithm 4) and
discussed in the sequel. The latter is relevant because of the exponential
growth of the noise level with the multiplication depth. The authors in [18]
do propose a fully-homomorphic matrix multiplication algorithm to perform
these operations in a more efficient manner. On the other hand, for the
stopping rule, determining whether an encrypted value is larger than another
encrypted value or even a plain-text without decrypting both values is
directly not feasible, but complex approaches to implement comparisons have
appeared in [17] for BFV and [8] for CKKS.

Given the challenge to implement the stopping rule in an HE setup, we propose
that the HE version of the AGD algorithm is slightly modified:

  * •

Instead of specifying the tolerance $\epsilon$, the user fixes the number of
iterations $N$. Given $N$ and the convergence rate of the algorithms we can
infer the attainable tolerance (Table I).

  * •

The user may hand in the initial estimate $\mathtt{x_{0}}$, although this is
not necessary.

The HE versions of the the AGD and GD still follow an iterative procedure.
These take the form of (Algorithm 1) and (Algorithm 2) respectively and are
very similar to the usual AGD and GD algorithms. The main difference is the
use of HE arithmetic operators and the special MMult matrix multiplication
procedure.

Algorithm 1 HE AGD for an unconstrained QP

1:function
HEagdQP($\mathtt{Q},\mathtt{p},d,\lambda_{\min},\lambda_{\max},\mathtt{x_{0}},N$)

2: $\kappa\leftarrow\frac{\lambda_{\max}}{\lambda_{\min}}$

3: $\mathtt{x_{-}}\leftarrow\mathtt{x_{0}}$

4: $\mathtt{y_{-}}\leftarrow\mathtt{x_{0}}$

5: $\eta\leftarrow\frac{-1}{\lambda_{\max}}$

6: for $t=0$ to $N-1$ do

7: $\mathtt{y_{+}}\leftarrow\mathtt{x_{-}}\mathord{\leavevmode\hbox
to8.9pt{\vbox to8.9pt{\pgfpicture\makeatletter\hbox{\hskip
0.75pt\lower-0.75pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\pgfsys@setlinewidth{1.5pt}\pgfsys@invoke{ } {}{{}}{}
{}{}{}\pgfsys@moveto{3.69899pt}{0.0pt}\pgfsys@lineto{3.69899pt}{7.398pt}\pgfsys@stroke\pgfsys@invoke{
} {}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{3.69899pt}\pgfsys@lineto{7.398pt}{3.69899pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\textsc{MMult}(\mathtt{Q},\mathtt{x_{-}},d,\eta)\mathord{\leavevmode\hbox
to8.9pt{\vbox to8.9pt{\pgfpicture\makeatletter\hbox{\hskip
0.75pt\lower-0.75pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\pgfsys@setlinewidth{1.5pt}\pgfsys@invoke{ } {}{{}}{}
{}{}{}\pgfsys@moveto{3.69899pt}{0.0pt}\pgfsys@lineto{3.69899pt}{7.398pt}\pgfsys@stroke\pgfsys@invoke{
} {}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{3.69899pt}\pgfsys@lineto{7.398pt}{3.69899pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\eta\odot\mathtt{p}$

8:
$\mathtt{x_{+}}\leftarrow\biggl{(}1+\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\biggr{)}\odot\mathtt{y_{+}}\mathord{\leavevmode\hbox
to7.19pt{\vbox to1.5pt{\pgfpicture\makeatletter\hbox{\hskip 0.75pt\lower
1.52637pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\pgfsys@setlinewidth{1.5pt}\pgfsys@invoke{ } {}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{2.27637pt}\pgfsys@lineto{5.69089pt}{2.27637pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\odot\mathtt{y_{-}}$

9: $\textsc{Relinearize}(\mathtt{x_{+}})$

10: $\mathtt{y_{-}}\leftarrow\mathtt{y_{+}}$

11: $\mathtt{x_{-}}\leftarrow\mathtt{x_{+}}$

12: end for

13: return $\mathtt{x_{+}}$

14:end function

Algorithm 2 HE GD for an unconstrained QP

1:function
HEgdQP($\mathtt{Q},\mathtt{p},d,\lambda_{\min},\lambda_{\max},\mathtt{x_{0}},N$)

2: $\kappa\leftarrow\frac{\lambda_{\max}}{\lambda_{\min}}$

3: $\mathtt{x_{-}}\leftarrow\mathtt{x_{0}}$

4: $\eta\leftarrow\frac{-2}{\lambda_{\min}+\lambda_{\max}}$

5: for $t=0$ to $N-1$ do

6: $\mathtt{x_{+}}\leftarrow\mathtt{x_{-}}\mathord{\leavevmode\hbox
to8.9pt{\vbox to8.9pt{\pgfpicture\makeatletter\hbox{\hskip
0.75pt\lower-0.75pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\pgfsys@setlinewidth{1.5pt}\pgfsys@invoke{ } {}{{}}{}
{}{}{}\pgfsys@moveto{3.69899pt}{0.0pt}\pgfsys@lineto{3.69899pt}{7.398pt}\pgfsys@stroke\pgfsys@invoke{
} {}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{3.69899pt}\pgfsys@lineto{7.398pt}{3.69899pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\textsc{MMult}(\mathtt{Q},\mathtt{x_{-}},d,\eta)\mathord{\leavevmode\hbox
to8.9pt{\vbox to8.9pt{\pgfpicture\makeatletter\hbox{\hskip
0.75pt\lower-0.75pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\pgfsys@setlinewidth{1.5pt}\pgfsys@invoke{ } {}{{}}{}
{}{}{}\pgfsys@moveto{3.69899pt}{0.0pt}\pgfsys@lineto{3.69899pt}{7.398pt}\pgfsys@stroke\pgfsys@invoke{
} {}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{3.69899pt}\pgfsys@lineto{7.398pt}{3.69899pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\eta\odot\mathtt{p}$

7: $\mathtt{x_{+}}\leftarrow\mathtt{x_{-}}$

8: end for

9: return $\mathtt{x_{+}}$

10:end function

The only difference between the two algorithms is the presence of the extra
two $\odot$ and one $\mathord{\leavevmode\hbox to7.19pt{\vbox
to1.5pt{\pgfpicture\makeatletter\hbox{\hskip 0.75pt\lower 1.52637pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\pgfsys@setlinewidth{1.5pt}\pgfsys@invoke{ } {}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{2.27637pt}\pgfsys@lineto{5.69089pt}{2.27637pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}$ operations on line
8 of the (Algorithm 1) when compared to (Algorithm 2). These operations are
intrinsic to the accelerated gradient method as the method uses additional
past information to update to the next step.

###  IV-B Matrix multiplication seen differently

Because we would like to study gradient descent methods, let us consider a
simple matrix multiplication. Halevi and Shoup [16] introduced a method
(Algorithm 3) to evaluate an arbitrary linear transformation on encrypted
vectors. They exploit the diagonal encoding of a matrix to easily express the
matrix-vector multiplication by combining rotation and constant multiplication
operations

Algorithm 3 Halevi-Shoup LinTrans algorithm

function LinTrans($\mathtt{c},U$)

$n\leftarrow\dim(U)$

$\mathtt{cU}\leftarrow\mathtt{c}\odot u_{0}$

for $l=1$ to $n-1$ do

$\mathtt{cU}\leftarrow\mathtt{cU}\mathord{\leavevmode\hbox to8.9pt{\vbox
to8.9pt{\pgfpicture\makeatletter\hbox{\hskip 0.75pt\lower-0.75pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\pgfsys@setlinewidth{1.5pt}\pgfsys@invoke{ } {}{{}}{}
{}{}{}\pgfsys@moveto{3.69899pt}{0.0pt}\pgfsys@lineto{3.69899pt}{7.398pt}\pgfsys@stroke\pgfsys@invoke{
} {}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{3.69899pt}\pgfsys@lineto{7.398pt}{3.69899pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\textsc{Rot}(\mathtt{c},l)\odot
u_{l}$

end for

$\textsc{Relinearize}(\mathtt{cU})$

return $\mathtt{cU}$

end function

In [18] the authors elaborated further on the Halevi-Shoup method and proposed
a new matrix multiplication scheme (JKLS) that allows for a ciphered-matrix
multiplication that uses only one cipher-text per matrix following a row-
ordering encoding $A\rightarrow a$. Although convenient, the JKLS algorithm
needs 2 $\odot$ operations [18].

We propose a modified version of a matrix multiplication algorithm with 1 less
$\odot$ multiplication step:

$\displaystyle a_{k}$ $\displaystyle=V_{k}\odot a\ ,b_{k}=W_{k}\odot b\ ,\
k=0,\ldots,d-1$ $\displaystyle ab$ $\displaystyle=\sum_{k=0}^{d-1}a_{k}\bullet
b_{k}$

with:

$V_{k}(d\cdot i+j,l)=\begin{cases}1&\text{if }l=d\cdot i+[i+j+k]_{d}\\\
0&\text{otherwise}\\\ \end{cases}$ $W_{k}(d\cdot
i+j,l)=\begin{cases}1,&\text{if }l=d\cdot[i+j+k]_{d}+j\\\
0,&\text{otherwise}\\\ \end{cases}$

where $[\cdot]_{d}$ is a shortcut for $\cdot$ modulo $d$.

$V_{0}$$=$,$V_{1}$$=$, $\ldots$, $V_{3}$$=$$W_{0}$$=$,$W_{1}$$=$, $\ldots$,
$W_{3}$$=$

Figure 1: Matrix multiplication - $V_{k}$ and $W_{k}$ examples for $d=4$ Methodology |  No. of Ciphertexts |  Complexity |  Mult. depth |  Relineariza-tions
---|---|---|---|---
Halevi-Shoup |  $d$ |  $\mathcal{O}(d^{2})$ |  $1\bullet$ |  $1$
JKLS |  $1$ |  $\mathcal{O}(d)$ |  $1\bullet+2\odot$ |  $3$
Our work |  $1$ |  $\mathcal{O}(d)$ |  $1\bullet+1\odot$ |  $2$
TABLE II: Comparison of the different matrix multiplication algorithms.

The matrices $V_{k}$ and $W_{k}$ are permuting the row-encoded matrices $A$
and $B$ respectively such that the matrix multiplication algorithm as we know
can be implemented with element-wise multiplication and additions. Table II
summarizes the complexity differences of our method compared to Halevi-Shoup
and JKLS and Algorithm 4 the implementation of our methodology. With this
reduction of 1 $\odot$ operation in the matrix the multiplication algorithm we
are able to perform 9 and 6 iterations on GD and AGD respectively, as opposed
to 6 and 4 iterations if we were using the JKKS multiplication scheme.

Algorithm 4 HE Matrix Multiplication (MMult)

function MMult($\mathtt{A},\mathtt{B},d,a$)

$\mathtt{AB}\leftarrow\textsc{CipherText()}$

for $k=0$ to $d-1$ do

$\mathtt{A_{k}}\leftarrow\textsc{LinTrans}(\mathtt{A_{0}},V_{k}(a))$

$\mathtt{B_{k}}\leftarrow\textsc{LinTrans}(\mathtt{B_{0}},W_{k}(1))$

$\mathtt{AB_{k}}\leftarrow\mathtt{A_{k}}\bullet\mathtt{B_{k}}$

$\textsc{Relinearize}(\mathtt{AB_{k}})$

$\mathtt{AB}\leftarrow\mathtt{AB}\mathord{\leavevmode\hbox to8.9pt{\vbox
to8.9pt{\pgfpicture\makeatletter\hbox{\hskip 0.75pt\lower-0.75pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\pgfsys@setlinewidth{1.5pt}\pgfsys@invoke{ } {}{{}}{}
{}{}{}\pgfsys@moveto{3.69899pt}{0.0pt}\pgfsys@lineto{3.69899pt}{7.398pt}\pgfsys@stroke\pgfsys@invoke{
} {}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{3.69899pt}\pgfsys@lineto{7.398pt}{3.69899pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\mathtt{AB_{k}}$

end for

return $\mathtt{AB}$

end function

###  IV-C Extension to other QP problems

Extension to other QP problems is feasible. For instance, linear equality
constraints could be handled by converting the problem to an unconstrained QP,
or by solving primal-dual methods. These approaches sound completely viable
but are subject to the multiplication depth limitations on the HE circuit. In
other words, extra arithmetic operations can take place, but at the cost of
reducing the number of maximum iterations. Linear inequality constraints are
not directly supported, but an approach would be to decrypt and re-encrypt at
each iteration (not actually a practical solution).

##  V Numerical analysis

The HE resource requirements are directly proportional to the capacity of the
encryption circuit. The larger the circuit’s capacity, the larger the
computing memory and computational power required at each arithmetic
operation. Given our computing resources and the parameters of the Microsoft
SEAL [30], the largest circuit we can implement, in the CKKS scheme, has a
multiplication depth of 18. At each iteration, AGD and GD have a
multiplication depth of 3 and 2, resulting in a cap of 6 and 9 steps for AGD
and GD respectively. Even though we could not implement longer iterations due
to limitations on our computational resources, we believe that the results
would apply in that case too.

We start by running both algorithms with initial condition $x_{0}\neq x^{*}$
for the same matrix $Q$ and decrypt the outcomes at every iteration. Figure 2
shows that at each iteration the solution gets closer to the optimal $x^{*}$.

Figure 2: Decrypted HE-AGD steps for 100 repetitions with a 2-by-2 matrix with
$\kappa=2$, optimal value at $x^{*}=(1,1)$ initial condition $x_{0}=(3,3)$. A
similar behavior is observed for HE-GD.

Figure 3: Box plot of tolerance $f(x)-f(x^{*})$ for AGD (left) and GD (right)
with matrices of different dimensions (x axis) and different $\kappa$ (colored
bars).

As discussed in Section II AGD exhibits a superior convergence rate, hence it
allows meeting a given convergence (in terms of optimality) tolerance with
fewer iterations. However, in case of encryption, the computational limits
imposed by the allowable depth of the encryption circuit introduces a trade-
off, as AGD involves more arithmetic operations compared to GD (see Section
IV-A). As such, it might be computationally impossible to perform the number
of iterations needed by AGD to meet a given tolerance. We investigate this
trade-off numerically, and show that the preferred method depends on the
condition number $\kappa$ of the quadratic matrix $Q$.

To analyze this trade-off numerically we generate QP instances with condition
number $\kappa$ ranging from 1.5 to 50. To collect numerical statistics on the
effect of $\kappa$, for each $\kappa$ we generate 100 sets of randomly
generated symmetric positive-definite matrices $Q$ of dimension 2, 4 and
8333Higher dimensions are also feasible, and this is independent of the
multiplication depth limits., and associated random vectors $p$. We solve each
QP instance via the HE-GD and HE-AGD methods with an initial condition $x_{0}$
such that $\left\lVert x_{0}-x^{*}\right\rVert_{2}=1$. Our goal is to
investigate which algorithm achieves better tolerance values (in terms
distance to the optimal value) at the last iteration allowed by the
encryption’s circuit depth. The latter is iteration 6 for AGD and iteration 9
for GD. Figure 3 illustrates the distribution of the tolerance $f(x)-f(x^{*})$
(optimality gap of the returned solution $x$ from the optimal cost $f(x^{*}$)
for AGD (top) and GD (bottom) with matrices of different dimensions (x axis)
and different values of the condition number $\kappa$ (color code).

Table III highlights the important observations stemming from Figure 3. In
particular, GD profits from the extra iterations and achieves better tolerance
(getting closer to the optimal) values when $\kappa\leq 5$ (upper table). Yet,
AGD outperforms GD in cases where $\kappa>5$ (lower table) although with worse
tolerance values (that is, further form the optimal). This shows that the
limitation on multiplication depth is a major issue when handling matrices
with high values of $\kappa$.

$d$ $\kappa$ | 1.5 | 2 | 3 | 5
---|---|---|---|---
2 | $3\cdot 10^{-9}$ | $4\cdot 10^{-9}$ | $3\cdot 10^{-7}$ | $5\cdot 10^{-5}$
4 | $1\cdot 10^{-8}$ | $1\cdot 10^{-8}$ | $8\cdot 10^{-8}$ | $1\cdot 10^{-5}$
8 | $6\cdot 10^{-8}$ | $4\cdot 10^{-8}$ | $7\cdot 10^{-8}$ | $5\cdot 10^{-6}$

$d$ $\kappa$ | 10 | 20 | 50
---|---|---|---
2 | $7\cdot 10^{-3}$ | $2\cdot 10^{-3}$ | $5\cdot 10^{-3}$
4 | $2\cdot 10^{-4}$ | $8\cdot 10^{-4}$ | $2\cdot 10^{-3}$
8 | $6\cdot 10^{-5}$ | $2\cdot 10^{-4}$ | $9\cdot 10^{-4}$

TABLE III: Comparison of AGD (6th iteration) against GD (9th iteration) for
matrices of different sizes and conditional values $\kappa$. The upper and
lower table indicate the range of values for $\kappa$ for which the GD and the
AGD algorithm, respectively, are preferable in terms of returning a solution
closer to the optimal one. The numbers represent the median tolerance level
$f(x)-f(x^{*})$ out of the 100 repetitions corresponding to different matrices
$Q$.

The code running the numerical examples presented here
(https://github.com/f2cf2e10/agd-he) used our own Python wapper of the
Microsoft SEAL [30] C++ library (https://github.com/f2cf2e10/pSEAL). We used
an Intel Xeon E5-1620 with 24GB of RAM machine running Debian 11.

##  VI Conclusion

In this paper we studied, implemented, and analyzed both gradient and
accelerated gradient descent algorithms to solve QP problems in a HE fashion.
We evaluated different encryption schemes (BFV, BGV, and CKKS) and argued that
CKKS is the only suitable scheme for implementing gradient descent methods as
it allows for freedom in the choice of the algorithm’s step size. In our
implementation, AGD takes an extra multiplication operation at each step when
compared to GD. As a result, AGD cannot run for as many steps as GD for an
encryption circuit with the same security parameters and channel capacity
(i.e. multiplication depth). We demonstrate that the condition number of the
matrix of the quadratic term of the objective function plays an important role
in determining which algorithm is preferred. For higher values of the
condition number, AGD is preferred, as it converges faster to the solution,
even if performing fewer iterations. Whereas for lower values of the condition
number of the matrix GD performs better due to the extra iterations. In
addition, we proposed a new HE matrix multiplication algorithm that is more
efficient, in terms of multiplication depth, than other algorithms proposed in
the literature. These observations have been verified by means of a numerical
investigation. Yet, there are still many outstanding challenges in HE version
of iterative numerical procedures. For instance, solving for constrained
problems is in general a challenge, given the extra operations required to
project into the constrained set. To solve these issues, further work should
focus on optimizing for the multiplication depth and perhaps explore
alternative ways to encode matrices in the HE scheme.

## References

  * [1] Andreea B. Alexandru, Konstantinos Gatsis, Yasser Shoukry, Sanjit A. Seshia, Paulo Tabuada, and George J. Pappas.  Cloud-based quadratic optimization with partially homomorphic encryption.  IEEE Transactions on Automatic Control, 66(5):2357–2364, 2021.
  * [2] Andreea B. Alexandru, Anastasios Tsiamis, and George J. Pappas.  Encrypted distributed lasso for sparse data predictive control, 2021.
  * [3] Zvika Brakerski.  Fully homomorphic encryption without modulus switching from classical gapsvp.  In Reihaneh Safavi-Naini and Ran Canetti, editors, Advances in Cryptology – CRYPTO 2012, pages 868–886, Berlin, Heidelberg, 2012. Springer Berlin Heidelberg.
  * [4] Zvika Brakerski, Craig Gentry, and Vinod Vaikuntanathan.  (leveled) fully homomorphic encryption without bootstrapping.  In Proceedings of the 3rd Innovations in Theoretical Computer Science Conference, ITCS ’12, page 309–325, New York, NY, USA, 2012. Association for Computing Machinery.
  * [5] Sébastien Bubeck.  Convex optimization: Algorithms and complexity.  Found. Trends Mach. Learn., 8(3–4):231–357, November 2015.
  * [6] Jung Hee Cheon, Kyoohyung Han, Hyuntae Kim, Junsoo Kim, and Hyungbo Shim.  Need for controllers having integer coefficients in homomorphically encrypted dynamic system.  In 2018 IEEE Conference on Decision and Control (CDC), pages 5020–5025, 2018.
  * [7] Jung Hee Cheon, Andrey Kim, Miran Kim, and Yongsoo Song.  Homomorphic encryption for arithmetic of approximate numbers.  In Tsuyoshi Takagi and Thomas Peyrin, editors, Advances in Cryptology – ASIACRYPT 2017, pages 409–437, Cham, 2017. Springer International Publishing.
  * [8] Jung Hee Cheon, Dongwoo Kim, Duhyeong Kim, Hun Hee Lee, and Keewoo Lee.  Numerical method for comparison on homomorphically encrypted numbers.  Cryptology ePrint Archive, Paper 2019/417, 2019.  https://eprint.iacr.org/2019/417.
  * [9] T. Elgamal.  A public key cryptosystem and a signature scheme based on discrete logarithms.  IEEE Transactions on Information Theory, 31(4):469–472, 1985.
  * [10] Junfeng Fan and Frederik Vercauteren.  Somewhat practical fully homomorphic encryption.  Cryptology ePrint Archive, Report 2012/144, 2012.  https://eprint.iacr.org/2012/144.
  * [11] Farhad Farokhi, Iman Shames, and Nathan Batterham.  Secure and private cloud-based control using semi-homomorphic encryption.  IFAC-PapersOnLine, 49(22):163–168, 2016.  6th IFAC Workshop on Distributed Estimation and Control in Networked Systems NECSYS 2016.
  * [12] Farhad Farokhi, Iman Shames, and Nathan Batterham.  Secure and private control using semi-homomorphic encryption.  Control Engineering Practice, 67:13–20, 2017.
  * [13] Craig Gentry.  Fully homomorphic encryption using ideal lattices.  In Proceedings of the Forty-First Annual ACM Symposium on Theory of Computing, STOC ’09, page 169–178, New York, NY, USA, 2009. Association for Computing Machinery.
  * [14] Craig Gentry.  Computing arbitrary functions of encrypted data.  Commun. ACM, 53(3):97–105, March 2010.
  * [15] Craig Gentry, Amit Sahai, and Brent Waters.  Homomorphic encryption from learning with errors: Conceptually-simpler, asymptotically-faster, attribute-based.  In Ran Canetti and Juan A. Garay, editors, Advances in Cryptology – CRYPTO 2013, pages 75–92, Berlin, Heidelberg, 2013. Springer Berlin Heidelberg.
  * [16] Shai Halevi and Victor Shoup.  Faster homomorphic linear transformations in helib.  Cryptology ePrint Archive, Report 2018/244, 2018.  https://ia.cr/2018/244.
  * [17] Ilia Iliashenko and Vincent Zucca.  Faster homomorphic comparison operations for bgv and bfv.  Cryptology ePrint Archive, Paper 2021/315, 2021.  https://eprint.iacr.org/2021/315.
  * [18] Xiaoqian Jiang, Miran Kim, Kristin Lauter, and Yongsoo Song.  Secure outsourced matrix computation and application to neural networks.  In Proceedings of the 2018 ACM SIGSAC Conference on Computer and Communications Security, CCS ’18, page 1209–1222, New York, NY, USA, 2018. Association for Computing Machinery.
  * [19] Andrey Kim, Antonis Papadimitriou, and Yuriy Polyakov.  Approximate homomorphic encryption with reduced approximation error.  Cryptology ePrint Archive, Report 2020/1118, 2020.  https://ia.cr/2020/1118.
  * [20] Junsoo Kim, Dongwoo Kim, Yongsoo Song, Hyungbo Shim, Henrik Sandberg, and Karl H. Johansson.  Comparison of encrypted control approaches and tutorial on dynamic systems using learning with errors-based homomorphic encryption.  Annual Reviews in Control, 54:200–218, 2022.
  * [21] Junsoo Kim, Chanhwa Lee, Hyungbo Shim, Jung Hee Cheon, Andrey Kim, Miran Kim, and Yongsoo Song.  Encrypting controller using fully homomorphic encryption for security of cyber-physical systems.  IFAC-PapersOnLine, 49(22):175–180, 2016.  6th IFAC Workshop on Distributed Estimation and Control in Networked Systems NECSYS 2016.
  * [22] Junsoo Kim, Hyungbo Shim, and Kyoohyung Han.  Dynamic controller that operates over homomorphically encrypted data for infinite time horizon.  IEEE Transactions on Automatic Control, 68(2):660–672, 2023.
  * [23] Shane Kosieradzki, Yingxin Qiu, Kiminao Kogiso, and Jun Ueda.  Rewrite rules for automated depth reduction of encrypted control expressions with somewhat homomorphic encryption.  In 2022 IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM), pages 804–809, 2022.
  * [24] Yongwoo Lee, Joonwoo Lee, Young-Sik Kim, HyungChul Kang, and Jong-Seon No.  High-precision and low-complexity approximate homomorphic encryption by error variance minimization.  Cryptology ePrint Archive, Report 2020/1549, 2020.  https://ia.cr/2020/1549.
  * [25] Vadim Lyubashevsky, Chris Peikert, and Oded Regev.  On ideal lattices and learning with errors over rings.  Cryptology ePrint Archive, Paper 2012/230, 2012.  https://eprint.iacr.org/2012/230.
  * [26] Pascal Paillier.  Public-key cryptosystems based on composite degree residuosity classes.  In Jacques Stern, editor, Advances in Cryptology — EUROCRYPT ’99, pages 223–238, Berlin, Heidelberg, 1999. Springer Berlin Heidelberg.
  * [27] Moritz Schulze Darup, Andreea B. Alexandru, Daniel E. Quevedo, and George J. Pappas.  Encrypted control for networked systems: An illustrative introduction and current challenges.  IEEE Control Systems Magazine, 41(3):58–78, 2021.
  * [28] Moritz Schulze Darup, Adrian Redder, and Daniel E. Quevedo.  Encrypted cooperative control based on structured feedback.  IEEE Control Systems Letters, 3(1):37–42, 2019.
  * [29] Moritz Schulze Darup, Adrian Redder, Iman Shames, Farhad Farokhi, and Daniel Quevedo.  Towards encrypted mpc for linear constrained systems.  IEEE Control Systems Letters, 2(2):195–200, 2018.
  * [30] Microsoft SEAL (release 3.7).  https://github.com/Microsoft/SEAL, September 2021.  Microsoft Research, Redmond, WA.