# Privacy-Preserving Distributed Machine Learning Made Faster

Zoe L. Jiang, Jiajing Gu, Hongxiao Wang, Yulin Wu, Junbin Fang, Siu-Ming Yiu,
Wenjian Luo, Xuan Wang

###### Abstract

With the development of machine learning, it is difficult for a single server
to process all the data. So machine learning tasks need to be spread across
multiple servers, turning the centralized machine learning into a distributed
one. However, privacy remains an unsolved problem in distributed machine
learning. Multi-key homomorphic encryption is one of the good candidates to
solve the problem. However, the most recent result about Multi-key homomorphic
encryption scheme (MKTFHE) only supports NAND gate. Although it is Turing
complete, it requires efficient encapsulation of NAND gate to further support
mathematical calculation. This paper designs and implements a series of
operations on positive and negative integers accurately. First, we design
basic bootstrapped gates with the same efficiency as that of the NAND gate.
Second, we construct practical $k$-bit complement mathematical operators based
on our basic binary bootstrapped gates. The constructed created can perform
addition, subtraction, multiplication and division on both positive and
negative integers. Finally we demonstrated the generality of the designed
operators by achieving a distributed privacy-preserving machine learning
algorithm, i.e. linear regression with two different solutions. Experiments
show that the operators we designed are practical and efficient.

###### keywords:

Privacy-preserving machine learning , Distributed machine learning , Multi-key
fully homomorphic encryption , Practical complement operator

##  1 Introduction

In the big data era, it is necessary to spread the machine learning tasks
across multiple servers and transform centralized systems into distributed
ones [1]. These distributed systems present new challenges and especially
privacy remains an unsolved challenge in it [2].

Privacy computing is a kind of techniques which perform data computation
without specified information leakage. Homomorphic encryption is a special
form of encryption that permits users to perform computations on encrypted
data without first decrypting it, which has great practical implications in
the outsourcing of private computations. It is called fully homomorphic
encryption (FHE) if both addition and multiplication are supported. The
security of FHE is usually guaranteed by the hard problems on lattices, which
is also considered as one of the post-quantum cryptographic schemes [3].
However, computing on encrypted data using such FHE schemes will increase the
noise in the ciphertext. If the noise in the ciphertext grows beyond some
threshold, it will result in wrong decryption [4]. This problem can be solved
by a bootstrapping technique, first proposed by Gentry [5]. Bootstrapping can
reduce the noise in ciphertext to refresh it so that the homomorphic
evaluation can be performed constantly.

Obviously, _single-key_ fully homomorphic encryption only allows server to
perform addition and multiplication on data encrypted by same key. _Multi-key_
fully homomorphic encryption (MKFHE) was proposed to circumvent this
shortcoming [6]. It enables users to encrypt their own data under their own
keys, while homomorphic evaluations can be performed on encrypted data
directly at server side without decryption. It avoids the possibility that
user and server conspire to steal the data of other users. Therefore, MKFHE
indeed realizes secure multi-party computation with untrusted party. While
Chen et al. [7] developed the library Multi-key fully homomorphic encryption
over torus (MKTFHE) for NAND gate.

However, MKTFHE library only provides multi-key homomorphic NAND gates.
Perform multi-key homomorphic computation using only homomorphic NAND gates is
a complex and time-consuming work. As a result, a simple multi-key homomorphic
NAND gate is not practical to be used directly. In order to better apply
multi-key homomorphic encryption on privacy computing, we design and implement
a series of practical multi-key homomorphic mathematical operators.

In this study, we make the following contributions :

  1. 1.

We designed a series of fundamental multi-key bootstrapped gates based on
MKTFHE, including multi-key bootstrapped AND gate, multi-key bootstrapped OR
gate, multi-key bootstrapped NOR gate, multi-key bootstrapped XOR gate, multi-
key bootstrapped XNOR gate, and multi-key NOT gate (without gate
bootstrapping). Experiment results show that the proposed fundemental multi-
key bootstrapped gates are more efficient than directly joining multi-key
bootstrapped NAND gates to build basic binary gates.

  2. 2.

We designed a series of multi-key homomorphic operators based on the basic
multi-key bootstrapped gates we constructed. The operators we designed include
$k$-bit complement array integer adder, $k$-bit complement array integer
subtractor, $k$-bit complement array integer multiplier and $k$-bit complement
array integer divider, which can perform addition, subtraction, multiplication
and division on arbitrary bits integers in both positive and negative. We
construct the complement array adder linearly and used a structure similar to
the 4-bit array integer multiplier in [8] to construct our $k$-bit complement
multiplier. In the similar way we construct the $k$-bit complement array
integer divider. The subtractor is constructed with a adder and a simple
multiplier. The time of the adder and the subtractor grows linearly with the
bits of input numbers, while the time of the multiplier and divider grows
quadratically. The time of the divider grows almost linearly with the number
of layers in divider.

  3. 3.

We train linear regression model by utilizing our proposed multi-key
homomorphic operators. Taking linear regression as an example, we implement a
whole multi-key fully homomorphic machine learning scheme. Experimental
results show that the practicability and generality of our homomorphic
operator and the training time grows linearly with the number of participants.

The rest of this paper is organized as follows. Section 2 discusses the
related research work. In Section 3, it clarifies the notation and review some
constructions related. Section 4 describes the $k$-bit complement array
operators with implementation based on MKTFHE and the distributed machine
learning model we evaluate, followed by the performance analysis and
experimental results of our implementation in Section 5. Section 6 concludes
the paper with future work.

##  2 Related Work

Previous research efforts have made a number of contributions to the
development of privacy-preserving machine learning. Some work has focused on
securely outsourcing the training of ML models to the cloud, typically by
using homomorphic encryption techniques [9]. In 2016, Aono et al. [10]
proposed a secure system for protecting the training data in logistic
regression via homomorphic encryption. In 2018, Crawford et al. [11] built a
system that uses fully-homomorphic encryption to approximate the coefficients
of a logistic-regression model built from genomic data, while Kim et al. [12,
13] presented a method to train a logistic regression model without
information leakage. However, these scheme are all based on simple-key fully
homomorphic encryption, which are not suitable for distributed machine
learning models.

In 2012, Graepel et al. [14] first proposed that it is possible to perform
machine learning algorithm on encrypted data by homomorphic encryption scheme.
Then, several single-key homomorphic encryption schemes are used for machine
learning prediction [15, 16, 17, 18, 19] or machine learning training [20, 21,
22]. A few machine learning algorithms has adapted multi-key homomorphic
encryption and demonstrated the premising future [23].

Multi-key fully homomorphic encryption was first proposed based on the NTRU
assumption by Lopez et al. [6] at 2012, which is intended to apply to on-the-
fly multiparty computation. In 2015, Clear et al. [24] constructed the first
LWE-based MKFHE scheme, and improved by Mukherjee and Wichs [25] in 2016.
These schemes are single-hop MKFHE schemes, which means all the participants
must be known in advance. This problem was solved by Peikert et al. [26] and
Brakerski et al. [27] in 2016 by constructing multi-hop MKFHE schemes.
However, these schemes are impractical and without implementation.

The first implementation of the MKFHE scheme was achieved by Chen et al. [7]
in 2019, named MKTFHE. In their scheme, a multi-key bootstrapped NAND gate is
described. This scheme was improved by Lee and Park [28] in 2019, realizing
the distributed decryption in MKTFHE. However, using multi-key bootstrapped
NAND gates to construct homomorphic encryption operation directly has
disadvantages such as low practicability, high complexity and error-prone
construction process. At the same time, the complex structure lead to
expensive optimization.

##  3 Preliminaries

###  3.1 Notation

The rest of the paper uses the following notations. $\mathbb{T}$ denotes the
real Torus $\mathbb{R}/\mathbb{Z}$, the set of real numbers modulo $1$.
$\mathbb{T}_{N}[X]$ denotes $\mathbb{R}[X]/(X^{N}+1)$ mod $1$. $k$ represents
the number of parties. TLWE is used to denote the (scalar) binary learning
with error problem over Torus, while for the ring mode, we use the notation of
TRLWE. $params$ represents the parameter sets used in TFHE (fully homomorphic
encryption over torus) scheme, while $mkparams$ represents the parameter sets
used in MKTFHE scheme and our scheme.

###  3.2 MKTFHE

MKTFHE scheme is the multi-key version of TFHE scheme. TFHE, constructed by
Chillotti et al. [29, 30, 31], is a fast fully homomorphic encryption (FHE)
scheme over the torus, which generalizes and improves the FHE based on GSW
[32] and its ring variants. In the TFHE scheme, bootstrapped binary gates are
designed to represent the functions developers need.

The main idea of TFHE is to bootstrap after every binary gate evaluation to
refresh the ciphertext in order to make it usable for the following
operations, resulting in that arbitrarily deep circuits can be homomorphically
evaluated. That is to say, a single parameter set allows the sever to evaluate
any function. The entire homomorphic evaluation of the circuits will take time
proportional to the number of the binary gates used or, if parallelism is
involved, to the number of the circuit layers.

The message space of TFHE bootstrapping gates is $\mathbb{T}$. A TLWE
ciphertext $(a,b)\in T^{n+1}$ encrypted a message $\mu\in\mathbb{T}$ with
noise parameter $\alpha$.

In the TFHE scheme, the homomorphic evaluation of a binary gate is achieved
with operations between TLWE samples and a gate-bootstrapping just after that
(except the bootstrapping NOT gate, which does not need bootstrapping). By
using this approach, all the basic gates can be evaluated with a single
bootstrapping process (GB):

  * 1.

$\mathsf{TFHE.BootsNAND}(c_{1},c_{2})=\mathsf{GB}((0,\frac{5}{8})-c_{1}-c_{2})$

  * 2.

$\mathsf{TFHE.BootsAND}(c_{1},c_{2})=\mathsf{GB}((0,-\frac{1}{8})+c_{1}+c_{2})$

  * 3.

$\mathsf{TFHE.BootsOR}(c_{1},c_{2})=\mathsf{GB}((0,\frac{1}{8})+c_{1}+c_{2})$

  * 4.

$\mathsf{TFHE.BootsXOR}(c_{1},c_{2})=\mathsf{GB}(2\cdot(c_{1}-c_{2}))$

  * 5.

$\mathsf{TFHE.NOT}(c)=(0,\frac{1}{4})-c$

The TFHE scheme has the advantages of fast bootstrapping, efficient
homomorphic logic circuit evaluation, and so on. Its multi-key version, named
MKTFHE [7], was constructed by Chen et al. in 2019. MKTFHE is the first
attempt in the literature to implement an MKFHE scheme in codes.

In the MKTFHE scheme, the ciphertext length increases linearly with the number
of users, and a homomorphic NAND gate with bootstrapping is given. The MKTFHE
scheme is comprised of the following algorithms.

  * 1.

$\mathsf{MKTFHE.Setup}(1^{\lambda})$: The cloud equipment take a security
parameter $\lambda$ as input, and output the public parameter set $mkparams$.

  * 2.

$\mathsf{MKTFHE.KeyGen}(mkparams)$: Each party generates its keys
independently. First sample the TLWE secret key $sk_{i}$. Then the algorithm
set public key $pk_{i}$, bootstrapping key $BK_{i}$ and key-switching key
$KS_{i}$.

  * 3.

$\mathsf{MKTFHE.SymEnc}(\mu)$: This algorithm encrypts an input bit
$\mu\in\\{0,1\\}$, and returns a TLWE ciphertext with the scaling factor
$\frac{1}{4}$. The output ciphertext $c=(b,a)\in\mathbb{T}^{n+1}$ satisfies
$b+<a,s>\approx\frac{1}{4}\mu$.

  * 4.

$\mathsf{MKTFHE.SymDec}(c,\\{sk_{i}\\})$: Input a TLWE ciphertext
$c=(b,a_{i},...,a_{k})$ and a set of secret keys $\\{sk_{i}\\}$, and return
the message $\mu$ which minimizes
$|b+\sum_{i=1}^{k}<a_{i},s_{i}>-\frac{1}{4}m|$.

  * 5.

$\mathsf{MKTFHE.MKKeySwitch}(c,\\{KS_{i}\\}_{i\in[k]})$: Input an expa-nded
ciphertext $c$ corresponding to $t=(t_{1},...,t_{k})$ and a sequence of key-
switching keys from $t_{i}$ to $s_{i}$ and returns an encryption of the same
message under $s=(s_{1},...,s_{k})$.

  * 6.

$\mathsf{MKTFHE.NAND}(c_{1},c_{2},\\{pk_{i},BK_{i},KS_{i}\\}_{i\in[k]})$: Take
two TLWE ciphertext as input. This algorithm first extends two input
ciphertexts and evaluate the NAND gate homomorphically on encrypted bits. Then
the algorithm evaluates the decryption circuit of the TLWE ciphertext and run
the $\mathsf{MKTFHE.MKKeySwitch}(c,\\{KS_{i}\\}_{i\in[k]})$ algorithm.
Finally, the algorithm outputs a TLWE ciphertext encrypting $m=m_{1}\barwedge
m_{2}$.

###  3.3 Complement operators

Assume that there are two $k$-bit inputs $a_{in}$ and $b_{in}$, a $k$-bit
complement array integer adder can be constructed in following steps.

  1. 1.

Construct a semi-adder with two XOR gates: Input two addends $a$ and $b$ and
the carry $c$, return $out=(c\bigoplus(a\bigoplus b))$.

  2. 2.

Construct a carrier with three AND gates and two OR gates: Input two addends
$a$ and $b$ and the carry $c$, return $c=(a\wedge b)\vee(a\wedge
c)\vee(b\wedge c)$.

  3. 3.

Construct a $1$-bit adder with a semi-adder and a carrier: Input two addends
$a$ and $b$ and the carry $c$, return $out=(c\bigoplus(a\bigoplus b))$ and
$c=(a\wedge b)\vee(a\wedge c)\vee(b\wedge c)$.

  4. 4.

Construct a $k$-bit complement array integer adder with k $1$-bit adder: Input
two addends $a_{in}$ and $b_{in}$, set the carry $c[0]=0$, compute
$out[i]=(c[i]\bigoplus(a_{in}[i]\bigoplus b_{in}[i]))$,
$cout[i+1]=(a_{in}[i]\wedge b_{in}[i])\vee(a_{in}[i]\wedge
c[i])\vee(b_{in}[i]\wedge c[i])$ with $i=\\{0,1,...,k-1\\}$. Return the result
$out$.

Figure 1: The construction of $k$-bit complement array integer multiplier.The
index of the adder stands for the type of it (0-adder, $1$-adder or 2-adder).

Assume that there are two $k$-bit inputs $a_{in}$ and $b_{in}$, a $k$-bit
complement array integer multiplier can be constructed in following steps.

  1. 1.

Prepare three kinds of $1$-bit adder for complement multiplier:

    1. (a)

0-adder: The same as $1$-bit adder above.

    2. (b)

$1$-adder: Input two addends $a$ and $b$ and the carry $c$, return
$out=\overline{(c\bigoplus(\overline{a}\bigoplus b))}$ and
$out=(\overline{a}\wedge b)\vee(\overline{a}\wedge c)\vee(b\wedge c)$.

    3. (c)

2-adder: Input two addends $a$ and $b$ and the carry $c$, return
$out=(c\bigoplus(\overline{a}\bigoplus\overline{b}))$ and
$out=\overline{(\overline{a}\wedge\overline{b})\vee(\overline{a}\wedge c)}$
$\overline{\vee(\overline{b}\wedge c)}$.

  2. 2.

Construct a $k$-bit complement array integer multiplier according to the
following rules, as shown in Figure 1.

    1. (a)

Arrange the adders in $k$ rows and $k-1$ columns. The input $a$ and $b$ of the
adder in $i$ row and $j$ column are $a_{in}[j]\wedge b_{in}[i-1]$ and
$a_{in}[j-1]\wedge b_{in}[i]$. The inputs $c$ of the adders are $0$ in the
first line while the $cout$ of adders above in other lines.

    2. (b)

The highest bit of each addend (that is to say, $a_{in}[k-1]$ and
$b_{in}[k-1]$) has the weight $-1$.

    3. (c)

If two input addends have no weight then use 0-adder.

    4. (d)

If one of the input addends has the weight $-1$ after AND gate or $1$-adder
then use $1$-adder.

    5. (e)

If one of the addends is an output of a $1$-adder and the other addend has the
weight $-1$ after AND gate, use 2-adder.

    6. (f)

If one of the addends is an output of a 2-adder, use 2-adder.

Obviously a subtractor can be constructed by an adder and a multiplier. Assume
that there are two $k$-bit inputs $a_{in}$ and $b_{in}$, a $k$-bit complement
array integer subtractor can be calculated by: $out=a_{i}n+(-1)\times b_{i}n$.

Assume that there are two $k$-bit inputs $a_{in}$ and $b_{in}$, a $k$-bit
complement array integer divider [33] can be constructed in following steps.

  1. 1.

Prepare the Controlled Adder/subtractor (CAS) for the divider. We use three
bootstrapped XOR gates, two bootstrapped OR gates and two bootstrapped AND
gates to construct the CAS. The CAS takes $a_{in}[i]$, $b_{in}[i]$,
$c_{in}[i]$ and $p$ as input, and output $out[i]$ and the carrier $c[i+1]$.
The construction of a CAS is shown in Figure 2. The CAS will perform addition
if $p=0$ and perform subtraction if $p=1$.

Figure 2: The construction of CAS

  2. 2.

Then design a absolute value array divider with the CAS above. It takes the
$2k$-bit dividend $a_{in}$ and the $k$-bit divisor $b_{in}$ as input, and
output the quotient $q$ and the remainder $r$.

  3. 3.

XOR gates are used to decide the sign bit of the quotient to realize the
division of both positive and negative integers.

  4. 4.

In order to facilitate the operation, the input and output of the divider is
unified to complement format. As a result, a compensation device is needed
while processing the input and the output. To realize the complement, we first
perform XOR gates on the sign bit and each bit, and then add the sign bit to
the result, and finally get the complement code.

  5. 5.

Construct the $k$-bit complement array divider with a absolute value array
divider, a XOR gate and two compensation device. The $k$-bit complement array
divider takes the $2k$-bit complement dividend $a_{in}$ and the $k$-bit
complement divisor $b_{in}$ as input, and output the complement quotient $q$
and the complement remainder $r$. The structure of the $k$-bit complement
array divider is shown in Figure 3. For simplicity, we set $k=3$in this
figure.

Figure 3: The construction of $k$-bit complement array divider

###  3.4 Linear regression

In unary linear regression, there is a calculation formula to obtain the
global optimal solution during training the model. According to the definition
of linear regression, we can evaluate the accuracy of the model. Suppose we
observe $n$ data pairs and call them $\\{(x_{i},y_{i}),i=1,...,n\\}$. We
implement unary linear regression and utilize its calculation formula to train
model:

$(\omega^{*},b^{*})=argmin_{(\omega,b)}\sum_{i=1}^{m}(f(x_{i})-y_{i})^{2}$ (1)

By deriving the above formula and making the derivative zero, we can get the
calculation formula.

$\omega=\frac{\sum_{i=1}^{m}y_{i}(x_{i}-\bar{x})}{\sum_{i=1}^{m}x_{i}^{2}-\frac{1}{m}(\sum_{i=1}^{m}x_{i})^{2}},\bar{x}=\frac{1}{m}\sum_{i=1}^{m}x_{i}$
(2)

$b=\frac{1}{m}\sum_{i=1}^{m}(y_{i}-\omega x_{i})$ (3)

Where $\omega$ is the slope of the line and $b$ is the y-intercept. $\bar{x}$
stands for the average of $x_{1},x_{2},\cdots,x_{m}$.

In most linear regression, there is often no calculation formula to help us
obtain the global optimal solution, but the common practice is to utilize the
Gradient descent (GD) method to obtain the local optimal through multiple
iterations. The iterative formula of the classical gradient descent method is
as follows.

$\theta_{j}^{t+1}=\theta_{j}^{t}-\alpha\frac{\partial
J(\theta_{j}^{t})}{\theta_{j}^{t}}$ (4)

$J(\theta)=\frac{1}{2m}\sum_{i=1}^{m}(h_{\theta}(x_{i})-y_{i})^{2}$ (5)

Where $\alpha$ is the learning rate, $t$ is the number of iterations, $\theta$
is the parameter of model, $J(\theta)$ is the gradient of $\theta$, $j$ is the
index of parameters to be solved, $m$ is the number of training data and
$h_{\theta}$ stands for the prediction model.

##  4 The proposed privacy-preserving distributed machine learning

Although MKTFHE scheme provides bootstrapped NAND gates, it is time-consuming
to use NAND gates directly for privacy computing on Industrial Internet. We
expand NAND gates to other fundamental gates to better support mathematical
operators. At the same time, we construct an adder, a subtractor, a multiplier
and a divider, so that users can perform mathematical operators infinitely.

###  4.1 Homomorphic fundamental gates

In previous works, the other fundamental gates (except the NAND gate) are
constructed by joining the multi-key bootstrapped NAND gates together.
Although only the multi-key bootstrapped NAND gate is described in [8], any
arbitrary binary bootstrapped gate can be evaluated in the same way. Our
scheme uses a similar approach in MKTFHE to directly evaluate the following
basic gates in a multi-key version.

  * 1.

$\mathsf{MKBootsAND}(c_{1},c_{2},\\{pk_{i}\\}_{i\in[k]})$: Given ciphertexts
$c_{1}$ and $c_{2}$, this algorithm first extends two input ciphertexts into
$c_{1}^{{}^{\prime}}$ and $c_{2}^{{}^{\prime}}$ under the same key. Then it
evaluates the AND gate by computing:

$c^{{}^{\prime}}=(0,-\frac{1}{8})+c_{1}+c_{2}$ (6)

homomorphically on encrypted bits. Then this algorithm evaluates the
decryption circuit of the ciphertext $c^{{}^{\prime}}$ to bootstrap it.
Finally, this algorithm run the $\mathsf{MKTFHE.MKKey}$
$\mathsf{Switch(c,{KS_{i}}_{i\in[k]})}$ algorithm and outputs a TLWE
ciphertext encrypting $m=m_{1}\wedge m_{2}$.

  * 2.

$\mathsf{MKBootsOR}(c_{1},c_{2},\\{pk_{i}\\}_{i\in[k]})$: Given ciphertexts
$c_{1}$ and $c_{2}$, this algorithm first extends two input ciphertexts into
$c_{1}^{{}^{\prime}}$ and $c_{2}^{{}^{\prime}}$ and evaluates the OR gate by
computing:

$c^{{}^{\prime}}=(0,\frac{1}{8})+c_{1}+c_{2}$ (7)

homomorphically on encrypted bits. Then this algorithm bootstraps it, run the
$\mathsf{MKTFHE.MKKeySwitch(c,{KS_{i}}_{i\in[k]})}$ algorithm and outputs a
TLWE ciphertext encrypting $m=m_{1}\vee m_{2}$.

  * 3.

$\mathsf{MKNOT}(c,mkparams)$: Take ciphertext $c$ as input, this algorithm
evaluates the NOT gate by computing:

$c^{{}^{\prime}}=(0,\frac{1}{4})-c$ (8)

homomorphically. This computation will not increase the noise in ciphertext or
change the key so that no bootstrapping and key-switching process is needed.

  * 4.

$\mathsf{MKBootsNOR}(c_{1},c_{2},\\{pk_{i}\\}_{i\in[k]})$: Given ciphertexts
$c_{1}$ and $c_{2}$, this algorithm first extends two input ciphertexts into
$c_{1}^{{}^{\prime}}$ and $c_{2}^{{}^{\prime}}$ and evaluates the OR gate by
computing:

$c^{{}^{\prime}}=(0,\frac{1}{8})-c_{1}-c_{2}$ (9)

homomorphically on encrypted bits. Then this algorithm bootstraps it, run the
$\mathsf{MKTFHE.MKKeySwitch(c,{KS_{i}}_{i\in[k]})}$ algorithm and outputs a
TLWE ciphertext encrypting $m=m_{1}\bar{\vee}m_{2}$.

  * 5.

$\mathsf{MKBootsXOR}(c_{1},c_{2},\\{pk_{i}\\}_{i\in[k]})$: Given ciphertexts
$c_{1}$ and $c_{2}$, this algorithm first extends two input ciphertexts and
evaluates the XOR gate by computing:

$c^{{}^{\prime}}=2\cdot(c_{1}-c_{2})$ (10)

homomorphically. Then run the bootstrapping and key-switching algorithm, and
outputs a TLWE ciphertext encrypting $m=m_{1}\bigoplus m_{2}$.

  * 6.

$\mathsf{MKBootsXNOR}(c_{1},c_{2},\\{pk_{i}\\}_{i\in[k]})$: Given ciphertexts
$c_{1}$ and $c_{2}$, this algorithm first extends two input ciphertexts and
evaluates the XOR gate by computing:

$c^{{}^{\prime}}=\frac{1}{4}-2\cdot(c_{1}-c_{2})$ (11)

homomorphically. Then run the bootstrapping and key-switching algorithm, and
output a TLWE ciphertext encrypting $m=m_{1}\bigodot m_{2}$.

First of all, a mapping function $ModToT$ is required to transform a message
bit $m$ to an element on torus, like $(0,\frac{1}{8})$ in Equation (6), such
that all of the following computations are executed on torus. The goal of
Equation (6) is to compute $m_{1}\wedge m_{2}$. Let us look into it more in
detail. Then call $LweSub$ to calculate $c^{{}^{\prime}}$ by subtraction.
Similarly in Equation (7), after mapping message element to torus element,
$LweAdd$ is called to calculate $c^{{}^{\prime}}$. As shown in Equation (10),
$LweSub$ and $LweMul$ should be called. $NOT$ is the simplest operation which
only needs subtraction of $LweSub$. Luckily, $ModToT$ and $LweSub$ have been
provided in the library of MKTFHE. We have to design and implement two
algorithms $LweAdd$ and $LweMul$ by ourselves, to achieve addition and
multiplication of two TLWE ciphertexts, independently.

###  4.2 Two important components used in bootstrapped gates

This subsection illustrates two components used in bootstrapped gates,
$LweAdd$ and $LweMul$, which are missed in [7]. The idea of $LweAdd$ is to
execute addition on the $n+1$-dimension vectors $c_{1}$ and $c_{2}$ element by
element. Remember to refresh the current variance by adding the variance of
ciphertext $c_{1}$ and the variance of ciphertext $c_{2}$ to get the variance
of ciphertext $c$, as shown in Algorithm 1. As shown in Equation (4), $LweMul$
is only required in $XOR$. Such multiplication is dot multiplication, which
can be expressed by addend $\\{a_{1}\\}_{i,j}$ plus $k$-time
$\\{a_{2}\\}_{i,j}$. Remember to refresh the current variance by adding the
variance of ciphertext $c_{1}$ and the variance of ciphertext $c_{2}$ times
$k^{2}$ to get the variance of ciphertext $c$, as shown in Algorithm 2.

Algorithm 1 Multi-key TLWE ciphertexts addition ($LweAdd$)
---
Input: | MKTFHE parameter set, number of the participants $p$ and two MKTLWE ciphertexts
|  $c_{1}=(a_{1},b_{1})$ and $c_{2}=(a_{2},b_{2})$
Output: | A MKTLWE ciphertext $c=(a,b)=c_{1}+c_{2}$
1. | Extract the dimension of the lattice $n$ in parameter set.
2. | For $i=1$ to $p$ do
3. |  For $j=1$ to $n$ do
4. |  Compute $\\{a\\}_{i,j}=\\{a_{1}\\}_{i,j}+\\{a_{2}\\}_{i,j}$
5. |  End for
6. | End for
7. | Compute $b=b_{1}+b_{2}$
8. | Refresh the current variance

Algorithm 2 Multi-key TLWE ciphertexts multiplication ($LweMul$)
---
Input: | MKTFHE parameter set, number of the participants $p$, the coefficient $k$ and two MKTLWE
| ciphertexts $c_{1}=(a_{1},b_{1})$ and $c_{2}=(a_{2},b_{2})$
Output: | A MKTLWE ciphertext $c=(a,b)=c_{1}+k\cdot c_{2}$
1. | Extract the dimension of the lattice $n$ in parameter set.
2. | For $i=1$ to $p$ do
3. |  For $j=1$ to $n$ do
4. |  Compute $\\{a\\}_{i,j}=\\{a_{1}\\}_{i,j}+k\cdot\\{a_{2}\\}_{i,j}$
5. |  End for
6. | End for
7. | Compute $b=b_{1}+k\cdot b_{2}$
8. | Refresh the current variance

However, two components used in bootstrapped gates above, $LweAdd$ and
$LweMul$, which are missed in [7]. The idea of $LweAdd$ is to execute addition
on the $n+1$-dimension vectors $c_{1}$ and $c_{2}$ element by element.
Remember to refresh the current variance by adding the variance of ciphertext
$c_{1}$ and the variance of ciphertext $c_{2}$ to get the variance of
ciphertext $c$, as shown in Algorithm 1. As shown in Equation (4), $LweMul$ is
only required in $XOR$. Such multiplication is dot multiplication, which can
be expressed by addend $\\{a_{1}\\}_{i,j}$ plus $k$-time $\\{a_{2}\\}_{i,j}$.
Remember to refresh the current variance by adding the variance of ciphertext
$c_{1}$ and the variance of ciphertext $c_{2}$ times $k^{2}$ to get the
variance of ciphertext $c$, as shown in Algorithm 2.

We summarize the operations between two LWE ciphertexts used in each basic
binary gate in Table 1.

Table 1: Operations used in each basic gate in our scheme Gate | $BS$ | $Add$ | $Sub$ | $Mul$ | $ModToT$
---|---|---|---|---|---
AND | $\surd$ | $\surd$ |  |  | $\surd$
OR | $\surd$ | $\surd$ |  |  | $\surd$
NOT |  |  | $\surd$ |  | $\surd$
NAND | $\surd$ |  | $\surd$ |  | $\surd$
NOR | $\surd$ |  | $\surd$ |  | $\surd$
XOR | $\surd$ |  | $\surd$ | $\surd$ |
XNOR | $\surd$ |  | $\surd$ | $\surd$ | $\surd$

###  4.3  $k$-bit homomorphic adder

After achieving various basic gate operations in subsection 4.1, it is still a
gap between such gates and complicated machine learning functions or private
computing directly. In this subsection, we construct a $k$-bit complement
array integer adder based on the MKTFHE scheme. We choose complement
representation instead of true form or inverse code, for the complement can
compute both positive and negative numbers in the same way. A $k$-bit
complement array integer adder can be constructed in the following four steps.
The main idea was described in subsection 3.3. The only difference is that
bootstrapping $BS$ is necessary after each computation to decrease noise. The
$k$-bit complement array integer adder is shown in Algorithm 3.

Algorithm 3 $k$-bit complement array integer adder
---
Input: | MKTFHE parameter set, two ciphertexts $c_{1}=\mathsf{MKTF}$\- $\mathsf{HE.SymEnc}(\mu_{1})$ and
|  $c_{2}=\mathsf{MKTFHE.SymEnc}(\mu_{2})$ and the public keys of all
participants
Output: | Ciphertext $c=\mathsf{MKTFHE.SymEnc}(\mu_{1}+\mu_{2})$
1. | Set the first carry $c_{c}[0]=0$
2. | Encrypt addends and carry with public keys
3. | For $i=1$ to $k-1$ do
4. |  Compute the $out[i]$ and $c_{c}[i]$
5. | End for

  1. 1.

Construct a homomorphic semi-adder with two bootstrapped XOR gates: Input two
ciphertexts $c_{a}$ and $c_{b}$ and the encrypted carry $c_{c}$, return
$out=BS(c_{c}\bigoplus(BS(c_{a}\bigoplus c_{b})))$, where $BS$ stands for
bootstrapping process.

  2. 2.

Construct a homomorphic carrier with three bootstrapped AND gates and two
bootstrapped OR gates: Input two ciphertexts $c_{a}$ and $c_{b}$ and the
encrypted carry $c_{c}$, return $cout=BS(BS(BS(c_{a}\wedge c_{b})\vee
BS(c_{a}\wedge c_{c}))\vee BS(c_{b}\wedge c_{c}))$.

  3. 3.

Construct a homomorphic $1$-bit adder by merging a homomorphic semi-adder and
a homomorphic carrier: Input two ciphertexts $c_{a}$ and $c_{b}$ and the
encrypted carry $c_{c}$, return $out=BS(c_{c}\bigoplus(BS(c_{a}\bigoplus
c_{b})))$ and $cout=BS(BS(BS(c_{a}\wedge c_{b})\vee BS(c_{a}\wedge c_{c}))\vee
BS(c_{b}\wedge c_{c}))$.

  4. 4.

Construct a homomorphic $k$-bit complement array integer adder with k
homomorphic $1$-bit adder: Input two ciphertexts $c_{a}$ and $c_{b}$, set the
encrypted carry $c_{c}[0]=0$, compute
$out[i]=BS(c_{c}[i]\bigoplus(BS(c_{a}[i]\bigoplus c_{b}[i])))$ and
$c_{c}[i+1]=BS(BS(BS(c_{a}[i]\wedge c_{b}[i])\vee BS(c_{a}[i]\wedge
c_{c}[i]))\vee BS(c_{b}[i]\wedge c_{c}[i]))$, where $i=\\{0,1,...,$k$-1\\}$.
Return the result $out$.

In $k$-bit complement array integer adder, the output carry of the last
$1$-bit adder is abandoned because of the capture of the complement.

Figure 4: The construction of $k$-bit homomorphic complement array integer
adder. GB stands for gate-bootstrapping.

###  4.4  $k$-bit homomorphic multiplier

We also construct a $k$-bit complement array integer multiplier based on the
MKTFHE scheme. A $k$-bit complement array integer multiplier can be
constructed in the following steps. The main idea was described in subsection
3.3. The only difference is that bootstrapping $BS$ is necessary after each
computation to decrease noise. A $k$-bit complement array integer multiplier
is shown in Algorithm 4.

Algorithm 4 $k$-bit complement array integer multiplier
---
Input: | MKTFHE parameter set, two ciphertexts $c_{1}=\mathsf{MKTF}$-$\mathsf{HE.SymEnc}(\mu_{1})$ and
|  $c_{2}=\mathsf{MKTFHE.SymEnc}(\mu_{2})$ and the public keys of all
participants.
Output: | A ciphertext $c=\mathsf{MKTFHE.SymEnc}(\mu_{1}\times\mu_{2})$
1. | Set the first carry $c_{c}[0]=0$
2. | Encrypt numbers and carry with public keys
3. | For $i=1$ to $k-2$ do
4. |  For $j=0$ to $k-2$ do
5. |  Compute the AND gate.
6. |  If $j=0$, compute the 2-homadder
7. |  Elseif $i+j\geq k-1$, compute the $1$-homadder
8. |  Else compute the 0-homadder
9. |  End for
10. |  Get the $(j+k)$th bit of the result
11. | End for
12. | For $j=0$ to $k-1$ do
13. |  Compute the $j$th bit of the output by 2-homadder.
14. | End for

  1. 1.

Prepare three kinds of homomorphic $1$-bit adder:

    1. (a)

0-homadder: The same as the homomorphic $1$-bit adder above.

    2. (b)

$1$-homadder: Input two ciphertexts $c_{a}$ and $c_{b}$ and the encrypted
carry $c_{c}$, return
$out=\overline{BS(c_{c}\bigoplus(BS(\overline{c_{a}}\bigoplus}$
$\overline{c_{b})))}$ and $c_{out}=BS(BS(BS(\overline{c_{a}}\wedge c_{b})\vee
BS(\overline{c_{a}}\wedge c_{c}))\vee BS(c_{b}\wedge c_{c}))$.

    3. (c)

2-homadder: Input two ciphertexts $c_{a}$ and $c_{b}$ and the encrypted carry
$c_{c}$, return $out=BS(c_{c}\bigoplus(BS(\overline{c_{a}}\bigoplus$
$\overline{c_{b}})))$ and $c_{out}=\overline{BS(BS(BS(\overline{c_{a}}\wedge}$
$\overline{\overline{c_{b}})\vee BS(\overline{c_{a}}\wedge c_{c}))}$
$\overline{\vee BS(\overline{c_{b}}\wedge c_{c}))}$.

  2. 2.

Construct a $k$-bit homomorphic complement array integer multiplier according
to the rules in Section 3.4 but use homomorphic adder instead of the previous
adder.

We use Algorithm 4 to evaluate $k$-bit complement array integer multiplier
homomorphically, where $k$ is the number of bits of the addends. In $k$-bit
complement array integer multiplier, the output carry of the last $1$-bit
adder is considered as the highest bit of the result.

###  4.5  $k$-bit homomorphic subtractor and divider

We construct a $k$-bit complement array integer subtractor by using a $k$-bit
complement array integer adder and a $k$-bit complement array integer
multiplier. A $k$-bit complement array integer subtractor is shown in
Algorithm 5.

We also construct a $k$-bit complement array integer divider based on the
MKTFHE scheme. The main idea was described in subsection 3.3. The only
difference is that bootstrapping $BS$ is necessary after each computation to
decrease noise. A $k$-bit complement array integer divider can be constructed
in the following steps:

  1. 1.

Prepare the homomorphic CAS for the divider. We use three homomorphic
bootstrapped XOR gates, two homomorphic bootstrapped OR gates and two
homomorphic bootstrapped AND gates to construct the homomorphic CAS. The
homomorphic CAS takes ciphertexts $c_{a}[i]$, $c_{b}[i]$, $c_{c}[i]$ and
$c_{p}$ as input, and output $c_{out}[i]$ and the carrier $c_{c}[i+1]$
satisfying $c_{out}[i]=BS(c_{a}[i]\bigoplus BS(BS(c_{b}[i]\bigoplus
c_{p})\bigoplus c_{c}[i]))$ and $c_{c}[i+1]=BS(BS(BS(c_{a}[i]\vee
c_{c}[i])\wedge BS(c_{b}[i]\bigoplus c_{p}))\vee BS(c_{a}[i]\wedge
c_{c}[i]))$. The homomorphic CAS will perform addition if the message of
$c_{p}$ is $0$ and perform subtraction if the message of $c_{p}$ is $1$.

  2. 2.

Then design the homomorphic absolute value array divider with the homomorphic
CAS above. It takes the $2k$-bit dividend $c_{a}$ and the $k$-bit divisor
$c_{b}$ as input, and output the quotient $c_{q}$ and the remainder $c_{r}$.

  3. 3.

Homomorphic XOR gates are used to decide the sign bit of the ciphertext of
quotient to realize the homomorphic division of both positive and negative
integers.

  4. 4.

A homomorphic compensation device is also designed to realize the complement
division.

  5. 5.

Construct the $k$-bit homomorphic complement array divider with a homomorphic
absolute value array divider, a homomorphic XOR gate and two homomorphic
compensation device. The $k$-bit homomorphic complement array divider takes
the ciphertexts of $2k$-bit complement dividend $a_{in}$ and the $k$-bit
complement divisor $b_{in}$ as input, and output the ciphertexts of complement
quotient $q$ and the complement remainder $r$.

Algorithm 5 $k$-bit complement array integer subtractor
---
Input: | MKTFHE parameter set, two ciphertexts $c_{1}=\mathsf{MKTF}$-$\mathsf{HE.SymEnc}(\mu_{1})$ and
|  $c_{2}=\mathsf{MKTFHE.SymEnc}(\mu_{2})$ and the public keys of all
participants.
Output: | A ciphertext $c=\mathsf{MKTFHE.SymEnc}(\mu_{1}-\mu_{2})$
1. | Change every bit of $c_{2}$ and add $1$ to calculate $-c_{2}$
2. | Set the first carry $c_{c}[0]=0$
3. | Encrypt addends and carry with public keys
4. | For $i=1$ to $k-1$ do
5. |  Compute the $out[i]$ and $c_{c}[i]$
6. | End for

Algorithm 6 $k$-bit complement array integer divider
---
Input: | MKTFHE parameter set, two ciphertexts $c_{1}=\mathsf{MKTF}$-$\mathsf{HE.SymEnc}(\mu_{1})$ and
|  $c_{2}=\mathsf{MKTFHE.SymEnc}(\mu_{2})$ and the public keys of all
participants.
Output: | Two ciphertext $q=\mathsf{MKTFHE.SymEnc}(\mu_{1}/\mu_{2})$ and $q=\mathsf{MKTFHE.SymEnc}(\mu_{1}mod\mu_{2})$
1. | Encrypt the $p_{0}$ with public keys
2. | Change the inputs into complement format with compensation device.
3. | Calculate the first later of divider.
4. | For $i=1$ to $k-1$ do
5. |  For $j=0$ to $k-1$ do
5. |  Compute the CAS unit.
6. |  End for
7. | End for

We use Algorithm 6 to evaluate $k$-bit complement array integer divider
homomorphically, where $k$ is the number of bits of the divisor.

In summary, this section first expands NAND gates to other six basic gates,
AND, OR, NOT, NOR, XOR and XNOR. In order to compute these gates, two
important components, $LweAdd$ and $LweMul$ are then proposed. At last, we
construct the adder, subtractor, multiplier and divider. With these
homomorphic operators, we can evaluate arbitrary polynomials without
decryption.

###  4.6 Multi-key homomorphic distributed linear regression scheme

We first implement a multi-key fully homomorphic encryption scheme in which
the secret key is distributed among the parties, while the corresponding
collective public key $pk$ is known to all of them. Thus, each party can
independently compute on ciphertexts encrypted under $pk$ but all parties have
to collaborate to decrypt a ciphertext. This enables the participants to train
a collectively encrypted model, that cannot be decrypted as long as one
participant is honest and refuses to participate in the decryption. Then we
implement two methods: calculation formula method and GD method to train and
evaluate the linear regression model.

Calculation formula method: The advantage of this method is that the global
optimal solution can be obtained without multiple iterative operation which
can reduce the computational cost of time and space. We first estimate the
size of the input data, and then select homomorphic operators with appropriate
bits to calculate the above formula. Since our operators are integer
operators, we choose to directly round off the input data under the condition
of ensuring accuracy. We use algorithm 7 to train the linear regression model
under multi-key homomorphic encryption.

GD method: Considering that our homomorphic operators only support integer
calculation, and the learning rate $\alpha$ is usually a floating-point number
less than 1, we rewrite the GD iterative formula in Subsection 3.4 to zoom the
$\alpha$ to integers. We bring the linear regression model into the origin GD
iterative formula, the zooming multiple is $n$, and the practical integer GD
iterative formula is below:

$b^{\prime}=b\cdot n-\alpha\cdot
n\cdot\frac{2}{m}\cdot\sum_{i=1}^{m}[y_{i}\cdot n-(\frac{\omega}{n}\cdot
x_{i}-\frac{b}{n})]$ (12)

$\omega^{\prime}=\omega\cdot n-\alpha\cdot
n\cdot\frac{2}{m}\cdot\sum_{i=1}^{m}[y_{i}\cdot n-(\frac{\omega}{n}\cdot
x_{i}-\frac{b}{n})]$ (13)

$loss=\frac{1}{m}\cdot\sum_{i=1}^{m}[y_{i}-(\omega\cdot x_{i}+b)/n]^{2}$ (14)

After we get the rewritten iterative formula, we estimate the size of the
input data to select the zooming multiple, iterative times and homomorphic
operators with appropriate numbers of bits.

We use Algorithm 8 to train the linear regression model under multi-key
homomorphic encryption.

Algorithm 7 Calculation formula method in multi-key fully homomorphic linear
regression
---
Input: | MKTFHE parameter set, two ciphertexts $c_{x_{n}}=\mathsf{MKTF}$\- $\mathsf{HE.SymEnc}(x_{n})$ and
|  $c_{y_{n}}=\mathsf{MKTFHE.SymEnc}(y_{n})$ and the public keys of all
participants
Output: | Two ciphertext $c_{\omega}=\mathsf{MKTFHE.SymEnc}(\omega)$ and $c_{b}=\mathsf{MKTFHE.SymEnc}(b)$
1. | Reconstruct the public keys to bootstrapping key
2. | Extend the sing-key ciphertext to multi-key ciphertext
3. | Calculate the average value $c_{\bar{x}}$ of $c_{x}$, following the Equation (1)
4. | Calculate the $c_{\omega}$ by $c_{\bar{x}}$, following the Equation (2)
5. | Calculate the $c_{b}$ by $c_{\omega}$, following the Equation (3)

Algorithm 8 GD method in multi-key fully homomorphic linear regression
---
Input: | MKTFHE parameter set, zooming multiple $n$, iterative times $d$, two sets of ciphertexts
|  $c_{x_{n}}=\mathsf{MKTF}$\- $\mathsf{HE.SymEnc}(x_{n})$ and
$c_{y_{n}}=\mathsf{MKTFHE.SymEnc}(y_{n})$ and the public keys of all
| participants
Output: | Two ciphertext $c_{\omega}=\mathsf{MKTFHE.SymEnc}(\omega)$ and $c_{b}=\mathsf{MKTFHE.SymEnc}(b)$
1. | Reconstruct the public keys to bootstrapping key
2. | Extend the sing-key ciphertext to multi-key ciphertext
3. | Given iteration initial value to $\omega$ and $b$
4. | For $i=1$ to $d$ do
5. |  Calculate the GD iterative equation (13) (14)
6. |  Set $\omega=\omega^{\prime}$ and $b=b^{\prime}$
7. | End for

##  5 Implementation and experiments

The test environment is Ubuntu 18.04 operation system, with Intel Xeon Gold
5220 CPU and 46512MiB memory.

###  5.1 Experiments on homomorphic basic gates

The MKTFHE scheme provides only homomorphic bootstrapped NAND gates. Our
scheme constructs other basic homomorphic bootstrapped gates in the same way
as the bootstrapped NAND gate in MKTFHE. We compared our basic bootstrapped
gates with gates built by the NAND gates in the same test environment.
Experiment shows that our basic bootstrapped gates are as efficient as the
bootstrapped NAND gate in the MKTFHE scheme, which is much more efficient than
the gates with the same function constructed by only bootstrapped NAND gates.

Table 2: Expression of our basic bootstrapped gates and the average time of each gate in both our scheme and naïve scheme. Gate | Naïve scheme(s) | Our scheme(s)
---|---|---
AND | 0.238671 | 0.238008
OR | 0.584701 | 0.23558
NOT | 0.236541 | 2e-06
NAND | 0.235327 | 0.239638
NOR | 0.470561 | 0.236153
XOR | 0.71246 | 0.2354
XNOR | 0.710777 | 0.23588

The result in Table 2 shows that our scheme is more efficient than the naïve
scheme, which connects the bootstrapped NAND gates to realize the operations.
We have improved in efficiency, especially in bootstrapped OR gate,
bootstrapped NOR gate, NOT gate, bootstrapped XOR gate and bootstrapped XNOR
gate: the time cost by a single bootstrapped XOR gate reduced by about 67%
while the NOT gate is evaluated almost instantaneously for no bootstrapping is
needed.

###  5.2  Experiments on $k$-bit homomorphic adder

We construct a $k$-bit homomorphic complement integer adder based on the basic
bootstrapped gates we designed above. We create two participants for our
experiment. Each of them creates its addend and secret key individually. For
simple, we set $k=1,2,...,8$.

The server input multi-key parameters, the ciphertexts encrypted under
different keys by two participants, and the public keys of two participants.
Then the server extends the ciphertext, evaluates the $k$-bit adder
homomorphically, and returns the ciphertext to the participants. The
participants decrypt the ciphertext and get the result.

According to the structure of $k$-bit homomorphic adder, while computing the
addition two $k$-bit addends, $k$ $1$-bit adder is required. As shown in
subsection 4.2, it requires $5$ bootstrapped gates to construct a $1$-bit
adder. That is to say, $5k$ basic bootstrapped gates are needed to construct a
$k$-bit homomorphic adder. Obviously, the cost of the homomorphic evaluation
grows linearly with the number of bootstrapped gates. So, the cost of the
adder grows linearly with the bits of input numbers.

Table 3: The average time of $k$-bit homomorphic adder. Bits | Time(s) | Num of gates | Time/gate(s)
---|---|---|---
$1$-bit | 1.19515 | 5 | 0.23903
$2$-bit | 2.39738 | 10 | 0.239738
$3$-bit | 3.60485 | 15 | 0.240323
$4$-bit | 4.80614 | 20 | 0.240307
$5$-bit | 6.01884 | 25 | 0.240754
$6$-bit | 7.25483 | 30 | 0.241828
$7$-bit | 8.43608 | 35 | 0.241031
$8$-bit | 9.69468 | 40 | 0.242367
Figure 5: The cost of $k$-bit homomorphic complement array integer adder

From Table 3 and Figure 5, we can learn that the number of binary bootstrapped
gates grows linearly with the number of the bits of addends. As evaluation is
the most expensive step in our scheme, we can say that the cost of the $k$-bit
homomorphic complement adder grows linearly with the bits of input numbers.

###  5.3  Experiments on $k$-bit homomorphic multiplier

To construct a $k$-bit homomorphic multiplier, we construct 3 kinds of
homomorphic adder: 0-homadder, 1-homadder and 2-homadder. We also construct a
$k$-bit homomorphic complement integer multiplier based on three kinds of
adders and the basic bootstrapped gates we designed above. The structure of 3
kinds of homomorphic adder in section 4.3 shows that it requires $7$ basic
bootstrapped gates to construct a homomorphic adder: $5$ gates in $1$-bit
adder and two AND gates to compute the input. It can be learned from the
structure of the $k$-bit homomorphic multiplier in section 3.3 that it
requires $k\times(k-1)$ homomorphic adders to construct a $k$-bit homomorphic
multiplier. That is to say, it requires $7\times k\times(k-1)$ basic
bootstrapped gates for a $k$-bit homomorphic multiplier. So, the cost of the
multiplier grows quadratically with the bits of input numbers. Note that there
is parallel optimization in the multiplier, so the average time consumed by a
gate is much less than that consumed by a single gate.

We create two participants for our experiment, each of them creates its
complement number and secret key individually. For simplicity, we set
$k=2,3,...,8$. The server input multi-key parameters, the ciphertexts
encrypted under different keys by two participants, and the public keys of two
participants. Then the server extends the ciphertext, evaluates the $k$-bit
multiplier homomorphically, and returns the ciphertext to the participants.
The participants decrypt the ciphertext and get the result.

Table 4: The average time of each step of experiments on $k$-bit homomorphic adder, as well as the number of gates used in multiplier. bit of mul | Time(s) | Gates | Time/gate(s)
---|---|---|---
2-bit mul | 1.14368 | 14 | 0.081691429
3-bit mul | 2.66844 | 42 | 0.063534286
4-bit mul | 4.83101 | 84 | 0.057512024
5-bit mul | 7.43843 | 140 | 0.053131643
6-bit mul | 10.762 | 210 | 0.051247619
7-bit mul | 14.6522 | 294 | 0.049837415
8-bit mul | 19.1052 | 392 | 0.048737755
Figure 6: The cost of $k$-bit homomorphic complement array integer multiplier

From Table 4 and Figure 6, we can learn that the number of binary bootstrapped
gates grows quadratically with the number of the bits of addends. Similarly,
we can say that the cost of the $k$-bit homomorphic complement multiplier
grows quadratically with the bits of input numbers.

###  5.4  Experiments on $k$-bithomomorphic subtractor and divider

We construct a $k$-bit homomorphic complement integer subtractor based on the
basic bootstrapped gates we designed above. We create two participants for our
experiment. Each of them creates its addend and secret key individually. For
simple, we set $k=1,2,...,8$.

According to the structure of $k$-bit homomorphic subtractor, while computing
the $k$-bit subtraction, $k$ $1$-bit adder is required. As shown in subsection
4.2, it requires $5$ bootstrapped gates to construct a $1$-bit adder. At the
same time, a bootstrapped XOR gate is used to compute the negation. That is to
say, $6k$ basic bootstrapped gates are needed to construct a $k$-bit
homomorphic adder. Obviously, the cost of the homomorphic evaluation grows
linearly with the number of bootstrapped gates. So, the cost of the subtractor
grows linearly with the bits of input numbers.

Table 5: The average time of experiments on $k$-bit homomorphic subtractor bit of sub | Time(s) | Gates | Time/gate(s)
---|---|---|---
1 | 1.99209 | 6 | 0.332015
2 | 2.90969 | 12 | 0.242474167
3 | 4.45103 | 18 | 0.247279444
4 | 5.84264 | 24 | 0.243443333
5 | 7.42889 | 30 | 0.247629667
6 | 8.93639 | 36 | 0.248233056
7 | 10.35 | 42 | 0.246428571
8 | 11.9013 | 48 | 0.24794375
Figure 7: The cost of $k$-bit homomorphic complement array integer subtractor

Experiments show that the cost of $k$-bit homomorphic subtractor grows
linearly with the bit of minuend.

We also construct a $k$-bit homomorphic complement integer divider in the
similar way. We create two participants for our experiment. Each of them
creates its addend and secret key individually. For simple, we set
$k=1,2,...,8$.

According to the structure of $k$-bit homomorphic divider, while computing the
$k$-bit division, a homomorphic absolute value array divider, a homomorphic
XOR gate and two homomorphic compensation device are required. There are $7$
homomorphic gates in each hom-CAS while there are $k^{2}$ hom-CASs in each
homomorphic absolute value array divider. It takes $2k$ homomorphic gates to
construct a compensation device. As a result, $7k^{2}+2k+1$ basic bootstrapped
gates are needed to construct a $k$-bit homomorphic divider. Note that there
are more optimization in the divider than in the multiplier, and the CAS
almost run in parallel, so a $k$-bit divider is faster than a $k$-bit
multiplier, and the cost of the divider grows almost linearly with the number
of layers.

Table 6: The average time of experiments on $k$-bit homomorphic divider bit of div | Time(s) | Layers | Time/layer(s)
---|---|---|---
1 | 1.483 | 1 | 1.483
2 | 4.09858 | 2 | 2.04929
3 | 7.08719 | 3 | 2.362396667
4 | 10.0022 | 4 | 2.50055
5 | 13.1937 | 5 | 2.63874
6 | 15.965 | 6 | 2.660833333
7 | 19.1003 | 7 | 2.728614286
8 | 22.2938 | 8 | 2.786725
Figure 8: The cost of $k$-bit homomorphic complement array integer divider

###  5.5 Experiments on multi-key homomorphic distributed linear regression

We implement a multi-key fully homomorphic encryption scheme in which the
secret key is distributed among the parties, while the corresponding
collective public key $pk$ is known to all of them. Thus, each party can
independently compute on ciphertexts encrypted under $pk$ but all parties have
to collaborate to decrypt a ciphertext. This enables the participants to train
a collectively encrypted model, that cannot be decrypted as long as one
participant is honest and refuses to participate in the decryption.

Our multi-key fully homomorphic linear regression scheme is composed by
multiple parties which can be divided by three types of entities:
participants, cloud server, and decryption party, we take two participants as
examples. the whole scheme are shown in the Figure 9, the steps are as
follows:

  1. 1.

Participants have the need to outsource computing, and each participant offers
its own part of data for model training. During the step of data encryption,
all participants call $\mathsf{MKTFHE.KeyGen}(mkparams)$ to generate their
symmetric key, asymmetric key bootstrapping key and key-switching key
independently, then call $\mathsf{MKTFHE.SymEnc}(\mu)$ to encrypt the data.
And finally upload the ciphertext of input data and public key to the cloud
server.

  2. 2.

Cloud server is usually composed of one or more high-performance servers, and
doesn’t have their own data. After receiving the ciphertext data from
participants, the cloud server will use the participants’ key to generate a
new bootstrapping key and extend the single-key ciphertext of the input data
to multi-key ciphertext. Then cloud server uses our designed homomorphic
operators to train the linear regression model. After the training is
finished, the result of the model will be sent to the decryption party.

  3. 3.

When the decryption party receive the ciphertext of the trained model, they
will unite all the participants to call
$\mathsf{MKTFHE.SymDec}(c,\\{sk_{i}\\})$ to decrypt the result.

Figure 9: Multi-key fully homomorphic linear regression scheme

We show our setting of the parameters in the Table 7, following the notation
of MKTFHE library. The achievement estimated security level of our scheme is
110-bit while the dimension of the TLWE problem is $k=1$.

Table 7: Parameter sets of multi-key homomorphic linear regression

LWE-$n$ | LWE-$\alpha$ | LWE-$B^{\prime}$ | LWE-$d^{\prime}$ | RLWE-$N$ | RLWE-$\beta$ | RLWE-$B$ | RLWE-$d$
---|---|---|---|---|---|---|---
560 | $3.05\times 10^{-5}$ | $2^{2}$ | 8 | 1024 | $3.72\times 10^{-9}$ | $2^{9}$ | 3

We train a multi-key fully homomorphic linear regression model in two methods
by using our proposed homomorphic operators. The input data is composed of
several sets of linear data and some random noise. Considering the size of the
data, we chose 8-bit homomorphic operators to train the model in formula
method and 16-bit homomorphic operators in GD method. In addition, the
learning rate and the zooming multiple of the GD method is 0.001 and 10,000.
We make several experiments and find that the model will converge by
performing up to 10 iterations. We also create $m$ participants in our
experiment. Each of them generates their own secret key to encrypt their own
data. For simple, we set $m=2,4,8$, and the $3,4,5,6$ columns record the time
of each step.

Table 8: The result of experiments on multi-key fully homomorphic linear
regression

method | Num of Parties | KeyGen(s) | Ciphertext extension(s) | Training(s) | Evaluation(s)
---|---|---|---|---|---
Formular | $k=2$ | $1.984$ | $0.0008$ | $227.467$ | $37.131$
Formular | $k=4$ | $4.002$ | $0.0016$ | $445.593$ | $68.968$
Formular | $k=8$ | $8.916$ | $0.0035$ | $869.148$ | $138.259$
GD | $k=2$ | $2.025$ | $0.0008$ |  $1133.33$/iter | $190.644$
GD | $k=4$ | $4.008$ | $0.0015$ |  $2092.17$/iter | $372.353$

In Table 8, the result shows that: the running time including key
reconstruction time, extension ciphertext time, training time and evaluation
time all grow linearly with the numbers of participants; due to the large
zooming multiple of the GD method, the homomorphic operators with larger bits
are selected. At the same time, the GD method needs multiple iterations, and
the calculation process is more complex, so the running time including
training time and evaluation time is much more than that of the formula
method.

Note that there is more optimization in our linear regression scheme such as
running in parallel. Besides, with our homomorphic operators, we can implement
more complex multi-key fully homomorphic machine learning scheme.

##  6 Conclusion and Discussion

###  6.1 Conclusion

The scheme proposes a series of $k$-bit homomorphic complement operators based
on MKTFHE, which narrows the gap between the original NAND gate and
complicated machine learning functions. Experiment shows that the cost of the
adder and the subtractor grows linearly with the bits of input numbers and the
cost of the multiplier grows quadratically. Meanwhile, the cost of the divider
grows almost linearly with the number of layers in it.

To narrow the gap between multi-key bootstrapped NAND gates and multi-key
homomorphic mathematical operators, we construct other basic bootstrapped
gates (AND, OR, NOT, NOR, XOR, and XNOR) in the same way as the bootstrapped
NAND gate in MKTFHE and with the same efficiency as the NAND gate. Experiment
shows that constructing basic binary gates in this way is much more efficient
than building them by directly joining the bootstrapped NAND gates together,
especially when constructing XOR gate and NOT gate. Then we construct a
$k$-bit complement adder, a $k$-bit complement subtractor, a $k$-bit
complement multiplier and a $k$-bit complement divider based on our basic
binary bootstrapped gates. Finally, we train distributed linear regression
model by utilizing our proposed multi-key homomorphic operators. The operators
we designed can be directly used to achieve privacy-preserving distributed
machine learning schemes in distributed communication system.

###  6.2 Discussion

The plaintext space of the MKTFHE scheme is $\\{0,1\\}$, which is smaller than
the plaintext space in other homomorphic encryption schemes like BGV [34] or
CKKS [35, 36]. Extending the message space on torus may be helpful. It has
been experimented that integers or fixed-point numbers can also be mapped to a
ring on the torus [37, 38] as well as performing homomorphic evaluations. If
the message space of MKTFHE is magnified to integers or fixed-point numbers,
it may reduce the overhead in time and space in MKTFHE.

Common reference strings (CRS) are needed in this scheme for all the
participants, and the computing server has to know the multi-key parameters in
advance. If the CRS can be removed in MKTFHE, the scheme will no longer need a
trusted third party to generate the CRS, making the scheme more secure.

We described only multi-key homomorphic linear regression model in our scheme,
but more machine learning schemes (such as logistic regression, SVM, etc) can
be evaluated in the same way, for the basic mathematical operators are
provided in our scheme.

##  7 Acknowledgement

This work is supported by National Natural Science Foundation of China (No.
61872109), Shenzhen Basic Research Project, China (No. JCYJ20180507183624136),
Shenzhen Basic Research Project, China (No. JCYJ20200109113405927) and
National Science and Technology Major Project Carried on by Shenzhen, China
(No. CJGJZD20200617103000001).

## References

  * [1] J. Verbraeken, M. Wolting, J. Katzy, J. Kloppenburg, T. Verbelen, J. S. Rellermeyer, A survey on distributed machine learning, ACM Computing Surveys (CSUR) 53 (2) (2020) 1–33.
  * [2] Q. Jia, L. Guo, Z. Jin, Y. Fang, Preserving model privacy for machine learning in distributed systems, IEEE Transactions on Parallel and Distributed Systems 29 (8) (2018) 1808–1822.  doi:10.1109/TPDS.2018.2809624.
  * [3] J. Kolberg, P. Drozdowski, M. Gomez-Barrero, C. Rathgeb, C. Busch, Efficiency analysis of post-quantum-secure face template protection schemes based on homomorphic encryption, in: BIOSIG 2020 - Proceedings of the 19th International Conference of the Biometrics Special Interest Group, Vol. P-306 of LNI, 2020, pp. 175–182.
  * [4] J. L. H. Crawford, C. Gentry, S. Halevi, D. Platt, V. Shoup, Doing real work with FHE: the case of logistic regression, in: Proceedings of the 6th Workshop on Encrypted Computing & Applied Homomorphic Cryptography, WAHC@CCS 2018, 2018, pp. 1–12.
  * [5] C. Gentry, Fully homomorphic encryption using ideal lattices, in: Proceedings of the 41st Annual ACM Symposium on Theory of Computing, STOC 2009, 2009, pp. 169–178.
  * [6] A. López-Alt, E. Tromer, V. Vaikuntanathan, On-the-fly multiparty computation on the cloud via multikey fully homomorphic encryption, in: Proceedings of the 44th Symposium on Theory of Computing Conference, STOC 2012, 2012, pp. 1219–1234.
  * [7] H. Chen, I. Chillotti, Y. Song, Multi-key homomorphic encryption from TFHE, in: Advances in Cryptology - ASIACRYPT 2019 - 25th International Conference on the Theory and Application of Cryptology and Information Security, Vol. 11922 of Lecture Notes in Computer Science, 2019, pp. 446–472.
  * [8] K. Hwang, Global and modular two’s complement cellular array multipliers, IEEE Trans. on Computers 28 (4) (1979) 300–306.
  * [9] D. Froelicher, J. R. Troncoso-Pastoriza, A. Pyrgelis, S. Sav, J. S. Sousa, J.-P. Bossuat, J.-P. Hubaux, Scalable privacy-preserving distributed learning, arXiv preprint arXiv:2005.09532.
  * [10] Y. Aono, T. Hayashi, L. Trieu Phong, L. Wang, Scalable and secure logistic regression via homomorphic encryption, in: Proceedings of the Sixth ACM Conference on Data and Application Security and Privacy, 2016, pp. 142–144.
  * [11] J. L. Crawford, C. Gentry, S. Halevi, D. Platt, V. Shoup, Doing real work with fhe: The case of logistic regression, in: Proceedings of the 6th Workshop on Encrypted Computing & Applied Homomorphic Cryptography, 2018, pp. 1–12.
  * [12] A. Kim, Y. Song, M. Kim, K. Lee, J. H. Cheon, Logistic regression model training based on the approximate homomorphic encryption, BMC medical genomics 11 (4) (2018) 23–31.
  * [13] M. Kim, Y. Song, S. Wang, Y. Xia, X. Jiang, et al., Secure logistic regression based on homomorphic encryption: Design and evaluation, JMIR medical informatics 6 (2) (2018) e8805.
  * [14] T. Graepel, K. E. Lauter, M. Naehrig, ML confidential: Machine learning on encrypted data, in: Information Security and Cryptology - ICISC 2012, Vol. 7839, 2012, pp. 1–21.
  * [15] R. Gilad-Bachrach, N. Dowlin, K. Laine, K. E. Lauter, M. Naehrig, J. Wernsing, Cryptonets: Applying neural networks to encrypted data with high throughput and accuracy, in: Proceedings of the 33nd International Conference on Machine Learning, ICML 2016, Vol. 48 of JMLR Workshop and Conference Proceedings, 2016, pp. 201–210.
  * [16] F. Bourse, M. Minelli, M. Minihold, P. Paillier, Fast homomorphic evaluation of deep discretized neural networks, in: Advances in Cryptology - CRYPTO 2018, Vol. 10993 of Lecture Notes in Computer Science, 2018, pp. 483–512.
  * [17] F. Boemer, Y. Lao, R. Cammarota, C. Wierzynski, ngraph-he: a graph compiler for deep learning on homomorphically encrypted data, in: Proceedings of the 16th ACM International Conference on Computing Frontiers, CF 2019, 2019, pp. 3–13.
  * [18] Y. Tian, L. Njilla, J. Yuan, S. Yu, Low-latency privacy-preserving outsourcing of deep neural network inference, IEEE Internet Things J. 8 (5) (2021) 3300–3309.
  * [19] R. Dathathri, O. Saarikivi, H. Chen, K. Laine, K. E. Lauter, S. Maleki, M. Musuvathi, T. Mytkowicz, CHET: an optimizing compiler for fully-homomorphic neural-network inferencing, in: Proceedings of the 40th ACM SIGPLAN Conference on Programming Language Design and Implementation, PLDI 2019, 2019, pp. 142–156.
  * [20] H. Chen, R. Gilad-Bachrach, K. Han, Z. Huang, A. Jalali, K. Laine, K. E. Lauter, Logistic regression over encrypted data from fully homomorphic encryption, IACR Cryptol. ePrint Arch. 2018 (2018) 462.
  * [21] M. Kim, Y. Song, S. Wang, Y. Xia, X. Jiang, Secure logistic regression based on homomorphic encryption, IACR Cryptol. ePrint Arch. 2018 (2018) 74.
  * [22] J. H. Cheon, D. Kim, Y. Kim, Y. Song, Ensemble method for privacy-preserving logistic regression based on homomorphic encryption, IEEE Access 6 (2018) 46938–46948.
  * [23] H. Chen, W. Dai, M. Kim, Y. Song, Efficient multi-key homomorphic encryption with packed ciphertexts with application to oblivious neural network inference, in: Proceedings of the 2019 ACM SIGSAC Conference on Computer and Communications Security, CCS 2019, 2019, pp. 395–412.
  * [24] M. Clear, C. McGoldrick, Multi-identity and multi-key leveled FHE from learning with errors, in: Advances in Cryptology - CRYPTO 2015, Vol. 9216 of Lecture Notes in Computer Science, 2015, pp. 630–656.
  * [25] P. Mukherjee, D. Wichs, Two round multiparty computation via multi-key FHE, in: Advances in Cryptology - EUROCRYPT 2016, Vol. 9666 of Lecture Notes in Computer Science, 2016, pp. 735–763.
  * [26] C. Peikert, S. Shiehian, Multi-key FHE from lwe, revisited, in: Theory of Cryptography - 14th International Conference, TCC 2016-B, Vol. 9986 of Lecture Notes in Computer Science, 2016, pp. 217–238.
  * [27] Z. Brakerski, R. Perlman, Lattice-based fully dynamic multi-key FHE with short ciphertexts, in: Advances in Cryptology - CRYPTO 2016, Vol. 9814 of Lecture Notes in Computer Science, 2016, pp. 190–213.
  * [28] H. Lee, J. Park, On the security of multikey homomorphic encryption, in: Cryptography and Coding - 17th IMA International Conference, IMACC 2019, Vol. 11929 of Lecture Notes in Computer Science, 2019, pp. 236–251.
  * [29] I. Chillotti, N. Gama, M. Georgieva, M. Izabachène, Faster fully homomorphic encryption: Bootstrapping in less than 0.1 seconds, in: Advances in Cryptology - ASIACRYPT 2016, Vol. 10031 of Lecture Notes in Computer Science, 2016, pp. 3–33.
  * [30] I. Chillotti, N. Gama, M. Georgieva, M. Izabachène, Faster packed homomorphic operations and efficient circuit bootstrapping for TFHE, in: Advances in Cryptology - ASIACRYPT 2017, Vol. 10624 of Lecture Notes in Computer Science, Springer, 2017, pp. 377–408.
  * [31] I. Chillotti, N. Gama, M. Georgieva, M. Izabachène, TFHE: fast fully homomorphic encryption over the torus, J. Cryptol. 33 (1) (2020) 34–91.
  * [32] C. Gentry, A. Sahai, B. Waters, Homomorphic encryption from learning with errors: Conceptually-simpler, asymptotically-faster, attribute-based, in: Advances in Cryptology - CRYPTO 2013, Vol. 8042 of Lecture Notes in Computer Science, 2013, pp. 75–92.
  * [33] H. T. Vergos, An efficient bist scheme for non-restoring array dividers, in: 10th Euromicro Conference on Digital System Design Architectures, Methods and Tools (DSD 2007), IEEE, 2007, pp. 664–667.
  * [34] Z. Brakerski, C. Gentry, V. Vaikuntanathan, (leveled) fully homomorphic encryption without bootstrapping, ACM Trans. Comput. Theory 6 (3) (2014) 13:1–13:36.
  * [35] J. H. Cheon, A. Kim, M. Kim, Y. S. Song, Homomorphic encryption for arithmetic of approximate numbers, in: Advances in Cryptology - ASIACRYPT 2017, Vol. 10624 of Lecture Notes in Computer Science, 2017, pp. 409–437.
  * [36] J. H. Cheon, K. Han, A. Kim, M. Kim, Y. Song, Bootstrapping for approximate homomorphic encryption, in: Advances in Cryptology - EUROCRYPT 2018, Vol. 10820 of Lecture Notes in Computer Science, 2018, pp. 360–384.
  * [37] C. Boura, N. Gama, M. Georgieva, D. Jetchev, CHIMERA: combining ring-lwe-based fully homomorphic encryption schemes, J. Math. Cryptol. 14 (1) (2020) 316–338.
  * [38] C. Boura, N. Gama, M. Georgieva, Chimera: a unified framework for b/fv, TFHE and HEAAN fully homomorphic encryption and predictions for deep learning, IACR Cryptol. ePrint Arch. 2018 (2018) 758.


# Deterministic Preparation of Non-Gaussian Quantum States: Applications in
Quantum Information Protocols

Anindya Banerji1 1-Centre for Quantum Technologies, National University of
Singapore, Singapore.  Graciana Puentes2,3 2-Departamento de Fsica, Facultad
de Ciencias Exactas y Naturales, Universidad de Buenos Aires, Ciudad
Universitaria, 1428 Buenos Aires, Argentina
3-CONICET-Universidad de Buenos Aires, Instituto de Fsica de Buenos Aires
(IFIBA), Ciudad Universitaria, Buenos Aires, Argentina.

###### Abstract

We report a scheme for deterministic preparation of non-Gaussian quantum
states on-demand. In contrast to probabilistic approaches for preparation of
non-Gaussian quantum states, conditioned on photon subtraction or addition, we
present a scheme that can prepare non-Gaussian quantum states on-demand, by
applying a unitary transformation which removes the Gaussianity of measurement
statistics of field quadratures, namely a quadrature rotation via transmission
through a beam-splitter, using a two-mode photon-number squeezed state as
input. The resulting state exhibits a quantum vortex structure in quadrature
space, confirming its non-Gaussian nature. Such non-Gaussian quantum state
also reveals increased entanglement content, as quantified by the Logarithmic
Negativity and the Wigner function negative volume, therefore displaying high
potential for applications in quantum information protocols, in particular for
applications in entanglement distillation schemes.

††preprint: APS/123-QED

##  I Introduction

Gaussian quantum states are defined as the states for which measurement
statistics of field quadratures are Gaussian. As such, they can be fully
described by their mean field and covariance matrix. In the complete space of
states of Continuous Variable (CV) systems, Gaussian quantum states play a key
role, namely of all quantum states with a given covariance matrix, Gaussian
states have the least entanglement and the highest entropy [1, 2, 3]. From a
theoretical stand, Gaussian quantum states provide for a standard framework
for quantum information theory. On an experimental level, CV quantum
information has long been promoted due to the capability of on-demand
generation of large entangled states using either time-frequency modes [4, 6,
5, 9, 7, 8] or spatial modes [10, 12, 11].

Despite the experimental and theoretical advantages of Gaussian states, they
present major limitations in the context of quantum technologies: all Gaussian
measurements of such states can be efficiently simulated [15, 14, 13]. Pioneer
works on CV quantum computing argue that a non-Gaussian operation, meaning an
operation removing the Gaussian statistics of the states, is required for
implementation of a universal CV quantum computer [16]. Later works laying the
groundwork for CV quantum computation have left the question of such non-
Gaussian operation somewhat open [19, 17, 18]. Common schemes based on cubic
phase gates, turn out to be particularly hard to implement in an experimental
setup [21]. Furthermore, these protocols require highly non-Gaussian states to
encode information [20]. In spite the fact that such quantum states could
serve as a non-Gaussian resource for implementation of non-Gaussian gates
[23], they remain notoriously challenging to prepare. Despite the practical
problems involved with the non-Gaussian regime, non-Gaussian states are
expected to provide for enhanced entanglement content, and eventually it is
expected to be necessary to venture into non-Gaussian territory to reach a
quantum computational advantage, in the CV regime [22].

Experimentally, the preparation of a set of modes in a non-Gaussian quantum
state is generally much harder than the preparation of their Gaussian
counterparts. In essence, it suffices to apply a non-Gaussian unitary
operation to create a quantum state with non-Gaussian statistics of field
quadratures. In practice, such non-Gaussian unitary transformations are hard
to come by, meaning that often different preparation techniques are required.
There are two main approaches to reach non-Gaussianity. The first approaches
are probabilistic, and rely in performing non-Gaussian measurements on a
Gaussian state conditioned to a certain measurement outcome, such is the case
of photon-subtraction or photon-addition techniques. The second approach
concentrates on deterministic methods, which rely on the implementation of
non-Gaussian unitary transformations on-demand, meaning a unitary
transformation that removes the Gaussianity of the quantum states in a
deterministic fashion.

In this article, we introduce a deterministic method for preparation non-
Gaussian quantum states on demand, by implementation of a unitary
transformation, namely transmission through a beam-splitter, on an initial
two-mode photon-number squeezed state. The transformed state results in a
superposition of two-mode photon-subtracted Fock states, and therefore it
displays non-Gaussian statistics, more specific a non-Gaussian probability
distribution in quadrature space. In particular, the transformed state consist
of a quantum vortex structure, meaning a quantum state whose wave-function
takes the general form [24]:

$\psi(x,y)=(x-iy)^{m}e^{(x^{2}-y^{2})/2\sigma^{2}},$ (1)

where $m$ is an integer. It is well known that quantum vortex states are a
particular class of non-Gaussian quantum states [24, 25], and are therefore
amenable to applications requiring to venture into the non-Gaussian domain.
Moreover, we show that the non-Gaussian quantum state we propose presents
enhanced entanglement content, as quantified by the enhanced Logarithmic
Negativity of their Wigner function probability distribution. We argue that
such non-Gaussian quantum states can find relevant applications in the context
of quantum information protocols requiring non-Gaussian resources, in
particular in entanglement distillation schemes.

The article is structured as follows: In Section II, we introduce our
technique for deterministic preparation of non-Gaussian quantum vortex states
on-demand. Second, in Section III, we present an exhaustive investigation of
the Wigner function representation of the resulting non-Gaussian quantum
vortex state. In Section IV, we present a complete study of the entanglement
content in the proposed non-Gaussian quantum states, demonstrating that the
non-Gaussian nature of the quantum vortex state enhances the quantum
correlations existing in the initial two-mode squeezed Gaussian state, and are
therefore suitable candidates for quantum information protocols. Next, in
Section V, we propose a concrete application of such non-Gaussian quantum
states in quantum distillation protocols. Finally, in Section VI we outline
our conclusions.

##  II Non-Gaussian Quantum Vortex states

Figure 1: Phase profile of resulting non-Gaussian quantum vortex state in
quadrature representation for a squeezing parameter r = 0.02, exploring the
impact of the photon-number ($N/2$ per mode) in the formation of vortices.
Insets correspond to amplitude plots The numerical results confirm creation of
$N/2$ vortices for N total input photons. (a) $N/2=3$, (b) $N/2=4$, (c)
$N/2=5$, (d) $N/2=6$, (e) $N/2=7$, (f) $N/2=8$, with $N/2$ input photons per
mode (figure reproduced under the terms and conditions of the Creative Commons
Attribution License [26]).

The non-Gaussian quantum vortex states that we deal with in this article are
generated from truncated two-mode photon-number squeezed states. In the Fock
representation, these states are written as

$|\psi\rangle=\frac{A}{\cosh{r}}\sum_{j=0}^{N}\left(\tanh{r}\right)^{j}|j\rangle_{a}|j\rangle_{b}$
(2)

where $(a,b)$ are mode labels, $r$ is the squeezing parameter and $A$ is an
additional normalization factor. Similar states are routinely generated in
spontaneous parametric down conversion (SPDC) processes where generally an
infinite series is considered with diminishing probability of the higher order
photon-number terms. In fact, the pump power in such processes is kept
sufficiently small to ensure significantly less chances of occurrence of the
higher order terms. In such scenario, Eq. 2 can be treated as the more
practical representation of SPDC output. The additional normalization factor
is significant in this case since it arises due to such truncation of the
infinite series. In what follows, we treat $(a,b)$ as distinct spatial modes.

Next, we want to study the effect of beam splitter transformation on Eq. 2.
The two spatial modes are directed to the two different inputs of the beam
splitter. We map the input modes $(a,b)$ to the output modes
$(a^{\prime},b^{\prime})$ of the beam splitter as follows

$\displaystyle a$ $\displaystyle\rightarrow$ $\displaystyle
a^{\prime}=\frac{1}{\sqrt{2}}\left(a-ib\right)$ (3) $\displaystyle b$
$\displaystyle\rightarrow$ $\displaystyle
b^{\prime}=\frac{1}{\sqrt{2}}\left(b-ia\right)$ (4)

This transforms the initial photon-number squeezed state to the following

$\displaystyle|\psi^{\prime}\rangle$ $\displaystyle=$
$\displaystyle\frac{A}{\cosh{r}}\sum_{j=0}^{N}\tanh{r}^{j}\frac{j!}{2^{N/2}}$
(5) $\displaystyle\times$
$\displaystyle\sum_{k,l}^{j}i^{l+k}C_{k,l}|j-(l-k)\rangle|j+(l-k)\rangle$

where $C_{k,l}=\frac{\sqrt{(j-l+k)!(j+l-k)!}}{k!(j-k)!l!(j-l)!}$. We call
$|\psi^{\prime}\rangle$ the non-Gaussian quantum vortex state. The transformed
state can be regarded as a superposition of photon-subtracted Fock states, due
to the imbalance in photon-number between the two modes, therefore its non-
Gaussian character is apparent. Since a beam splitter transformation is a
unitary process, there is no change in the total number of photons in the two
modes combined between the photon-number squeezed state $|\psi\rangle$ and the
non-Gaussian quantum vortex state. Rather, it results in a redistribution of
the photons between the two modes for each term of the summation in Eq. 2. The
nomenclature would become clear if we look at the associated quadrature
distribution. Using the relation

$\psi_{n}\left(x\right)=\langle
x|n\rangle=\frac{1}{\sqrt{2^{n}n!\sqrt{\pi}}}\exp\left(-\frac{x^{2}}{2}\right)H_{n}(x)$
(6)

where $H_{n}(x)$ is the Hermite polynomial of order $n$ and $|n\rangle$ is a
number state, and noting that product of Hermite polynomials
$H_{n}(x)H_{m}(y)$ can be written in terms of Laguerre polynomials
$L^{n-m}_{m}(x^{2}+y^{2})$, the quadrature distribution can be derived as
follows [26]

$\displaystyle\psi^{\prime}(x,y)$ $\displaystyle=$
$\displaystyle\frac{A^{\prime}}{\cosh{r}}\sum_{j=0}^{N}\sum_{m=0}^{j}\tanh{r}^{j}\frac{j!}{2^{N/2}}C_{k,l}$
(7) $\displaystyle\times$ $\displaystyle
L_{j-m}^{2m}\left(x^{2}+y^{2}\right)\exp\left(-\frac{x^{2}+y^{2}}{2}\right)$

where $m=l-k$. This illustrates the fact the state after the beam splitter
transformation consists of a finite of quantum vortices of order $2m$. This
justifies the nomenclature. As a further illustration, we showcase a few cases
of quantum vortex states in Fig. 1. It is quite evident from the amplitude
plots that these states exhibit a departure from their initial Gaussian-like
characteristics.

##  III Wigner Function Representation of Non-Gaussian Quantum Vortex States

In the previous section, we have outlined the method to generate the non-
Gaussian quantum vortex states and looked at their associated quadrature
distribution. In this section, we want to study their phase-space distribution
to better understand the non-Gaussian structure inherent in these states. In
order to do that, we will use the Wigner distribution function. The Wigner
distribution is a quasi-probability distribution function that is real, non-
singular and produces accurate quantum mechanical operator averages. It is
also characterised by a definite marginal distribution. A remarkable advantage
of using the Wigner function is that it can reveal both non-classicality and
non-Gaussianity of quantum states.

Figure 2: Two dimensional slice of the Wigner function of the non-Gaussian
quantum vortex state for $N=6$ and different values of $r$. The non-Gaussian
nature is evident from the negative regions in the above plots. It is to be
noted that the negativity arises and then increases with increasing $r$. (a),
(b) and (c) are projections of the Wigner function at the plane
$\left(y=0,p_{x}=0\right)$ while (d), (e) and (f) are projections of the
Wigner function at the plane $\left(x=0,p_{y}=0\right)$. Inset are the
corresponding contour plots for each of the cases. Details are in the text.

For a single mode quantum state $|n\rangle$ in the Fock space representation,
the associated Wigner function is defined as

$W\left(x,p_{x}\right)=\frac{2}{\pi}\left(-1\right)^{n}L_{n}(4q^{2})e^{-2q^{2}}$
(8)

where $q^{2}=x^{2}+p_{x}^{2}$ and $L_{n}(4q^{2})$ is again Laguerre polynomial
of order $n$. $\left(x,p_{x}\right)$ are the quadrature variables. The Wigner
function of the two-mode state $|n\rangle_{a}|m\rangle_{b}$ is then simply the
product of the corresponding Wigner functions as follows

$\displaystyle W\left(x,p_{x},y,p_{y}\right)$ $\displaystyle=$
$\displaystyle\frac{4}{\pi^{2}}\left(-1\right)^{n+m}$ (9)
$\displaystyle\times$ $\displaystyle
L_{n}(4q_{a}^{2})e^{-2q_{a}^{2}}L_{m}(4q_{b}^{2})e^{-2q_{b}^{2}}$

where $q_{a}^{2}=x^{2}+p_{x}^{2}$, $q_{b}^{2}=y^{2}+p_{y}^{2}$ while
$\left(x,p_{x}\right)$ and $\left(y,p_{y}\right)$ are the quadrature variables
for modes $a$ and $b$ respectively.

Figure 3: Two dimensional slice of the Wigner function of the non-Gaussian
quantum vortex state for a fixed value of $r$ and different value of $N$. The
non-Gaussian nature is again evident from the negative regions in the above
plots. (a), (b) and (c) are projections of the Wigner function at the plane
$\left(y=0,p_{x}=0\right)$ while (d), (e) and (f) are the projections of the
Wigner function at the plane $\left(x=0,p_{y}=0\right)$. Inset are the
corresponding contour plots. Details are in the text.

In our case, the non-Gaussian quantum vortex state is obtained by performing
an unitary operation on an initial state in photon-number superposition. Using
its Fock space representation of Eq. (5), after some manipulation, we can
write the Wigner function in the following form

$W=\frac{A^{\prime\prime}}{\cosh{r}^{2}}\sum_{j=0}^{N}\sum_{m=0}^{j}\tanh{r}^{2*j}\frac{j!^{2}}{2^{N}}|C_{k,l}|^{2}W_{j,m}$
(10)

where $A^{\prime\prime}$ is the normalization term and $W_{j,m}$ are defined
as

$\displaystyle W_{j,m}$ $\displaystyle=$
$\displaystyle\frac{(-1)^{2*j}}{\pi}L_{j+m}[4(Q_{0}+Q_{1})]$ (11)
$\displaystyle\times$ $\displaystyle
L_{j-m}[4(Q_{0}-Q_{1})]\exp{\left(-4Q_{0}\right)}$

where $Q_{0}=\frac{1}{4}\left(x^{2}+y^{2}+p_{x}^{2}+p_{y}^{2}\right)$ and
$Q_{1}=\frac{xp_{y}-yp_{x}}{2}$. We study projections of the Wigner function
of the non-Gaussian quantum vortex state on different planes in Fig. 2 and
Fig. 3. Since the Wigner function is 4-dimensional for a two-mode state, it is
impossible to graphically reproduce the entire structure.We therefore chose to
highlight only those planes which best illustrates the negative volume of the
Wigner function. As can be seen from the figures, the projections at the
planes $\left(y=0,p_{x}=0\right)$ and $\left(x=0,p_{y}=0\right)$ exhibit
maximum negativity. The negative volume increases with increasing $r$ and $N$
though, the squeezing parameter has greater influence. Even for high values of
$N$, there is no negative region for low values of $r$ as is shown in Fig, 2a
and 2d. Comparing this with Fig. 3a and 3d, we see that even for the case with
minimum number of photons, negative regions are present for sufficiently high
values of $r$. In order to better quantify these results, we compute the total
negative volume of the entire Wigner function below.

Now, as has been mentioned before, an important aspect of the Wigner function
is its usefulness in detecting nonclassicality and non-Gaussianity of quantum
states. In general, the presence of negative regions in the Wigner function is
accepted as a signature of nonclassicality. Additionally, the Wigner function
for a pure quantum state that is Gaussian is strictly positive [30]. This
means that the negativity of the Wigner function is a witness of non-
Gaussianity. This can be quantified with the help of the negativity volume of
the Wigner function defined as

$\mathrm{NV}=\frac{1}{2}\left(\int_{V}|W|\text{d}V-1\right)$ (12)

where the integration is performed over the entire phase space. We have
performed a numerical integration to calculate the negative volume of the non-
Gaussian quantum vortex states. We study the results in Fig. 4.

Figure 4: Negativity volume of the Wigner function of the non-Gaussian quantum
vortex states for different values of $N$. The black markers are the
datapoints from numerical evaluation of Eq. 12. The continuous curves are
numerical best fits.

We see that the negative volume increases with increasing values of $r$ and
saturates to a maximum value. This has two related interpretations. It means
that increased squeezing leads to increases non-classicality which is already
known. In addition, it also signifies that an increase in the squeezing
parameter also leads to increased non-Gaussianity of the final state.

Figure 5: _Top row_ : Logarithmic negativity of the non-Gaussian quantum
vortex states and two-mode photon-number squeezed states for different values
of N in arbitrary units. _Bottom row_ : Ratio of the logarithmic negativities
for the two states better illustrates the increase in entanglement due to the
beam splitter transformation. Here
$R=\frac{\mathcal{L}\left(\rho\right)}{\mathcal{L}\left(\rho^{\prime}\right)}$
where $\rho^{\prime}$ is the density matrix of the two-mode photon-number
squeezed state before the beam splitter transformation. (a) and (d) $N=2$, (b)
and (e) $N=4$, (c) and (f) $N=6$.

We should also point out that the negativity increases with increasing number
of photons as is evident from the plots (red for $N=2$ and blue for $N=4$).
This means that the redistribution of photon number between the two modes
effected by the beam splitter transformation is successful in increasing the
nonclassical effects of the photon-number squeezed states. More interestingly,
it shows that the transformation also leads to the generation of significant
non-Gaussianity.

The above results serve as a pointer to investigate how the entanglement
between the two modes are influenced by this redistribution of photon numbers
which we study in the next section.

##  IV Entanglement Content of Non-Gaussian Quantum Vortex States

In this section we study the entanglement content of non-Gaussian quantum
vortex states. We use Logarithmic Negativity [31] which is easily computable
and a proven entanglement monotone. It is defined as

$\mathcal{L}\left(\rho\right)=\log_{2}\left(1+2\mathrm{N}\right)$ (13)

where $\rho$ is the density matrix of the quantum state and $\mathrm{N}$ is
the negativity measure [32], defined as

$\mathrm{N}=\frac{||\rho^{T}||-1}{2}$ (14)

Here $||.||$ denotes trace norm and $\rho^{T}$ is the partial transpose of the
density matrix of the quantum state with respect to one of the subsystems. The
negativity then is the absolute values of the sum of all negative eigenvalues
of the partially transposed density matrix. In order to calculate
$\mathcal{L\left(\rho\right)}$ of the non-Gaussian quantum vortex state, we
first need to write the density matrix corresponding to
$|\psi^{\prime}\rangle$. It has the form

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle|\psi^{\prime}\rangle\langle\psi^{\prime}|$ $\displaystyle=$
$\displaystyle\frac{|A|^{2}}{\cosh^{2}r}\sum_{j,j^{\prime}}^{N}\tanh{r}^{j+j^{\prime}}\frac{j!j^{\prime}!}{2^{N}}\sum_{k,l,k^{\prime},l^{\prime}}^{j,j^{\prime}}i^{k+l}(-i)^{k^{\prime}+l^{\prime}}C_{k,l}C_{k^{\prime},l^{\prime}}$
$\displaystyle\times$ $\displaystyle|j-(l-k)\rangle\langle
j^{\prime}-(l^{\prime}-k^{\prime})|\otimes|j+(l-k)\rangle\langle
j^{\prime}+(l^{\prime}-k^{\prime})|$ (15)

The partial transpose of Eq. IV is obtained by standard techniques and then
diagonalized numerically to calculate the negative eigenvalues which are then
used to determine $\mathcal{L}(\rho)$. A similar process is followed to
calculate the logarithmic negativity of the two-mode photon-number squeezed
state of Eq. 2 for which the density matrix has the form

$\rho^{\prime}=\frac{|A|^{2}}{\cosh^{2}r}\sum_{j,j^{\prime}}^{N}\left(\tanh{r}\right)^{j+j^{\prime}}|j\rangle\langle
j^{\prime}|\otimes|j\rangle\langle j^{\prime}|$ (16)

We compare the values of $\mathcal{L}\left(\rho\right)$ and
$\mathcal{L}\left(\rho^{\prime}\right)$ as well as study their ratio $R$ as a
function of squeezing parameter $r$ in Fig. 5.

We see that the entanglement increases after the beam-splitter operation. This
is due to the fact that the uncertainty in photon number in each mode
increases due to the redistribution of photons between the two modes mediated
by the beam-splitter. The total entanglement increases with increasing number
of photons as well. Also, as $r$ increases, the entanglement increases sharply
before saturating to a maximum value. More importantly, entanglement vanishes
in the absence of squeezing since there are no photons present in either modes
in such situation.

##  V Applications of Non-Gaussian Quantum Vortex States in Quantum
Information Protocols

###  V.1 Entanglement Distillation Protocols

For finite-dimensional systems, the term entanglement distillation has been
linked to the notion that one can obtain highly entangled states by means of
local quantum operations and classical communication, by starting from a large
number of weakly entangled quantum states and ending with a smaller number of
more entangled ones. Such approaches also work as the basis for quantum
cryptographic schemes. An equivalent procedure should also exists for the
distillation of Gaussian states by means of local Gaussian operations and
classical communication [28]. At any rate, entanglement distillation aims at
producing more highly entangled states out of a situation where entanglement
is present only in a noisy form, presumably as a consequence of some lossy
quantum channel. Entanglement distillation can be regarded as a key element in
quantum repeater approaches, allowing for long-range entanglement distribution
in the presence of noise. In essence, it is possible to differentiate between
distillation protocols that involve several copies of an entangled state at
each step of the scheme, and local filtering approaches that take a single
specimen of a state and, under appropriate filtering, give rise to a more
entangled state. In the context of Gaussian operations, CV entanglement
distillation of neither kind is possible without the aid of non-Gaussian
operations, such as photon addition or subtraction [28].

###  V.2 Proposed Experimental Scheme

Figure 6: Proposed experimental scheme for implementation of entanglement
distillation protocol. (a) two-mode Gaussian entangled state produced in a
non-linear crystal (NLC) via non-collinear SPDC. (b) Deterministic non-
Gaussian quantum state prepared via transmission through a beam-splitter BS.
($a,b$) and ($a^{\prime},b^{\prime}$) label the input and output spatial modes
transmitted by the BS, respectively. Such unitary operation enhances the
entanglement content. (c) Entanglement distillation quantified by a partial
detection approach via weak homodyne detection using photon-number-resolving
(PNR) detectors and weak local oscillators (LO). Selected joint POVM elements
$\Pi_{ab}$ can provide for accurate bounds on the entanglement content via
Convex Optimization approaches [27] (see text for details).

In order to distill CV entanglement from Gaussian states, such as the two-mode
photon-number squeezed state described by in Section II, an operation that
removes the Gaussianity of the field quadrature statistics is required.
Examples of such non-Gaussian operations include the conditional subtraction
or addition of a photon [1-11]. However, one main limitation of such
conditional operations is their probabilistic nature. Here we propose a
deterministic approach for removing the Gaussian statistics of the initial
two-mode photon-number squeezed state, such deterministic unitary operation
consists of transmission through a standard beam-splitter (BS). It has
recently been demonstrated that such operation introduces photon-number
fluctuations [26], which result in the creation of quantum vortices in the
quadratures with non-Gaussian statistics, therefore creating a non-Gaussian
quantum state suitable for applications in quantum information protocols, such
as entanglement distillation.

The proposed experimental scheme for entanglement distillation is depicted in
Fig. 6. It consists of three main steps. The first step (Fig. 6 (a)) is the
preparation of the initial Gaussian entangled state, via non-collinear SPDC,
using a non-linear crystal (NLC). As anticipated in Section II, for a
sufficiently attenuated pump the resulting state can be approximated by a
truncated two-mode photon-number squeezed Gaussian (G) quantum state of the
form [27]:

$|\psi\rangle_{\mathrm{G}}=\frac{A}{\cosh{r}}\sum_{j=0}^{N}\left(\tanh{r}\right)^{j}|j\rangle_{a}|j\rangle_{b}$
(17)

where $A$ is a normalization constant and $r$ the squeezing parameter.

The second step (Fig. 6 (b)) consists of a non-Gaussian unitary operation,
meaning an operation that removes the Gaussian measurement statistic of field
quadratures in the initial state, resulting in a non-Gaussian (NG) quantum
state. As described in detail in Ref. [26], the beam-splitter introduces
photon-number fluctuations resulting in a quantum state with a binomial
photon-number distribution of the form:

$\displaystyle|\psi\rangle_{\mathrm{NG}}$ $\displaystyle=$
$\displaystyle\frac{D}{\cosh{r}}\sum_{j=0}^{N/2}A^{r,N}_{j}\times$
$\displaystyle\sum_{k=0}^{j}\sum_{l=0}^{j}B^{\phi}_{k,l}C^{Nj}_{kl}|j-(l-k),j+(l-k)\rangle.$

where $D$ is the normalization factor. Explicit expressions for the
coefficients $A^{r,N}_{j},B^{\phi}_{k,l},C^{Nj}_{kl}$ are given in Section II.
Such states exhibits a vortex structure in quadrature space (Fig. 1), and non-
Gaussian field quadrature statistics.

As reported in Section III and IV, the resulting non-Gaussian state displays
enhanced entanglement content, which is later on quantified by a partial
detection approach, based on Entanglement Witnesses and Convex Optimization
schemes [27]. The non-Gaussian quantum state is eventually routed towards a
weak homodyne detection station, using standard silver mirrors ($M$).

The final step (Fig.6 (c)) in the entanglement distillation protocol consists
of a partial detection approach, implemented by constructing suitable
Entanglement Witnesses, which in turn are built by selecting suitable joint
Positive Operator Valued Measures (POVM) $\Pi_{ab}=\Pi_{a}\otimes\Pi_{b}$ for
each detector $D_{a,b}$, where $a,b$ label each spatial mode transmitted by
the BS. Photons are detected using photon number resolving (PNR) detectors and
weak local oscillators. The POVM elements of such weak homodyne detectors have
been fully characterized in Ref [29]. The measurement outcomes provided by the
selected POVM elements can provide for accurate bounds on the entanglement
content via Convex Optimization approaches. A full description of such Convex
Optimization schemes is reported in Ref. [27].

##  VI Discussion

We presented a scheme for deterministic preparation of non-Gaussian quantum
states on-demand. In contrast to the standard probabilistic approaches for
preparation of non-Gaussian states, conditioned on photon subtraction or
addition, the scheme presented here can prepare non-Gaussian quantum states
on-demand by applying a deterministic unitary transformation which removes the
Gaussianity of quadrature statistics of the initial state, namely a quadrature
rotation via transmission through a standard beam-splitter using a photon-
number squeezed state as input. The resulting non-Gaussian quantum state
consists of a superposition of photon-subtracted Fock states, and it exhibits
a quantum vortex structure in quadrature space, thus confirming the non-
Gaussian character of measurement statistics of field quadratures, it also
reveals increased entanglement content, as quantified by the Logarithmic
Negativity and the Wigner function negative volume, therefore displaying high
potential for applications in quantum information protocols such as
entanglement distillation schemes.

##  VII Acknowledgements

The authors gratefully acknowledge Jens Eisert for helpful discussions. G.P.
acknowledges financial support via PICT Startup.

## References

  * [1] M. M. Wolf, G. Giedke, and J. I. Cirac, Phys. Rev. Lett. 96, 080502 (2006).
  * [2] A. S. Holevo, M. Sohma, and O. Hirota, Phys. Rev. A 59, 1820 (1999).
  * [3] M. Walschaer, arXiv:2104.12596 [quant-ph] (2021).
  * [4] J. Roslund, R. M. de Araujo, S. Jiang, C. Fabre, and N. Treps, ´ Nat Photon 8, 109 (2014).
  * [5] Y. Cai, J. Roslund, G. Ferrini, F. Arzani, X. Xu, C. Fabre, and N. Treps, Nat. Commun. 8, 15645 (2017).
  * [6] M. Chen, N. C. Menicucci, and O. Pfister, Phys. Rev. Lett. 112, 120505 (2014).
  * [7] M. V. Larsen, X. Guo, C. R. Breum, J. S. Neergaard-Nielsen, and U. L. Andersen, Science 366, 369 (2019).
  * [8] Z. Yang, M. Jahanbozorgi, D. Jeong, S. Sun, O. Pfister, H. Lee, and X. Yi, arXiv:2103.03380 [physics.optics] 2021.
  * [9] W. Asavanant, Y. Shiozawa, S. Yokoyama, B. Charoensombutamon, H. Emura, R. N. Alexander, S. Takeda, J.-i. Yoshikawa, N. C. Menicucci, H. Yonezawa, and A. Furusawa, Science 366, 373 (2019).
  * [10] M. Yukawa, R. Ukai, P. van Loock, and A. Furusawa, Phys. Rev. A 78, 012301 (2008).
  * [11] D. Barral, M. Walschaers, K. Bencheikh, V. Parigi, J. A. Levenson, N. Treps, and N. Belabas, Phys. Rev. Applied 14, 044025 (2020).
  * [12] X. Su, Y. Zhao, S. Hao, X. Jia, C. Xie, and K. Peng, Opt. Lett. 37, 5178 (2012).
  * [13] S. D. Bartlett, B. C. Sanders, S. L. Braunstein, and K. Nemoto, Phys. Rev. Lett. 88, 097904 (2002).
  * [14] J. E. Bourassa, R. N. Alexander, M. Vasmer, A. Patil, I. Tzitrin, T. Matsuura, D. Su, B. Q. Baragiola, S. Guha, G. Dauphinais, K. K. Sabapathy, N. C. Menicucci, and I. Dhand, Quantum 5, 392 (2021).
  * [15] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.- C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, P. Hu, X.-Y. Yang, W.-J. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N.-L. Liu, C.-Y. Lu, and J.-W. Pan,Science 370, 1460 (2020).
  * [16] S. Lloyd and S. L. Braunstein, Phys. Rev. Lett. 82, 1784 (1999).
  * [17] M. Gu, C. Weedbrook, N. C. Menicucci, T. C. Ralph, and P. van Loock, Phys. Rev. A 79, 062318 (2009).
  * [18] N. C. Menicucci, Phys. Rev. Lett. 112, 120504 (2014).
  * [19] N. C. Menicucci, P. van Loock, M. Gu, C. Weedbrook, T. C. Ralph, and M. A. Nielsen, Phys. Rev. Lett. 97, 110501 (2006).
  * [20] D. Gottesman, A. Kitaev, and J. Preskill, Phys. Rev. A 64, 012310 (2001).
  * [21] F. Arzani, N. Treps, and G. Ferrini, Phys. Rev. A 95, 052352 (2017).
  * [22] A. Mari and J. Eisert, Phys. Rev. Lett. 109, 230503 (2012)
  * [23] B. Q. Baragiola, G. Pantaleoni, R. N. Alexander, A. Karanjai, and N. C. Menicucci, Phys. Rev. Lett. 123, 200502 (2019).
  * [24] G. S. Agarwal, R. R. Puri, and R. P. Singh, Phys. Rev. A 56, 4207 (1997).
  * [25] G. S. Agarwal, New J. Phys. 13 073008 (2011).
  * [26] G. Puentes and A. Banerji, Frontiers in Physics 9, 690721 (2021).
  * [27] G. Puentes, A. Datta, A. Feito, J. Eisert, M. B. Plenio, and I. A Walmsley, New J. Phys. 12, 033042 (2010).
  * [28] D. E. Browne, J. Eisert, S. Scheel, and M. B. Plenio, Phys. Rev. A 67, 062320 (2003).
  * [29] G. Puentes, J. S Lundeen, M. PA Branderhorst, H. B Coldenstrodt-Ronge, B. J Smith, I. A Walmsley, Phys. Rev. Lett. 102 (8), 080404 (2009).
  * [30] R. L. Hudson, Rep. Math. Phys. 6, 249 (1974).
  * [31] M. B. Plenio, Phys. Rev. Lett. 95. 090503 (2005).
  * [32] G. Vidal and R. F. Werner, Phys. Rev. A 65 032314 (2002).