# Continuous-variable quantum repeaters based on bosonic error-correction and
teleportation: architecture and applications

Bo-Han Wu1, Zheshen Zhang2,4, Quntao Zhuang3,4 1 Department of Physics,
University of Arizona, Tucson, Arizona 85721, USA  2 Department of Materials
Science and Engineering, University of Arizona, Tucson, Arizona 85721, USA  3
Electrical and Computer Engineering, University of Arizona, Tucson, Arizona
85721, USA  4 J. C. Wyant College of Optical Sciences, University of Arizona,
Tucson, Arizona 85721, USA<EMAIL_ADDRESS>
(August 30, 2024)

###### Abstract

Quantum repeaters are essential ingredients for quantum networks that link
distant quantum modules such as quantum computers and sensors. Motivated by
distributed quantum computing and communication, quantum repeaters that relay
discrete-variable quantum information have been extensively studied; while
continuous-variable (CV) quantum information underpins a variety of quantum
sensing and communication application, a quantum-repeater architecture for
genuine CV quantum information remains largely unexplored. This paper reports
a CV quantum-repeater architecture based on CV quantum teleportation assisted
by the Gottesman-Kitaev-Preskill (GKP) code to significantly suppress the
physical noise. The designed CV quantum-repeater architecture is shown to
significantly improve the performance of CV quantum key distribution,
entanglement-assisted communication, and target detection based on quantum
illumination, as three representative use cases for quantum communication and
sensing.

  * Sep 2021

##  1 Introduction: file preparation and submission

Quantum networks [1-5] not only offer unconditional security in private-key
distributions [6-8], but also enable the establishment of entanglement across
multiple parties to endow quantum-enhanced capabilities. Photons are ideal
information carriers for long-haul quantum communications by virtue of their
robustness against environmental noise, but they are susceptible to loss
because, unlike classical information, quantum information cannot be
regenerated by amplifiers due to the quantum no-cloning theorem [9,10]. Such a
restriction places a fundamental rate-loss trade-off between entanglement-
distribution rate and transmission distance, which, in terms of the
distribution of bipartite entanglement, was formulated as the
Pirandola–Laurenza–Ottaviani–Banchi (PLOB) bound [11] and has been
subsequently generalized to end-to-end capacity of a general quantum network
[12].

To circumvent the rate-loss trade-off, a long-distance quantum link is divided
into shorter and less lossy links via introducing intermediate quantum
repeater (QR) nodes [13-17]. Based on the processing power at each node, QRs
are categorized into three generations (see Refs. [18-20], Refs. [15,21] and
Refs. [22,23]). The mainstream QR architectures have been dedicated to the
long-distance distribution of discrete-variable (DV) quantum states
[15,24-26], i.e., qubits, to link quantum computers, in analogy to sharing
digital information among classical computers. On the other hand, continuous-
variable (CV) quantum states, akin to analog information, underpins a variety
of quantum-enhanced sensing and communication capabilities including entangled
sensor networks [27-34], physical-layer quantum data classification [35,36],
quantum-illumination (QI) target detection [37-40] and ranging [41], and
entanglement-assisted (EA) communication [42-50]. Apart from a handful of
investigations for a few specific use cases [51,52], the QR architecture for
CV quantum states remains largely unexplored.

Figure 1: Scheme of $m$-relay repeaters based on CV error-correction protocol.
ENC: encoding. DEC: decoding. $L$ is the physical distance between Alice and
Bob. $L_{\Delta}$ is the inter-repeater spacing.

Quantum error correction (QEC) is an essential ingredient for QRs to reliably
relay quantum information. QEC for qubits has been well established to support
the development of fault-tolerant quantum computing [53,54]. QEC for QRs,
however, requires an additional framework to account for the infinite
dimensional Hilbert space that photons reside in. In this regard, bosonic QEC
[55] has emerged as a powerful paradigm to protect quantum information carried
on photons. To date, multiple single-mode bosonic codes, including the
binomial code [56], Schrödinger-cat-state codes [57], and _Gottesman-Kitaev-
Perskill_ (GKP) codes [58-63], have been proposed and experimentally produced
[64-72]. Most bosonic codes have been designed to protect qubits by encoding
them into bosonic modes. The more recent works of Rozpȩdek _et al._ [73] and
Fukui _et al._ [74,75] introduced GKP-formed qubit codes into the QR
architecture to transmit qubits, but a QR based on bosonic QEC to transmit CV
quantum information, which will significantly benefit a wide range of quantum-
enhanced applications, remains elusive.

This article proposes a CV QR architecture based on the recently developed
GKP-assisted QEC [76,77] combined with CV quantum teleportation [78,79] and is
organized as follows: Sec. 2 provides an overview on the QR architecture; Sec.
3 introduces the GKP-assisted QEC to the QR architecture. Finally, in Sec. 4,
the QEC QR is shown to boost the performance of EA communication, target
detection based on QI, and CV quantum key distribution (QKD).

##  2 Quantum-repeater architecture with bosonic quantum error correction

Fig. 1 illustrates the architecture for our CV QR based on the bosonic QEC
code [76]. Consider a quantum link comprising $m$ QR nodes. At the source,
Alice performs an _encoding_ operation on the message mode and an ancilla mode
and then transmits both modes to the first QR node through a quantum channel.
The QR node performs a _decoding_ operation on both received modes to correct
the accumulated errors incurred by the channel. Afterwards, _encoding_
operations are operated on the error-corrected message mode and an additional
ancilla mode; the two modes are subsequently transmitted to the next QR node
for decoding and encoding, until the message mode is finally decoded at Bob’s
terminal.

Note that here the quantum channels not only model the transmission via fiber
quantum links, but also takes into account some pre- and post-processing that
enhances the quantum information transmission. Each fiber link between two
nodes can be modeled as a bosonic pure-loss channel with the transmissivity
$\eta=10^{-\gamma L_{\Delta}/10}$, where $L_{\Delta}$ is the physical distance
between the two nodes, with an attenuation factor $\gamma=0.2$ decibels per
kilometer. With additional pre- and post-processing, we convert the pure-loss
link into two types of quantum channels, the amplified one-way channel (Sec.
2.1) and the quantum teleportation channel (Sec. 2.2). The effect of
transmitting the message and ancilla modes through the amplified one-way or
quantum teleportation channel is equivalent to adding to their quadratures
some additive noises of variance $\sigma^{2}_{A}$ or $\sigma^{2}_{T}$, instead
of the original pure-loss.

Figure 2: The scheme of (a) amplified one-way channel, and (b) teleportation
channel. CC: classical communication. HM: homodyne measurement. Amp:
amplification.

Figure 3: Variances of the additive noise for the protocols of amplified one-
way and teleportation channels. Amp: amplification. Tele: teleportation.

###  2.1 Amplified one-way channel

Sketched in Fig. 3(a), the amplified one-way channel introduced in the QR
architecture studied by Fukui _et al._ [74] applies a phase-insensitive
amplifier of gain $1/\eta$ before the pure-loss channel of transmissivity
$\eta$ induced by the fiber transmission. The variance of additive noise of
the amplified one-way channel is derived to be

$\sigma^{2}_{A}=1-\eta,$ (1)

i.e.,
$\braket{\hat{q}^{2}}_{\text{vac}}=\braket{\hat{p}^{2}}_{\text{vac}}=1/2$.
Because both the channel loss and the amplification add noise, the performance
of QEC is limited. To overcome the drawback of the amplified one-way channel,
we introduce the quantum teleportation channel below.

###  2.2 Quantum teleportation channel

CV quantum teleportation transmits CV quantum states from the source to the
destination with only local operations, classical communication (LOCC), and
shared CV entangled states. To implement a CV quantum teleportation channel in
the CV QR architecture, a two-mode squeezed vacuum (TMSV) source placed in the
middle of QR nodes, as shown in Fig. 3(b), generates entangled signal and
idler modes that are sent to two adjacent QR nodes through two pure-loss
channels, yielding a shared entangled state that is subsequently used to
teleport a CV quantum state between the two QR nodes. Earlier results of CV
quantum teleportation (e.g., Ref. [80]) showed that the teleportation channel
is equivalent to an additive thermal noise channel due to finite squeezing and
TMSV distribution loss. The variance of additive noise is

$\sigma^{2}_{T}=\sqrt{\eta}10^{-s/10}+\left(1-\sqrt{\eta}\right),$ (2)

where $s$ (i.e. unit dB) characterizes the squeezing level of TMSV.

Fig. 3 plots the additive noise of the amplified one-way channel (red) and the
teleportation channel (blue). Apparently, the inter-repeater spacing,
$L_{\Delta}$ is a crucial factor for determining the optimal transmission
protocol, and Fig. 3 implies there exists a minimal inter-repeater spacing
(MIRS)

$L^{*}_{\Delta}\equiv-2\left(\log_{10}\left[1-10^{-s/10}\right]\right)/\gamma,$
(3)

such that $\sigma^{2}_{T}<\sigma^{2}_{A}$, $\forall
L_{\Delta}>L^{*}_{\Delta}$.

##  3 GKP-error-correction code

Before proceeding to GKP-assisted QEC, we will first introduce the GKP ancilla
mode in Sec. 3.1 and the GKP-two-mode-squeezing code in Sec. 3.2.

###  3.1 The GKP state

A bosonic mode of, e.g., the photon or the phonon, encompasses the continuous
degrees of freedom in the position and momentum quadratures. Mathematically,
the quadratures, $\hat{q}$ and $\hat{p}$, are the normalized real and
imaginary parts of the annihilation operator $\hat{a}$,

$\hat{q}=\frac{1}{\sqrt{2}}\left(\hat{a}+\hat{a}^{\dagger}\right)\;,\;\;\hat{p}=\frac{1}{i\sqrt{2}}\left(\hat{a}-\hat{a}^{\dagger}\right),$
(4)

satisfying the commutation relation $\left[\hat{q},\hat{p}\right]=i$
($\hbar\equiv 1$ for simplicity). The GKP state is pure and stabilized by the
following CV analog of the Pauli-Z and Pauli-X operators:

$\hat{Z}=\hat{D}(0,\sqrt{2\pi})\;,\;\;\hat{X}=\hat{D}(\sqrt{2\pi},0),$ (5)

where
$\hat{D}\left(\alpha,\beta\right)=e^{i\left(\alpha\hat{p}-\beta\hat{q}\right)}.$
An ideal GKP state can be considered as the superposition of an infinite
number of position or momentum eigenstates along a grid, i.e.,

$\ket{\text{GKP}}\propto\sum_{n\in\mathbb{Z}}\ket{q=n\sqrt{2\pi}}\propto\sum_{n\in\mathbb{Z}}\ket{p=n\sqrt{2\pi}}.$
(6)

Figure 4: The Wigner functions in phase space of (a) ideal and (b) approximate
GKP states.

The Wigner function of the ideal GKP state is sketched in Fig. 4(a), where
each dot represents a Dirac delta function. A GKP state incorporates precise
information of both quadratures within a critical range without violating the
uncertainty principle. Precisely, the standard deviation of both quadrature
operators modulo $\sqrt{2\pi}$ are zero. Hence, both quadratures can be
measured simultaneously modulo $\sqrt{2\pi}$, rendering the GKP state perfect
to calibrate any other states encoded by the GKP codes. Nonetheless, ideal GKP
states are not normalizable and thus not physical. The consideration of
experimental feasibility calls for a CV QEC based on approximate GKP states,
as presented below.

The Wigner function of an approximate GKP state [71,72,81,82] is plotted in
Fig. 4(b). For an approximate GKP state, a series of Dirac delta functions in
Eq. (6) are replaced by a series of Gaussian packets weighted by a Gaussian
profile:

$\displaystyle\ket{\text{GKP}}\propto$
$\displaystyle\sum_{n\in\mathbb{Z}}e^{-\pi\sigma_{\text{G}}^{2}n^{2}}\int_{-\infty}^{\infty}e^{-\frac{\left(q-\sqrt{2\pi}n\right)^{2}}{2\sigma_{\text{G}}^{2}}}\ket{q}dq\;\;\propto$
$\displaystyle\sum_{n\in\mathbb{Z}}e^{-\pi\sigma_{\text{G}}^{2}n^{2}}\int_{-\infty}^{\infty}e^{-\frac{\left(p-\sqrt{2\pi}n\right)^{2}}{2\sigma_{\text{G}}^{2}}}\ket{p}dp.$
(7)

The linewidths of each Gaussian teeth is characterized by the squeezing
parameter $s^{(\text{G})}=-10\log_{10}\left[2\sigma^{2}_{\text{G}}\right]$
(i.e. unit dB). At $\sigma_{\text{G}}\ll 1$ so
$e^{-\pi\sigma_{\text{G}}^{2}n^{2}}\sim 1$, the Gaussian envelope can be
ignored so that the approximate GKP state approaches the ideal GKP state. The
approximate GKP considers an uncertainty
$\xi^{(\text{G})}_{q(p),2}\in\mathcal{N}\left(0,\sigma^{2}_{\text{G}}\right)$
on both quadratures of each tooth.

Figure 5: General architecture of CV QEC protocol. The light blue shaded area
denotes decoding (i.e. DEC), and the light yellow shaded area denotes encoding
(i.e. ENC).

###  3.2 GKP-two-mode-squeezing code

The CV QEC code that is assisted with GKP state refers to GKP code and were
developed to protect a bosonic mode by encoding it into multiple oscillator
modes. A few of GKP codes have already been well discussed in Ref. [76], such
as, GKP-two-mode-squeezing (GKP-TMS), GKP-repetition (GKP-R) and GKP-
squeezing-repetition (GKP-SR) codes, and, for consistency, the following QEC
protocols all refer to the GKP-TMS code. To exploit the GKP-TMS code in the CV
QR architecture, a QR node that entails an encoding operation and a decoding
operation is designed, as sketched in Fig. 5.

To correct the additive noise, which can be modeled as independent and
identically distributed (i.i.d.) Gaussian random displacements
$\left(\xi_{q,1},\xi_{p,1},\xi_{q,2},\xi_{p,2}\right)$ on the four quadratures
of the two modes, the _encoding_ process is carried out by a two-mode-
squeezing (TMS) gate,
$\hat{T}_{1,2}(g)=e^{\frac{g}{2}\left(\hat{a}_{1}\hat{a}_{2}-\hat{a}^{\dagger}_{1}\hat{a}^{\dagger}_{2}\right)}$,
where $g=\ln\left(\sqrt{G}+\sqrt{G-1}\right)$ is determined by $G$ (i.e.
$G\geq 1$), and $\hat{a}_{1}$ and $\hat{a}_{2}$ denote, respectively, the
annihilation operators of the involved bosonic message and ancilla mode (mode
1 and mode 2). The _decoding_ process entails three parts: inverse TMS
operation (i.e. $\hat{T}^{\dagger}_{1,2}$), estimating the channel-induced
noise by a quantum measurement and using displacement operations
$\hat{D}_{2}\left[-\bar{\xi}_{q,1},0\right]$ and
$\hat{D}_{3}\left[0,-\bar{\xi}_{p,1}\right]$ to compensate for the
displacement errors incurred by the channel on the message mode, where the
displacements depend on the measurement outcomes of the corresponding modes.
To perform the quantum measurement, one introduces an additional GKP ancilla
(mode 3). Two homodyne measurements on the prepared two ancilla modes (mode 2
and mode 3), implemented by a SUM gate beforehand, i.e.,
$\hat{S}_{2,3}=e^{-i\hat{q}_{2}\otimes\hat{p}_{3}}$). Here, $\bar{\xi}_{q,1}$
and $\bar{\xi}_{p,1}$ are the estimations of the displacement error
$\xi_{q,1}$ and $\xi_{p,1}$ acquired by measuring the ancila states in mode 2
and mode 3. The corrected message mode is subsequently encoded with a new GKP
ancilla at mode 2 generated at the present QR node, and both mode 1 and mode 2
are transmitted to the next QR node for decoding and encoding.

The displacement noise continuously accumulates on the message mode until it
arrives at Bob’s terminal. In a weak additive noise regime [77], the
displacement noise is approximately a Gaussian noise so the Wigner function of
the message mode can be fully derived based on the variance of displacement
noise. Let $L$ be the physical distance between Alice and Bob, the average
variances of the displacement noise for Bob’s received message mode are
derived as

$\displaystyle\Sigma^{2}_{QA}=$
$\displaystyle\left(L/L_{\Delta}\right)\Sigma^{2}_{Q}\left[\sigma^{2}_{A}\right]\;,\;\;\Sigma^{2}_{QT}=$
$\displaystyle\left(L/L_{\Delta}\right)\Sigma^{2}_{Q}\left[\sigma^{2}_{T}\right],$
(8)

over, respectively, the QEC amplified one-way and the QEC teleportation
channels, where $\sigma^{2}_{A(T)}$ is a number given by Eq. (1) (Eq. (2)).
Here, $\Sigma_{Q}^{2}$ is a function to calculate the variance of the
displacement noise (see Appendix A for the derivation of $\Sigma_{Q}^{2}$).

###  3.3 Fidelity Performances

This section compares the performances of CV QR with different types of
quantum channels from the choices of pre- and post-processing. We will focus
on the establishment of CV entanglement in the form of TMSV pairs between
Alice and Bob. The overall input-output relations are constructed as the
following channels: $\mathcal{T}^{D}_{L}$ for direct one-way transmission,
$\mathcal{T}^{QA}_{L}$ for QEC amplified one-way transmission, and
$\mathcal{T}^{QT}_{s,L}$ for QEC teleportation. In the three regimes, the GKP-
TMS code is optimized over $G$ for any given parameters of the inter-repeater
spacing $L_{\Delta}$, the squeezing parameter $s$ of the TMSV in quantum
teleportation, and the finite squeezing teeth of the GKP state
$s^{\text{(G)}}$.

$r$ | Message squeezing | $L$ | Spacing between Alice and Bob
---|---|---|---
$s$ | Teleportation squeezing | $L_{\Delta}$ | Inter-repeater spacing
$s^{\text{(G)}}$ | GKP teeth squeezing | $L^{*}_{\Delta}$ | MIRS s.t. $\sigma^{2}_{T}<\sigma^{2}_{A}$
Table 1: Definitions of related symbols.

To establish CV entanglement in the form of TMSV pairs, we focus on the
following scenario: Alice generates a TMSV state consisting of a pair of
modes, signal and idler, characterized by the squeezing level $r$. Alice
attempts to transmit the idler mode to Bob via a series of QRs while locally
retaining the signal mode. In doing so, Alice and Bob share a pair of noisy
TMSV. We will evaluate the performance of the QR in terms of the fidelity of
the established TMSV to the ideal TMSV. The symbols of related parameters are
summarized in Tab. 1.

The Uhlmann fidelity is a measure to quantify the similarity between two
density operators, $\hat{\rho}$ and $\hat{\sigma}$, defined as

$F\left(\hat{\sigma},\hat{\rho}\right)\equiv\left(\text{Tr}\left[\sqrt{\sqrt{\hat{\sigma}}\hat{\rho}\sqrt{\hat{\sigma}}}\right]\right)^{2}.$
(9)

This is used to quantify the deviation between the distributed TMSV state and
the original TMSV state.

The fidelities of direct one-way transmission (i.e. neither pre- nor post-
processing), QEC amplified one-way transmission, and QEC teleportation are
defined, respectively, as

$\displaystyle F_{{O}}\equiv$ $\displaystyle
F(\hat{\rho}_{\text{S},\text{I}},\hat{\rho}^{D}_{\text{S},\text{I}})\;,\;\;F_{A}\equiv$
$\displaystyle
F(\hat{\rho}_{\text{S},\text{I}},\hat{\rho}^{QA}_{\text{S},\text{I}})\;,\;\;F_{T}\equiv$
$\displaystyle
F(\hat{\rho}_{\text{S},\text{I}},\hat{\rho}^{QT}_{\text{S},\text{I}}),$ (10)

where

$\displaystyle\hat{\rho}^{{D}}_{\text{S,I}}=$
$\displaystyle\left(\mathcal{I}\otimes\mathcal{T}^{D}_{L}\right)\left[\hat{\rho}_{\text{S},\text{I}}\right]\;,\;\;\hat{\rho}^{QA}_{\text{S,I}}=$
$\displaystyle\left(\mathcal{I}\otimes\mathcal{T}^{QA}_{L}\right)\left[\hat{\rho}_{\text{S},\text{I}}\right]\;,\;\;\hat{\rho}^{QT}_{\text{S,I}}=$
$\displaystyle\left(\mathcal{I}\otimes\mathcal{T}^{QT}_{s,L}\right)\left[\hat{\rho}_{\text{S},\text{I}}\right].$
(11)

Here, $\mathcal{I}$ is the identity channel assuming ideal signal storage, and
$\hat{\rho}_{S,I}=\ket{\text{TMSV}}\bra{\text{TMSV}}_{\text{S,I}}$ is the
input TMSV state.

Figure 6: Fidelities of direct one-way transmission ($F_{O}$), QEC amplified
one-way transmission ($F_{A}$), and $m$-relay QEC teleportation ($F_{T}$)
based on ideal GKP state with repeater separation (a) $L_{\Delta}=0.25$ km and
(b) $L_{\Delta}=1$ km. $s^{(G)}\rightarrow\infty$, $r=15$ dB and
$s=\left\\{20,\infty\right\\}$ dB (i.e.
$L^{*}_{\Delta}=\left\\{0.44,0\right\\}$ km).

First, let us assume a perfect GKP state is available (i.e.
$s^{(\text{G})}\rightarrow\infty$) and plot the optimized fidelities in Fig.
6(a) and (b). Given that the teleportation squeezing is $s=20$ dB, we choose
$L_{\Delta}=250$ m to coincide with the optimal repeater separation that
Rozpȩdek _et al._ selected in their article [73].

The simulation result indicates that at an infinite teleportation squeezing
level, i.e. $s\rightarrow\infty$, $\sigma^{2}_{A}>\sigma^{2}_{T}$ always
holds, yielding $L^{*}_{\Delta}=0$; yet, infinite squeezing requires unbounded
energy and is therefore unphysical. With a practical finite teleportation
squeezing level, there is an associated non-zero MIRS. However, a shorter
inter-repeater spacing increases the density of QRs and the associated
resource overhead. In contrast, the QR protocol based on quantum teleportation
channels reduces the density of QRs while maintaining a high fidelity for the
transmitted quantum states by placing the TMSV source in the middle between
two QR nodes separated by a distance of $L_{\Delta}>L^{*}_{\Delta}$, as shown
in Fig. 6(b). The GKP-TMS code drastically improves the fidelity for the
transmitted quantum state in both channel scenarios, as compared to the direct
one-way transmission.

Assuming using imperfect GKP states in QEC, $F_{T}$s are plotted in Fig.
7(a)(b)(c)(d) as functions of $L$ and $s^{(\text{G})}$ while fixing $r=15$ dB,
corresponding to different $s$. Fig. 7 concludes that $s^{(\text{G})}\gtrsim
s\gtrsim r$ is required for effective QEC over quantum teleportation channels;
otherwise, under $s<r$, the additive noise caused by teleportation will add
too much noise to the transmitted quantum state while under
$s^{(\text{G})}<s$, the GKP state only increases the added noise because the
variance of GKP state is even larger than the noise to be corrected.

Figure 7: Fidelities of teleportation-enabled repeater ($F_{T}$) based on
imperfect GKP. Input TMSV is set $r=15$ dB and the inter-repeater separation
is $L_{\Delta}=1$ km. (a) $s=10$ dB ($L^{*}_{\Delta}=4.6$ km), (b) $s=15$ dB
($L^{*}_{\Delta}=1.4$ km), (c) $s=20$ dB ($L^{*}_{\Delta}=0.44$ km), (d)
$s=25$ dB ($L^{*}_{\Delta}=0.14$ km).

###  3.4 Concatenation of GKP-TMS code

Figure 8: Scheme of $m$-relay $k$-layer CV QEC repeaters.

Recent study has proven that concatenation of multiple layers of QEC would
substantially reduce the displacement noise comparing with only a single layer
code [77]. In a multi-layer QEC scheme, Alice, Bob and all repeaters prepare
$k$ GKP ancilla (i.e. $k\in\mathbb{N}$) to be encoded with a single message
state, shown in Fig. 8 and are used to decode teleported state at Bob’s
terminal. As the assumption before, the physical noise of single layer
encoding can be approximately Gaussian. This assumption can be justified if
the displacement noise is much less than unity [77].

Although the resources for implementing $m$-relay $k$-layer GKP-assisted QEC
are immense (i.e. in total, $2\left(m+1\right)k$ GKP ancilla modes need to be
prepared beforehand), the correction outcomes are remunerable. In Fig. 9, we
demonstrate the fidelities of the $m$-relay QEC QRs, that correspond to
different layers of QEC and it shows that the fidelities are significantly
improved. Albeit TMSV and GKP modes are never ideal in practice, concatenating
multi-layer QEC codes is an alternative approach to suppress the additive
noises of the channel, shown in Fig. 9. In Fig. 9, as $k\gtrsim 13$, the
endmost iterative noise almost converges to a finite value, which is
ultimately determined by $s^{(\text{G})}$.

Figure 9: Fidelities of $m$-relay $k$-layer ($k=1,2,13$) QEC teleportation.
The message squeezing is $r=15$ dB.

##  4 Applications

Preshared entanglement between distant parties underpins numerous quantum
applications. Nonetheless, establishing entanglement at a distance is impeded
by the loss of the entanglement-distribution channel. The proposed GKP-
assisted QEC can correct the Gaussian errors to enhanced the performance of a
multitude of applications, including EA communication, QI, and CV QKD. For
simplicity, we will set $s^{(\text{G})}=s$ in the following performance
analysis on the three applications assisted by the proposal QR protocol.

###  4.1 Entanglement-assisted communication

Figure 10: Entanglement-assisted communication (a) scheme and (b) the phase
encoding Holevo information normalized to classical capacity ($H/C$) at $L=25$
km for the cases of ideal, direct one-way transmission and $m$-relay $k$-layer
($k=1,13$) QEC teleportation.

The classical information rate over a thermal-loss channel is upper bounded by
the classical capacity [47,83], formulated as

$C=g\left(\kappa N_{S}+N_{B}\right)-g\left(N_{B}\right),$ (12)

where
$g\left(x\right)\equiv\left(x+1\right)\log_{2}\left(x+1\right)-x\log_{2}x$,
$N_{S}$ is the mean photon number of a signal mode, $\kappa$ is the
transmissivity of the channel, and $N_{B}$ is the mean photon number of
thermal-noise bath mode. EA communication is able to surpass the classical
capacity [47,84-86]. In an ideal EA communication scenario illustrated in Fig.
10(a), Alice performs phase encoding on the signal mode of a preshared TMSV
state and sends it to Bob over a very lossy and noisy channel, i.e.,
$\kappa\ll 1$ and $N_{B}\gg 1$. Bob then performs a joint measurement on the
received signal with the idler at hand.

However, building up preshared entanglement in real-world operational
environments hinges on lossy entanglement-distribution channels that degrade
the quality of the entanglement, holding back the advantage endowed by EA
communication. The proposed CV QR architecture opens a promising route toward
mitigating the loss arising from the entanglement-distribution channel.

The EA capacity normalized to the classical capacities are sketched as the
dashed black, dashed blue, solid blue, solid purple and solid magenta curves,
associated with different scenarios of entanglement sharing, in Fig. 10(b).
Over an extremely lossy and noisy communication channel, the asymptotic Holevo
capacity normalized to the classical capacity is given by

$\displaystyle H_{\text{Ideal}}/C$
$\displaystyle\approx\left(N_{S}+1\right)\log\left(1+1/N_{S}\right)\,\,,\,\,\,H_{D}/C\approx\eta\left(N_{S}+1\right)\log\left(1+1/\eta
N_{S}\right),$ (13) $\displaystyle H_{QT}/C$
$\displaystyle\approx\left(N_{S}+1\right)\log\left(1+1/\Sigma^{2}_{QT}\right)-N_{S}/\left(\Sigma^{2}_{QT}+\Sigma^{4}_{QT}\right),$

where $H_{\text{Ideal}}$, $H_{D}$ and $H_{QT}$ denote the Holevo information
associated with ideal preshared TMSV states, TMSV sharing via direct one-way
transmission and QEC teleportation-enabled QR. The QEC inevitably introduces
thermal noise, causing the EA Holevo information to saturate at weak
$N_{S}$’s. In this regime, teleportation is inferior to direct one-way
transmission in entanglement distribution. Conversely, as $N_{S}$ increases,
QEC teleportation-enabled QR starts to outperform the direct one-way
entanglement distribution approach. Under this parameter setting, we find that
the multi-layer encoding on finite squeezed TMSV and GKP states is more
powerful than single-layer encoding on infinitely squeezed TMSV and GKP
states.

###  4.2 Quantum illumination

Figure 11: (a) Scheme of quantum illumination. (b) The quantum Chernoff bounds
of error probability $P_{e}$ versus transmitted $M$ modes for CI and three QI
cases: ideal entanglement distribution, direct one-way, and $m$-relay
$k$-layer ($k=1,13$) QEC teleportation at $L=25$ km with $N_{S}=0.01$.

QI is a paradigm for quantum-enhanced target detection through a very lossy
and noisy environment [37,40,87-89]. Illustrated in Fig. 11(a), the QI
transmitter prepares TMSV states composed of entangled signal-idler mode
pairs. The idler modes are distributed to receiver over a distribution channel
while the signal modes are transmitted to interrogate a target residing in an
environment modeled as a thermal-loss channel. The QI receiver performs a
joint measurement on the transmitted signal embedded in a bright noise
background and the idler to infer the presence or absence of the target. Tan
et al. [37] showed that QI, with ideal equipment and the optimum quantum
receiver, achieves a 6-dB advantage in the error-probability exponent of the
quantum Chernoff bound (magenta curve in Fig. 11(b)) over that of classical
illumination (CI) based on the coherent-state transmitter and homodyne
receiver (cyan curve in Fig. 11(b)).

A practical challenge for QI lies in the requirement for high-fidelity quantum
memories used to match the propagation delay between the signal and idler
modes. At present, QI experiments [90] utilize low-loss optical fibers to
store the idler, which mimics the one-way entanglement distribution channel.
Due to the idler-storage loss, QI’s advantage over CI quickly diminishes, as
shown in the black dashed curve of Fig. 11(b). The proposed QR architecture
based on QEC and teleportation would constitute an effective approach to
mitigate the idler-storage loss. The blue dashed and solid curves in Fig.
11(b) depicts the simulation results for QI enhanced by QEC on the idler
modes, showing reduced error probabilities as compared to QI without QEC. Akin
to EA communication, in this case the multi-layer QEC with finite squeezing
outperforms the single-layer QEC with infinite squeezing.

###  4.3 CV quantum key distribution

CV-QKD enables two distant parties, Alice and Bob, to securely share a common
binary random key despite the adversary, Eve, mounts the optimal attack to
capture the communicated information [91-93]. Unlike its DV counterpart, CV-
QKD can be fully implemented with off-the-shelf telecommunication components
without resorting to single-photon detectors and is thus particularly
intriguing for real-world deployment. The security of CV-QKD protocols is
analyzed by upper bounding the accessible information to Eve assuming her
power is only limited by the laws of physics. Specifically, the secret-key
rate (SKR) for CV-QKD is given by

$R=-\frac{1}{2}\log_{2}{\left[e^{2}\epsilon\left(1+\epsilon\right)/4\right]},$
(14)

where $\epsilon$ quantifies the variance of overall additive excess noise. The
proposed QR architecture based on $m$-relay $k$-layer QEC mitigates the loss
of the quantum channel to boost the SKR, as shown in Fig. 12.

To further investigate the application of the QR architecture to CV-QKD, two
additional remarks on Fig. 12 are worth making. First, the SKR of the QR
architecture based on $k$-layer QEC and teleportation are below the PLOB bound
at $k\leq 9$, hindered by the accumulated noise introduced at the QR nodes.
Second, given $s=s^{(\text{G})}=25$ dB, the theoretical maximal distance of
the QR architecture based on 13-layer QEC and teleportation, as shown in the
purple curve of Fig. 12, reaches 596 km. We expect that the incorporation of
an additional DV QEC layer would suppress the residue noise and further extend
the CV-QKD distance [73,74].

Figure 12: The secret key rates per mode of $m$-relay $k$-layer
($k=1,7,10,13$) QEC teleportation-enabled QR.

##  5 Discussion and outlook

The QR architecture based on teleportation channels places an entanglement
source in the middle of two adjacent QR nodes. In contrast, the QR scheme
based on amplified one-way channels directly connects the adjacent nodes by
optical fibers. One may argue that adding an intermediate QR node in a
amplified one-way channel would surpass the performance of the teleportation-
based scheme. However, a full-scale QR node needs multiple GKP ancilla modes,
which consumes much more resources than the widely available TMSV source.

The combination of CV and DV QEC was recently proposed by Rozpȩdek _et al._
[73] and Fukui _et al._ [74]. Such a hybrid QEC scheme would allow the
proposed QR architecture based on $m$-relay $k$-layer QEC to be further
concatenated with a DV QEC code to drastically reduce the amount of residue
noise. As long as the CV errors after QEC are limited to a certain range, the
DV QEC will be able to correct these errors to maximize the fidelity of the
transmitted the quantum state.

##  6 Conclusions

In this article, we proposed a deterministic CV QR architecture based on GKP
states to enable the distribution of CV quantum states over long distances. It
is shown that the proposed QR architecture based on GKP QEC can significantly
reduce the additive errors caused by a lossy bosonic channel. In our study, we
showed that the QR architecture based on GKR QEC and teleportation outperforms
direct one-way transmission when the squeezing level is higher than 15 dB. The
proposed QR architecture is applied to improve the performance of EA
communication, QI and CV-QKD. Once GKP state with sufficient squeezing becomes
available, the proposed QR architecture will enable CV quantum states to be
faithfully transmitted over unprecedented distances, thereby making a large
stride forward in the development of quantum technology.

This research is supported by National Science Foundation Grant No.
ECCS-1920742, CCF-1907918, and No. EEC-1941583, and Defense Advanced Research
Projects Agency (DARPA) under Young Faculty Award (YFA) Grant No.
N660012014029.

##  Appendix A Full derivation of $\Sigma^{2}_{Q}$

The variances of both quadratures are identical and can be calculated through
the integration formula,

$\displaystyle\Sigma^{2}_{Q}\left[\sigma^{2}\right]=$
$\displaystyle\sum^{\infty}_{n=-\infty}\int_{-\infty}^{\infty}d\xi^{(\text{G})}_{q,2}\int_{-\infty}^{\infty}d\xi_{q,2}\int_{-\infty}^{\infty}d\xi_{q,1}\;\left[\frac{1}{\sqrt{4\pi}\sigma_{G}}e^{-\frac{\xi_{q,2}^{(\text{G})2}}{4\sigma^{2}_{G}}}\right]\left[\frac{1}{2\pi\sigma^{2}}e^{-\frac{2G-1}{2\sigma^{2}}\left(\xi^{2}_{q,1}+\xi^{2}_{q,2}\right)-\frac{2\sqrt{G(G-1)}\xi_{q,1}\xi_{q,2}}{\sigma^{2}}}\right]\times$
(15)
$\displaystyle\left(\xi_{q,1}-\bar{\xi}_{q,1}\right)^{2}\mathcal{U}\left(\xi_{q,2}+\xi^{(\text{G})}_{q,2}\in\left[\left(n-\frac{1}{2}\right)\sqrt{2\pi},\left(n+\frac{1}{2}\right)\sqrt{2\pi}\right]\right),$

where $\mathcal{U}$ is an indicator function (i.e.
$\mathcal{U}\left(\mathcal{S}\right)=1$, if $\mathcal{S}$ is true; otherwise,
$\mathcal{U}\left(\mathcal{S}\right)=0$), $\bar{\xi}_{q,1}$ is the estimation
of $\xi_{q,1}$. We infer $\bar{\xi}_{q,1}=c_{q}\,\xi_{q,2}$ and specify
$c_{q}$ by the minimum variance unbiased estimation (MVUE) [76]; hence,

$\bar{\xi}_{q,1}=-\frac{2\sqrt{G\left(G-1\right)}\sigma^{2}}{\left(2G-1\right)\sigma^{2}+2\sigma_{\text{G}}^{2}}R_{\sqrt{2\pi}}\left[\xi_{q,2}+\xi^{(\text{G})}_{q,2}\right],$
(16)

and, Eq. (15) becomes

$\displaystyle\Sigma^{2}_{Q}\left[\sigma^{2}\right]=$
$\displaystyle\sum^{\infty}_{n=-\infty}\int_{-\infty}^{\infty}d\xi^{(\text{G})}_{q,2}\int_{-\infty}^{\infty}d\xi_{q,2}e^{-\frac{\xi^{2}_{q,2}}{2\left(2G-1\right)\sigma^{2}}-\frac{\xi_{q,2}^{(\text{G})2}}{4\sigma^{2}_{G}}}\;\mathcal{U}\left(\xi_{q,2}+\xi^{(\text{G})}_{q,2}\in\left[\left(n-\frac{1}{2}\right)\sqrt{2\pi},\left(n+\frac{1}{2}\right)\sqrt{2\pi}\right]\right)$
(17)
$\displaystyle\times\left\\{\frac{\sigma}{\left[2\left(2G-1\right)\right]^{3/2}\pi\sigma_{G}}+\frac{\sqrt{2}G\left(G-1\right)\left[\left(2G-1\right)\left(n\sqrt{2\pi}-\xi^{(\text{G})}_{q,2}\right)\sigma^{2}+2\xi_{q,2}\sigma^{2}_{G}\right]^{2}}{\left(2G-1\right)^{5/2}\pi\sigma_{G}\left[\left(2G-1\right)\sigma^{2}+2\sigma^{2}_{G}\right]^{2}\sigma}\right\\}$
$\displaystyle=$
$\displaystyle\sum^{\infty}_{n=-\infty}\left\\{\frac{\sigma^{2}\left[8\left(G-1\right)Gn^{2}\pi\sigma^{2}+\left(2G-1\right)\sigma^{4}+4\left(1+2G\left(G-1\right)\right)\sigma^{2}\sigma^{2}_{G}+4\left(2G-1\right)\sigma^{4}_{G}\right]}{2\left[\left(2G-1\right)\sigma^{2}+2\sigma^{2}_{G}\right]^{2}}\right\\}$
$\displaystyle\times\left\\{\text{erfc}\left[\frac{\left(n-1/2\right)\sqrt{\pi}}{\sqrt{\left(2G-1\right)\sigma^{2}+2\sigma^{2}_{G}}}\right]-\text{erfc}\left[\frac{\left(n+1/2\right)\sqrt{\pi}}{\sqrt{\left(2G-1\right)\sigma^{2}+2\sigma^{2}_{G}}}\right]\right\\}.$

## References

## References

  * [1] H. J. Kimble, _The quantum internet_ , Nature 453, 1023(2008).
  * [2] J. Biamonte, M. Faccin, and M. De Domenico, _Complex networks from classical to quantum_ , Commun. Phys. 2, 53 (2019).
  * [3] S. Wehner, D. Elkouss, and R. Hanson, _Quantum internet:A vision for the road ahead_ , Science 362, (2018).
  * [4] W. Kozlowski and S. Wehner, Towards large-scale quan-tum networks, in _Proceedings of the Sixth Annual ACM International Conference on Nanoscale Computing and Communication_ pp. 1–7 (2019).
  * [5] J. Miguel-Ramiro, A. Pirker, and W. Dür, _Genuine quantum networks: superposed tasks and addressing_ , arXiv:2005.00020 (2020).
  * [6] C. H. Bennett and G. Brassard, _Quantum cryptography: Public key distribution and coin tossing_ , Theoretical Computer Science 560, 7 (2014).
  * [7] N. Gisin, G. Ribordy, W. Tittel, and H. Zbinden, _Quantum cryptography_ , Rev. Mod. Phys. 74, 145 (2002).
  * [8] N. Gisin and R. Thew, _Quantum communication_ , Nat. Photonics 1, 165 (2007).
  * [9] W. K. Wootters and W. H. Zurek, _A single quantum cannot be cloned_ , Nature 299, 802 (1982).
  * [10] M. Koashi and N. Imoto, _No-cloning theorem of entangled states_ , Phys. Rev. Lett 81, 4264 (1998).
  * [11] S. Pirandola, R. Laurenza, C. Ottaviani, and L. Banchi, _Fundamental limits of repeaterless quantum communications_ , Nat. Commun. 8, 1 (2017).
  * [12] S. Pirandola, _End-to-end capacities of a quantum communication network_ , Commun. Phys 2, 51 (2019).
  * [13] W. Dür, H.-J. Briegel, J. I. Cirac, and P. Z, _ler, Quan-tum repeaters based on entanglement purifica_ , Phys. Rev. A 59, 169 (1999).
  * [14] P. Van Loock, T. Ladd, K. Sanaka, F. Yamaguchi, K. Nemoto, W. Munro, and Y. Yamam, _Hybrid quantum repeater using bright coherent light_ , Phys. Rev. Lett. 96, 240501 (2006).
  * [15] L. Jiang, J. M. Taylor, K. Nemoto, W. J. Munro, R. V. Meter, and M. D. Lukin, _Quantum repeater with encoding_ , Phys. Rev. A 79, 032325 (2009).
  * [16] S. Santra, L. Jiang, and V. S. Malinov, _Quantum repeater architecture with hierarchically optimized memory buffer time_ , Quantum Science and Technology 4, 025010 (2019).
  * [17] Z.-D. Li, R. Zhang, X.-F. Yin, L.-Z. Liu, Y. Hu, Y.-Q.Fang, Y.-Y. Fei, X. Jiang, J. Zhang, L. Li, et al., _Experimental quantum repeater without quantum memory_ , Nat. Photonics 13, 644 (2019).
  * [18] H. J. Briegel, W. Dür, J. I. Cirac, and P. Zoller, _Quantum repeaters: The role of imperfect local operations in quantum communication_ , Phys. Rev. Lett. 81, 5932 (1998).
  * [19] W. Dür J., H. J. Briegel, I. Cirac, and P. Zoller, _Quantum repeaters based on entanglement purification_ , Phys. Rev. A 59, 169 (1999).
  * [20] P. Hilaire, E. Barnes, and S. E. Economou, _Resource requirements for efficient quantum communication using all-photonic graph states generated from a few matter qubits_ , Quantum 5, 397 (2021).
  * [21] W. Munro, K. Harrison, A. Stephens, S. Devitt, and K. Nemoto, _From quantum multiplexing to high-performance quantum network_ , Nucl. Phys. 4, 792 (2010).
  * [22] T. C. Ralph, A. J. F. Hayes, and A. Gilchrist, _Loss-tolerantoptical qubit_ , Phys. Rev. Lett. 95, 100501 (2005).
  * [23] P. Kok, W. J. Munro, K. Nemoto, T. C. Ralph, J. P. Dowl-ing, _, Linear optical quantum computingwith photonic qubi_ , Rev. Mod. Phys. 79, 135 (2007).
  * [24] S. Muralidharan, J. Kim, N. Lütkenhaus, M. D. Lukin, and L. Jiang, _Ultrafast and Fault-Tolerant Quantum Communication across Long Distances_ , Phys. Rev. Lett. 112, 250501 (2014).
  * [25] S. Muralidharan, C.-L. Zou, L. Li, J. Wen, and L. Jiang, _Overcoming erasure errors with multilevel systems_ , New J. Phys. 19, 013026 (2017).
  * [26] S. Muralidharan, C.-L. Zou, L. Li, and L. Jiang, _One-way quantum repeaters with quantum Reed-Solomon codes_ , Phys. Rev. A 97, 052316 (2018).
  * [27] W. Ge, K. Jacobs, Z. Eldredge, A. V. Gorshkov, and M. Foss-Fei, _One-way quantum repeaters with quantum Reed-Solomon codes_ , Phys. Rev. A 97, 052316 (2018).
  * [28] T. J. Proctor, P. A. Knott, and J. A. Dunningham, _Multiparameter estimation in networked quantum sensors_ , Phys. Rev. Lett. 120, 080501 (2018).
  * [29] Q. Zhuang, Z. Zhang, and J. H. Shapiro, _Distributed quantum sensing using continuous-variable multipartite entanglement_ , Phys. Rev. A 97, 032329 (2018).
  * [30] Z. Eldredge, M. Foss-Feig, J. A. Gross, S. L. Rolston,and A. V. Gorshkov, _Optimal and secure measurement protocols for quantum sensor networks_ , Phys. Rev. A 97, 042337 (2018).
  * [31] X. Guo, C. R. Breum, J. Borregaard, S. Izumi, M. V.Larsen, T. Gehring, M. Christandl, J. S. Neergaard-Nielsen, and U. L. Andersen, _Distributed quantum sensing in a continuous-variable entangled network_ , Nat. Phys. 16, 281 (2020).
  * [32] Y. Xia, W. Li, W. Clark, D. Hart, Q. Zhuang, and Z. Zhang, _Demonstration of a reconfigurable entangled radio-frequency photonic sensor network_ , Phys. Rev. Lett. 124, 150502 (2020).
  * [33] Z. Zhang and Q. Zhuang, _Distributed quantum sensing_ , Quantum Science and Technology 6, 043001 (2021).
  * [34] K. Qian, Z. Eldredge, W. Ge, G. Pagano, C. Monroe, J. V.Porto, and A. V. Gorshkov, _Heisenberg-scaling measurement protocol for analytic functions with quantum sensor networks_ , Phys. Rev. A 100, 042304 (2019).
  * [35] Q. Zhuang and Z. Zhang, _Physical-layer supervised learn-ing assisted by an entangled sensor network_ , Phys. Rev. X 9, 041023 (2019).
  * [36] Y. Xia, W. Li, Q. Zhuang, and Z. Zhang, _Quantum-enhanced data classification with a variational entangled sensor network_ , Phys. Rev. X 11, 021047 (2021).
  * [37] Si-Hui Tan et. al, _Quantum Illumination with Gaussian States_ , Phys. Rev. Lett. 101, 253601 (2008).
  * [38] Z. Zhang, S. Mouradian, F. N. Wong, and J. H. Shapiro, _Entanglement-enhanced sensing in a lossy and noisy environment_ , Phys. Rev. Lett. 114, 110506 (2015).
  * [39] Q. Zhuang, Z. Zhang, and J. H. Shapiro, _Optimum mixed-state discrimination for noisy entanglement-enhanced sensing_ , Phys. Rev. Lett. 118, 040801 (2017).
  * [40] J. H. Shapiro, _The Quantum Illumination Story_ , IEEE Aerospace and Electronic System Magazine 35, 8 (2020).
  * [41] Q. Zhuang, _Quantum ranging with gaussian entanglement_ , Phys. Rev. Lett. 126, 240501 (2021).
  * [42] C. H. Bennett and S. J. Wiesner, _Communication via one-and two-particle operators on einstein-podolsky-rosen states_ , Phys. Rev. Lett. 69, 2881 (1992).
  * [43] C. Bennett, P. Shor, J. Smolin, and A. Thapliyal, _Entanglement-assisted capacity of a quantum channel and the reverse shannon theorem_ , IEEE Trans. Inf. Theory 48, 2637 (2002).
  * [44] C. H. Bennett, P. W. Shor, J. A. Smolin, and A. V. Thap-liyal, _Entanglement-assisted classical capacity of noisy quantum channels_ , Phys. Rev. Lett. 83, 3081 (1999).
  * [45] A. S. Holevo, _On entanglement-assisted classical capacity_ , J. Math. Phys. 43, 4326 (2002).
  * [46] M.-H. Hsieh, I. Devetak, and A. Winter, _Entanglement-assisted capacity of quantum multiple-access channels_ , IEEE. Trans. Inf. Theory 54, 3078 (2008).
  * [47] H. Shi, Z. Zhang, and Q. Zhuang, _Practical route to entanglement-assisted communication over noisy bosonic channels_ , Phys. Rev. Applied 13, 034029 (2020).
  * [48] Q. Zhuang, _Quantum-enabled communication without a phase reference_ , Phys. Rev. Lett. 126, 060502 (2021).
  * [49] H. Shi, M.-H. Hsieh, S. Guha, Z. Zhang, and Q. Zhuang, _Entanglement-assisted multiple-access channels: capacity regions and protocol designs_ , arXiv:2101.12173 (2021).
  * [50] S. Hao, H. Shi, W. Li, J. H. Shapiro, _Entanglement-assisted communication surpassing the ultimate classical capacity_ , Phys. Rev. Lett. 126, 250501 (2021).
  * [51] K. P. Seshadreesan, H. Krovi, and S. Guha, _Continuous-variable quantum repeater based on quantum scissors and mode multiplexing_ , Phys. Rev. Research 2, 013310 (2020).
  * [52] T. C. Ralph and A. P. Lund, _Nondeterministic noiseless linear amplification of quantum systems_ , AIP Conference Proceedings. 1110, 155 (2009).
  * [53] P. W. Shor, _Fault-Tolerant Quantum Computation_ , arXiv:quant-ph/9605011 (1997).
  * [54] J. Preskill, _Fault-Tolerant Quantum Computation_ , arXiv:9712048v1 (1997).
  * [55] B. M. Terhal, J. Conrad, and C. Vuillot, _Towards scalable bosonic quantum error correction_ , Quantum Science and Technology 5, 043001 (2020).
  * [56] M. H. Michael, M. Silveri, R. T. Brierley, V. V. Albert, J. Salmilehto, L. Jiang, and S. M. Girvin, _Error-Correcting Codes for A Bosonic mode_ , Phys. Rev. X 6, 031006 (2016).
  * [57] P. T. Cochrane, G. J. Milburn, and W. J. Munro, _Macroscopically distinct quantum-superposition states as a bosonic code for amplitude damping_ , Phys. Rev. A 59, 2631 (1999).
  * [58] D. Gottesman, A. Kitaev, and J. Preskill, _Encoding a qubit in an oscillator_ , Phys. Rev. A 64, 012310 (2001).
  * [59] J. Guillaud and M. Mirrahimi, _Repetition Cat-Qubits: Fault-Tolerant Quantum Computation with Highly Reduced Overhead_ , Phys. Rev. X 9, 041053 (2019).
  * [60] S. Puri, L. St-Jean, J. A. Gross, A. Grimm, N. E. Frattini,P. S. Iyer, A. Krishna, S. Touzard, L. Jiang, A. Blais, S. T. Flammia, and S. M. Girvin, _Bias-preserving gates with stabilized cat qubits_ , Science Advances 6, 10.1126/sci-adv.aay5901 (2020).
  * [61] K. Noh, and C. Chamberland, _Fault-tolerant bosonic quantum error correction with the surface–Gottesman-Kitaev-Preskill code_ , Phys. Rev. A 101, 012316 (2020).
  * [62] V. V. Albert, _Performance and structure of single-mode bosonic codes_ , Phys. Rev. Lett. 97, 032346 (2018).
  * [63] Q. Zhuang, J. Preskill, and L. Jiang, _Distributed quantumsensing enhanced by continuous-variable error correction_ , New J. Phys. 22, 022001 (2020).
  * [64] Z. Leghtas et al., _Confining the state of light to a quantum manifold by engineered two-photon loss_ , Science 347, 853 (2015).
  * [65] N. Ofek et al., _Extending the lifetime of a quantum bit with error correction in superconducting circuits_ , Nature (London) 536, 441 (2016).
  * [66] S. Touzard et al., _Coherent Oscillations Inside a Quantum Manifold Stabilized by Dissipation_ , Phys. Rev. X 8, 021005 (2018).
  * [67] L. Hu et al., _Quantum error correction and universal gate set operation on a binomial bosonic logical qubit_ , Nat. Phys. 15, 503 (2019).
  * [68] P. Campagne-Ibarcq et al., _A stabilized logical quantum bit encoded in grid states of a superconducting cavity_ , arXiv:1907.12487 (2020).
  * [69] A. Grimm et al., _The Kerr-cat qubit: Stabilization, readout, and gates_ , Nature (London) 584, 205 (2020).
  * [70] C. Flühmann, V. Negnevitsky, M. Marinelli, and J. P. Home, _Sequential Modular Position and Momentum Measurements of a Trapped ion Mechanical Oscillator_ , Phys. Rev. X 8, 021001 (2018).
  * [71] C. Flühmann, T. L. Nguyen, M. Marinelli, V. Negnevitsky, K. Mehta, and J. P. Home, _Encoding a qubit in a trapped-ion mechanical oscillator_ , Nature (London) 566, 513 (2019).
  * [72] C. Flühmann and J. P. Home, _Direct characteristic-function tomography of quantum states of the trapped-ion motional oscillator_ , Phys. Rev. Lett. 125, 043602 (2020).
  * [73] F. Rozpȩdek, K. Noh, Q. Xu, S. Guha, and L. Jiang, _Quantum repeaters based on concatenated bosonic and discrete-variable quantum codes_ , npj Quantum Inf. 7, 102 (2021).
  * [74] K. Fukui, R. N. Alexander, and P. v. Loock, _All-Optical Long-Distance Quantum Communication with Gottesman-Kitaev-Preskill qubits_ , Phys. Rev. Research 3, 033118 (2021).
  * [75] K. Fukui and N. C. Menicucci, _An efficient, concatenated, bosonic code for additive Gaussian noise_ , arXiv:2102.01374v1 (2021).
  * [76] K. Noh, S. M. Girvin, and L. Jiang, _Encoding an Oscillator into Many Oscillators_ , Phys. Rev. Lett. 125, 080503 (2020).
  * [77] J. Wu and Q. Zhuang, _Continuous-variable error correction for general Gaussian noises_ , Phys. Rev. Appl. 15, 034073 (2021).
  * [78] S. L. Braunstein and H. J. Kimble, _Teleportation of Continuous Quantum Variables_ , Phys. Rev. Lett. 80, 869 (1998).
  * [79] S. Pirandola and M. Stefano, _Quantum teleportation with continuous variables: A survey_ , Laser Physics 16, 1418 (2006).
  * [80] J. Wu, C. Cui, L. Fan, and Q. Zhuang, _Deterministic microwave-optical transduction based on quantum teleportation_ , arXiv preprint arXiv:2106.14037 (2021).
  * [81] M. Eaton, R. Nehra and O. Pfister, _Non-Gaussian and Gottesman–Kitaev–Preskill state preparation by photon catalysis_ , New J. Phys. 21, 113034 (2019).
  * [82] Y. Shi, C. Chamberland, and A. Cross, _Fault-tolerant preparation of approximate GKP states_ , New J. Phys. 21, 093007 (2019).
  * [83] V. Giovannetti, R. García-Patròn, N. J. Cerf, and A. S. Holevo, _Ultimate classical communication rates of quantum optical channels_ , Nat. Photonics 8, 796 (2014).
  * [84] C. H. Bennett, P. W. Shor, J. A. Smolin, and A. V. Thapliyal, _Entanglement-assisted capacity of a quantum channel and the reverse shannon theorem_ , IEEE Trans. Inf. Theory 48, 2637 (2002).
  * [85] B. Schumacher and M. D. Westmoreland, _Sending classical information via noisy quantum channels_ , Phys. Rev. A 56, 131 (1997).
  * [86] A. S. Holevo, _The capacity of the quantum channel with general signal states_ , IEEE Trans. Inf. Theory 44, 269 (1998).
  * [87] Seth Lloyd, _Enhanced Sensitivity of Photodetection via Quantum Illumination_ , Science 321, 1463 (2008).
  * [88] S. Pirandola and S. Lloyd, _Computable bounds for the discrimination of Gaussian states_ , Phys. Rev. A 78, 012331 (2008).
  * [89] A. Karsa, G. Spedalieri, Q. Zhuang, and S. Pirandola, _Quantum illumination with a generic Gaussian source_ , Phys. Rev. Research. 2, 023414 (2020).
  * [90] Z. Zhang, S. Mouradian, F. N. C. Wong, and J. H. Shapiro, _Entanglement-enhanced sensing in a lossy and noisy environment_ , Phys. Rev. Lett. 114, 110506 (2015).
  * [91] J. Lodewyck et. al, _Quantum key distribution over 25 km with an all-fiber continuous-variable system_ , Phys. Rev. A 76, 042305 (2007).
  * [92] F. Grosshans et al., _Quantum key distribution using gaussian-modulated coherent states_ , Nature 421, 238 (2003).
  * [93] F. Grosshans and P. Grangier, _Continuous Variable Quantum Cryptography Using Coherent States_ , Phys. Rev. Lett. 88, 057902 (2002).


11institutetext: University of Science and Technology of China 11email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
# A Novel Skip Orthogonal List for Dynamic Optimal Transport Problem

Xiaoyang Xu  Hu Ding

###### Abstract

Optimal transportation is a fundamental topic that has attracted a great
amount of attention from machine learning community in the past decades. In
this paper, we consider an interesting discrete dynamic optimal transport
problem: can we efficiently update the optimal transport plan when the weights
or the locations of the data points change? This problem is naturally
motivated by several applications in machine learning. For example, we often
need to compute the optimal transportation cost between two different data
sets; if some change happens to a few data points, should we re-compute the
high complexity cost function or update the cost by some efficient dynamic
data structure? We are aware that several dynamic maximum flow algorithms have
been proposed before, however, the research on dynamic minimum cost flow
problem is still quite limited, to the best of our knowledge. We propose a
novel 2D Skip Orthogonal List together with some dynamic tree techniques.
Although our algorithm is based on the conventional simplex method, it can
efficiently complete each pivoting operation within $O(\absolutevalue{V})$
time with high probability where $V$ is the set of all supply and demand
nodes. Since dynamic modifications typically do not introduce significant
changes, our algorithm requires only a few simplex iterations in practice. So
our algorithm is more efficient than re-computing the optimal transportation
cost that needs at least one traversal over all the
$O(\absolutevalue{E})=O(\absolutevalue{V}^{2})$ variables in general cases.
Our experiments demonstrate that our algorithm significantly outperforms
existing algorithms in the dynamic scenarios.

##  1 Introduction

The discrete optimal transport problem involves finding the optimal transport
plan “$X$” that minimizes the total cost of transporting one weighted dataset
$A$ to another $B$, given a cost function “$C$” [25]. The datasets $A$ and $B$
represent the supply and demand node sets, respectively, and the problem can
be represented as a minimum cost flow problem by adding edges between $A$ and
$B$ to create a complete bipartite graph. The discrete optimal transport
problem finds numerous applications in areas such as image registration [15],
seismic tomography [22], and machine learning [34]. However, these
applications traditionally consider a static scenario where the weights of the
datasets and the cost function remain constant. Yet, many real-world
applications need to handle dynamic optimal transport problems:

  * •

Dataset Similarity. In data analysis, measuring the similarity between
datasets is a crucial task, and optimal transport has emerged as a powerful
tool for this purpose [2]. Real-world datasets are often dynamic, with data
points being replaced, weights adjusted, entries removed, or new data points
added over time. Therefore, addressing the concept of dynamic optimal
transport becomes essential.

  * •

Time Series Analysis. Optimal Transport can serve as a metric in time series
analysis [7]. The main intuition lies in the smooth transition of states
between time points in a time series. The smoothness implies the potential to
iteratively refine a new solution based on the previous one, circumventing the
need for a complete recomputation.

  * •

Neuroimage analysis [14, 18]. In the medical imaging applications, we may want
to compute the change trend of a patient’s organ (e.g., the MRI images of
human brain over several months), and the differences are measured by the
optimal transportation cost. Since the changes are often local and small, we
can apply our method to quickly update the cost over the period.

  * •

Domain adaptation. It has been shown that OT can be used to compute the
transformation from one dataset (the source domain) to another dataset (the
target domain) [12]. But the datasets may need to be updated over time (e.g.,
in social network the user’s data can be changed), and so we can apply our
method to quickly update the transformation between different domains.

Existing methods, such as the Sinkhorn algorithm [9] and the Network Simplex
algorithm [23], are not adequately equipped to handle this dynamic problem.
Upon any modification to the cost matrix, the Sinkhorn algorithm requires at
least one regularization of the solution matrix, while the Network Simplex
algorithm needs to traverse all edges at least once to ensure optimality.
Consequently, the time complexity of these algorithms on the dynamic model is
$\Omega(\absolutevalue{V}^{2})$.

Our novel algorithm and data structure present an efficient alternative,
offering an $O(s\absolutevalue{V})$ solution for handling evolving datasets,
where $s$ is determined by the magnitude of the modification. In practice, $s$
is much less than $n$, making our algorithm an invaluable tool for a wide
range of applications in data science and beyond.

###  1.1 Related Works

Exact Solutions. Over the years, researchers have developed linear program
minimum cost flow algorithms to address discrete optimal transport problems.
The simplex method by Danzig et al. [10] offers a solution for general linear
programs. Despite its worst-case exponential time complexity, Spileman and
Teng [30] showed that its average time complexity is polynomial through smooth
analysis. Cunningham [8] adapted the simplex method for minimum cost flow
problems. Orlin [23] enhanced the network simplex algorithm with cost scaling
techniques. Recently, Lee and Sidford’s algorithm [20] presented an algorithm
with a time complexity upper bound of
$\tilde{O}(\absolutevalue{E}\sqrt{\absolutevalue{V}}\log^{2}U)$, and Li et al.
[21] proposed a parallel algorithm.

Approximations. Sherman [27] proposed a $(1+\epsilon)$ approximation
algorithm. Pele and Werman [24] introduced the FastEMD algorithm that applies
classic algorithms on a heuristic sketch of the input graph. Cuturi [9] used
Sinkhorn-Knopp iterations to approximate the optimal transport problem by
adding regular entropies to the optimization target. Recent optimizations on
the Sinkhorn algorithm include Overrelaxing Sinkhorn by Thibault et al. [33]
and the Screening Sinkhorn algorithm by Alaya et al. [1].

Dynamic Trees. Our data structure incorporates techniques from dynamic trees,
which maintain and manipulate tree-like structures under specific operations.
Sleater and Tarjan [28] proposed the Link/Cut tree data structure for path
aggregate maintenance. Tarjan and Vishkin [32] introduced the Euler Tour
Technique for maintaining a tree as a cyclic ordered set of its Euler Tour
sequence. Tarjan [31] later used splay trees to maintain such information,
leading to the development of Euler Tour Trees.

Search Trees. Our data structure utilizes high-dimensional extensions of skip
lists to maintain a 2-dimensional Euler Tour sequence. Existing high-
dimensional data structures based on self-balanced binary search trees, such
as Bentley’s $k$-d tree [3], are not suitable as they do not support cyclic
ordered set maintenance. Pugh’s concept of skip lists [26], which are linked
lists with additional layers of pointers for element skipping, is adapted in
our context to form skip orthogonal lists.

###  1.2 Overview of Our Algorithm

Our algorithm for the dynamic optimal transportation problem employs two key
strategies:

First, dynamic optimal transport operations are reduced to simplex iterations.
Our technique, grounded on the Simplex method, operates by eliminating the
smallest cycle in the graph. We assume that modifications influence only a
small portion of the result, requiring only a few simplex iterations. However,
existing algorithms like the Network Simplex Algorithm perform poorly under
dynamic modifications as they require scanning all edges at least once to
ensure global minimum attainment.

Second, a data structure is provided for performing each simplex iteration
within expected linear time relative to the supply and demand set. Our data
structure, described in Section 3, employs the Euler Tour Technique. We adapt
skip lists to maintain the cyclic ordered set produced by the Euler Tour
Technique and introduce an additional dimension to create a Skip Orthogonal
List. This structure aids in maintaining information about $A\times
B\mapsto\mathbb{R}$, crucial for our pivoting step.

Figure 1: 1D Euler Tour Tree with skip list

Figure 2: 2D Euler Tour Tree with Skip Our Orthogonal List

The rest of the paper is organized as follows:

  * •

Section 3 will give a introduction on the data structure Skip Orthogonal List
we innovated, where subsection 3.1 explains how the data structure is
organized and subsection 3.2 uses cut operation as an example to demonstrate
one of the updates on this data structure.

  * •

Section 4 explains how to use our data structure to solve the dynamic optimal
transport problem, by using subsection 4.1 to show that dynamic optimal
transport model could be reduced to simplex iterations, while using subsection
4.2 to show how our data structure could be used to improve the performance of
each simplex iteration.

##  2 Preliminaries

###  2.1 Optimal Transport

Let $A$ and $B$ represent source and target point sets, respectively, with
discrete probability distributions $\alpha:A\mapsto\mathbb{R}_{\geq 0}$ and
$\beta:B\mapsto\mathbb{R}_{\geq 0}$, such that $\sum_{a\in
A}\alpha_{a}=\sum_{b\in B}\beta_{b}=1$. The cost function $C:A\times
B\mapsto\mathbb{R}$ includes entries $c_{ab}$ denoting the cost of
transporting a unit of mass from point $a\in A$ to point $b\in B$. The
discrete optimal transport problem can be formulated as (1).

$\displaystyle W(\alpha,\beta,C)\triangleq$ $\displaystyle\min_{X:A\times
B\mapsto\mathbb{R}_{\geq 0}}\sum_{a\in A}\sum_{b\in B}c_{ab}x_{ab}$ (1)
subject to $\displaystyle\begin{cases}\sum_{b\in B}x_{ab}=\alpha_{a}&\forall
a\in A\\\ \sum_{a\in A}x_{ab}=\beta_{b}&\forall b\in B\end{cases}$

Since Problem (1) is a network flow problem, it can be transformed into
Problem (2) by adding infinity-cost edges, as constraint (3) indicates, and
redefining point weights as (4) [25]. Note that in (2), the input weight $w$
must always satisfy $\sum_{v\in V}w_{v}=0$. Otherwise, the constraints cannot
be satisfied.

$\displaystyle W(w,C)\triangleq$ $\displaystyle\min_{X:V\times
V\mapsto\mathbb{R}_{\geq 0}}\sum_{u\in V}\sum_{v\in B}c_{uv}x_{uv}$ (2)
subject to $\displaystyle\sum_{v\in V}x_{uv}-\sum_{v\in
V}x_{vu}=w_{u}\quad\forall u\in V$ $c_{aa}=c_{ba}=c_{bb}=+\infty\quad\forall
a\in A,b\in B$ (3) $w_{u}=\begin{cases}\alpha_{a}&\text{if }u=a\in A\\\
-\beta_{b}&\text{if }u=b\in B\end{cases}$ (4)

For a basic transportation plan $X$, we notice that the basis variables always
form a spanning tree of the complete directed graph with self loops
$G(V,V^{2})$ [8]. Let the dual variables be $\pi:V\mapsto\mathbb{R}$,
satisfying constraint (5). We may then define the adjusted cost function
$C^{\pi}:V^{2}\mapsto\mathbb{R}$ as $C^{\pi}_{uv}\triangleq
C_{uv}+\pi_{v}-\pi_{u}$, where $C^{\pi}_{uv}$ represents the simplex
multipliers for the linear program [23].

$x_{uv}\text{ is a basic variable }\implies\pi_{u}-\pi_{v}=C_{uv}$ (5)

###  2.2 Euler Tour Technique

The Euler Tour Technique is a method for representing a tree $T(V,E)$ as a
cyclic ordered set of length
$\absolutevalue{V}+2\absolutevalue{E}=O(\absolutevalue{V})$ [31]. Given a tree
$T(V,E)$, we construct a directed graph $D_{T}(V,E^{\prime})$ as follows:

  * •

For each vertex $v\in V$, add the self-loop $(v,v)$ to $E^{\prime}$.

  * •

For each undirected edge $\quantity{u,v}\in E$, add two directed edges $(u,v)$
and $(v,u)$ to $E^{\prime}$.

Since the difference of In-Degree and Out-Degree of each vertex in $D_{T}$ is
0, $D_{T}$ must contain an Euler Tour;

###### Definition 1

Define the Euler Tour representation of tree $T$ as an arbitrary sequence of
Euler Tour of $D_{T}$ represented by edges. Denote $E(D_{T})$ as the circular
ordered set with regard to the Euler Tour.

Through Definition 1, we can reduce edge linking, edge cutting, sub-tree
weight modification and sub-tree weight querying to constant number of element
insertion, element deletion, range weight modification and range weight
querying on a circular ordered set [31]. We show in Section 2.1 that the
dynamic optimal transport could be reduced to the 2 dimensional version of
these four operations, where $V$ would be interpreted as the union of supply
node set and demand node set (i.e. dual variables of the OT problem), and $E$
would be explained as the set of edges that are potentially non-zero, i.e., a
feasible basis for the primal variables. Euler Tour Technique is crucial to
our work.

###  2.3 Orthogonal Lists and Skip Lists

Skip Lists are probabilistic data structures that extend a singly linked list
with forward links at various levels, improving search, insertion, and
deletion operations. Figure 2 demonstrates a skip list. Each level contains a
circular linked list, where lists at a higher level is a subset of lists at a
lower level and the bottom level contains all elements. Nodes on the same
level is painted by the same color and linked horizontally. Corresponding
elements in adjacent lists are connected by vertical pointers. We apply this
skipping technique to circular singly linked lists in our work.

Most self-balanced binary search trees support lazy propagation techniques,
allowing range modifications within $O(\log n)$ time, where $n$ is the
sequence length maintained by the tree [29]. This technique is commonly used
in dynamic trees for network problems [31].

A $k$-dimensional Orthogonal List has $k$ orthogonal forward links and reduces
to a standard linked list when $k=1$. Orthogonal lists, which can be singly,
doubly, or circularly linked, can maintain information mapped from the
Cartesian product of ordered sets, such as sparse tensors [6].

Figure 2 demonstrates an orthogonal list maintaining a $3\times 4$ matrix.
Each node has 2 forward links, denoted by row links (yellow) and column links
(green). Row links connect elements in each row into a circular linked list
horizontally and column links connect elements in each column into a circular
linked list vertically.

Figure 3: An Orthogonal List Example

##  3 Skip Orthogonal List

This section introduces our novel data structure Skip Orthogonal List. Section
4 uses this data structure for fast update in the dynamical optimal transport
problem.

Formally, with the help of a Skip Orthogonal List, we can maintain a forest
$T(V,E)$ with at most two trees, and a function
$C^{\pi}:V^{2}\mapsto\mathbb{R}$, supporting the following operations:

  * •

Cut. Given $\quantity{u,v}$, remove edge $\quantity{u,v}$ from the forest,
splitting a tree into two disjoint trees, provided that the forest consists of
exactly one tree and $\quantity{u,v}\in E$. Let the connected component
containing $u$ form the vertex set $V_{1}$, and the connected component
containing $v$ form the vertex set $V_{2}$.

  * •

Insert. Add a new node to the forest, provided that the forest consists of
exactly one tree. Let the original nodes form the vertex set $V_{1}$, and the
new node itself form the vertex set $V_{2}$.

  * •

Range Add. Given $x$, for each $(u,v)\in V^{2}$, update $c^{\pi}_{uv}$ as
equation (6), provided the forest contains exactly two trees.

$\displaystyle c^{\pi}_{uv}\leftarrow c^{\pi}_{uv}+\begin{cases}0&(u,v)\in
V_{1}\times V_{1}\\\ -x&(u,v)\in V_{1}\times V_{2}\\\ x&(u,v)\in V_{2}\times
V_{1}\\\ 0&(u,v)\in V_{2}\times V_{2}\end{cases},$ (6)

  * •

Link. Given $(u,v)$, add edge $(u,v)$ to the forest, connecting two disjoint
trees into a single tree, provided that $u$ and $v$ are disconnected.

  * •

Global Minimum Query. Return the minimum value of $C^{\pi}$ on the tree,
provided that the forest is connected.

For the remaining of the section, we will construct a data structure in
expected $O(\absolutevalue{V}^{2})$ amount of space where each operation can
be done within expected $O(\absolutevalue{V})$ amount of time. Section 3.1
shows the overall structure of the data structure and how to query this data
structure. Section 3.2 will show the cut operation based on this structure.
For other operations (linking, insertion, range adding), we have put it in the
appendix for due to the limited space.

###  3.1 Overall Structure

As Figure 5 shows, a Skip Orthogonal List is a layered collection of
Orthogonal lists, just like the fact that a skip list is a hierarchical view
of a list. Each layer has fewer elements than the one below it, and the
elements are evenly spaced out. The bottom layer has all the elements while
the top layer has the least. Formally, it can be defined as Definition 2.

Figure 4: Detailed Structure of Skip Orthogonal List

Figure 5: Node Structure of Skip Orthogonal List

###### Definition 2

Given a parameter $p$ and a cyclic ordered set $S$, a 2 Dimensional Skip
Orthogonal List is an infinite collection of 2 Dimensional Circular Orthogonal
Lists $\mathcal{L}=\quantity(L_{0},L_{1},\cdots)$ , where

  * •

$\quantity{h_{s}:s\in S}$ are $\absolutevalue{S}$ independent random variable.
The distribution is a geometric distribution with parameter $p$

  * •

For each $i\in\mathbb{N}$, let $L_{i}$ maintain all elements in $S_{i}^{2}$,
where $S_{i}$ is the cyclic sequence formed by $\quantity{s\in S\mid h_{s}\geq
i}$

Assume we use this data structure to maintain several information about
$E(D_{T})\times E(D_{T})$. Since
$\absolutevalue{E(D_{T})}=O(\absolutevalue{V})$ as discussed in Section 3,
from Definition 2, the expected space complexity is $O(\absolutevalue{V}^{2})$
(due to the space limit, we put our proof in appendix).

Now we augment this data structure to store some additional information for
range adding and global minimum query. We first adapt the definition of
”dominate” to the 2 dimensional case defined as Definition 3 . In Figure 7 and
Figure 7, the blue nodes dominate itself and all yellow nodes, while the red
nodes dominate every node in the orthogonal list.

Figure 6: Domination in Skip List

Figure 7: Domination in 2D Skip List

###### Definition 3

In a Skip Orthogonal List $\mathcal{L}$, $x_{i_{1}i_{2}}$ dominates
$x_{j_{1}j_{2}}$ at level $l$ if and only if $(x_{i_{1}i_{2}}=x_{j_{1}j_{2}}$
or the following 3 conditions are all satisfied:

  * •

$h_{x_{i_{1}}}\geq l$ and $h_{x_{i_{2}}}\geq l$

  * •

$h_{x_{i_{1}+1}}<l$ and $h_{x_{i_{1}+2}}<l$ and … and $h_{x_{j_{1}-1}}<l$ and
$h_{x_{j_{1}}}<l$

  * •

$h_{x_{i_{2}+1}}<l$ and $h_{x_{i_{2}+2}}<l$ and … and $h_{x_{j_{2}-1}}<l$ and
$h_{x_{j_{2}}}<l$

Similar techniques can be used to define Skip Orthogonal Lists of higher
dimensions, but for dynamic optimal transport operations, since we now only
need to maintain the adjusted cost function $C^{\pi}$, we only need a 2
dimensional instance for the circular ordered set $E(D_{T})$. We use
$x_{\tilde{u}\tilde{v}}$ to maintain $c^{\pi}_{uv}$ where $\tilde{u}$ and
$\tilde{v}$ denote the self loop $(u,u)$ and $(v,v)$ respectively and fill the
remaining entries of the Skip Orthogonal Lists as the identity for the global
minimum query, i.e. $+\infty$.

We now augment the Skip Orthogonal List defined in 2. For each node
$x_{i_{1}i_{2}}$ at orthogonal list $L_{l}$, beside the 2 forward links and 2
backward links, we add the following attributes as shown in Figure 5:

  * •

min to maintain the minimum amongst all nodes dominated by it if $l>0$

  * •

tag to maintain the lazy tag amongst for all nodes dominated by it if $l>0$

  * •

son to point to $x_{i_{1}i_{2}}$ at orthogonal list $L_{l-1}$ if $l>0$

  * •

parent to point to the node that dominates it if $L_{l+1}$ is not empty

###  3.2 Update (Cut)

In this subsection, we focus on the Update (Cut) algorithm for a 2D Skip
Orthogonal List as an example. Figure 9 and Figure 9 show the basic idea of
the cutting process.

Figure 8: Vertical View of Update (Cut)

Figure 9: 3D View of Update (Cut)

The algorithm can be crudely described as follows:

  1. 1.

Find the 2 rows and 2 columns representing directed edge $(u,v)$ and $(v,u)$
in $E(D_{T})$, i.e. the transparent nodes in Figure 9 or the gray nodes in
Figure 9.

  2. 2.

Push down the tag attribute of all nodes alongside the rows and columns, i.e.
the red nodes and transparent nodes in Figure 9. Only the nodes of which the
set of nodes dominated by them may change during step 3 need to push down.

  3. 3.

Cut the rows and columns, warping up the forward links and backward links of
points alongside, as Figure 9 demonstrates. This will cut the original Skip
Orthogonal List into 4 small lists.

  4. 4.

Update the min attribute of all nodes whose descendants were affected in step
3, i.e. the red nodes in Figure 9.

  5. 5.

Return the 4 small lists cutted out in step 3.

From the discussion above, we can see that this algorithm is an augmentation
of lazy propagation technique which is widely used in normal balanced tree,
and thus solves our problem. We can see that the time complexity is linear
corresponding to the number of nodes dominated by the red nodes and
transparent nodes in figure 9. Similar to the proof of space complexity, the
expected time complexity is $O(\absolutevalue{V})$ with high probability.

##  4 Our Dynamic Network Simplex Method

The simplex method performs simplex iterations on some initial feasible basis
until optimality is obtained [10]. Simplex iterations can refine the current
solution. In each simplex iteration, some variable with negative simplex
multiplier is selected by some strategy to perform a pivot subroutine. One
possible strategy is to pivot in the variable with the smallest simplex
multiplier. Section 4.1 focuses on defining dynamic optimal transport
operations and using simplex iterations to solve it, while Section 4.2 focuses
on analyzing the details in each simplex iteration to optimize it.

###  4.1 Dynamic Optimal Transport Operations

In an Optimal Transport problem, suppose the supply set $A$ and the demand set
$B$ are mapped to some points in some metric space $(\mathcal{X};d)$ by
function $f$, e.g. $X$ could be the Euclidean Space $\mathbb{R}^{k}$ while $d$
could be the Euclidean Distance $\ell_{2}$, while the cost function $C_{ab}$
is usually defined as $d(f_{a},f_{b})$. Let $w:V\mapsto\mathbb{R}$ denote the
weight function as equation (4) defined where $V=A\cup B$. A Dynamic Optimal
Transport algorithm supports the following 4 types of update and 1 type of
online query:

  * •

Space Position Modification. Select some supply or demand point $v\in V$ and
move $f_{v}$ to another point $x\in\mathcal{X}$. This usually results in
modification of an entire row or column in the cost matrix
$C\in\mathbb{R}^{A\times B}$.

  * •

Weight Modification. Select a pair of supply or demand point $u,v\in V$ with
some positive number $\Delta\in\mathbb{R}_{+}$. Let $w_{u}\leftarrow
w_{u}-\Delta$ while $w_{v}\leftarrow w_{v}+\Delta$.

  * •

Point Deletion. Select $v\in V$ with $w_{v}=0$. Remove $v$ from $V$.

  * •

Point Insertion. Select $v\notin V$ and $f_{v}\in\mathcal{X}$. Let $w_{v}=0$
and insert $v$ into set $V$

  * •

Query. Answer the current Optimal Transport plan and value.

These updates do not modify the sum of all points in supply and demand set
$V$, so $\sum_{v\in V}w_{v}\equiv 0$ and thus a transport plan always exists.
Therefore we can reduce these updates to updates on simplex basis:

  * •

Space Position Modification. The original optimal solution is primal feasible
but not primal optimal, i.e. not dual feasible. We perform primal simplex
method based on the original optimal solution. When moving point $v$, we first
update cost matrix $C$, dual variable $\pi$ and modified cost $C^{\pi}$ to
meet constraint (5). After that, we perform simplex iterations while the
minimum value of the adjusted cost $C^{\pi}$ is negative [25].

  * •

Weight Modification. The original optimal solution is dual feasible but not
primal feasible. We perform dual simplex iterations based on the original
optimal solution. Assume we try to decrease $w_{u}$ and increase $w_{v}$ by
$\Delta$, we send $\Delta$ amount of flow from $u$ to $v$ on the residual
network similar to the shortest path augmenting method [11, 25]: We send flow
on basic variables $X_{B}$. If some variable needs to be pivoted out, we pivot
in the variable with the smallest adjusted cost.

  * •

Point Deletion & Point Insertion. As the deleted/inserted point weighs $0$,
whether inserting or deleting the point does not influence our result. In
order to influence the Optimal Transport plan, it must be combined with
several Weight Modification operations for weight assignment. Therefore, we
can maintain a node pool to maintain all supply and demand nodes with 0
weight. Each Point Insertion operation takes some point in this pool and move
it to the correct place while each Point Deletion operation returns a node to
the pool.

Our solution updates the optimal transport plan as soon as an update takes
place, so we can answer the online query about optimal transport plan and
value. Since we are only modifying one or two nodes in the dataset,
intuitively the optimal transport plan should not change much, and we do not
need many simplex iterations to obtain optimality. Assume we need $s$ simplex
iterations, where we assume $s\ll\absolutevalue{V}$. Therefore the time
complexity for our algorithm is $O(s{\rm Time}_{\rm iteration})$ where ${\rm
Time_{iteration}}$ is the time of each simplex iteration.

###  4.2 Simplex Iteration Details

As Section 4.1 describes, dynamic operations on Optimal Transport could be
reduced to simplex iterations. In this section, we review on the operations
used in conventional network simplex algorithm and use data structure designed
in Section 3 to maintain $C^{\pi}$.

Conventional network simplex method is simplex method simplified by some graph
properties [8]. The steps of a (network) simplex iteration can be described
as:

  1. 1.

Select Variable to Pivot in. Select variables with the smallest adjusted cost
$C^{\pi}$ to pivot in. Assume $x_{i_{\rm in}j_{\rm in}}$ is to be pivoted in.

  2. 2.

Update Primal Solution. Send circular flow in the cycle formed by adding the
new variable to the current basis to perform the pivot operation until some
basic variable is pivoted out, i.e. some basic variable in the reverse
direction runs out of flow. Assume $x_{i_{\rm out}j_{\rm out}}$ is to be
pivoted out

  3. 3.

Update Dual Solution. Update dual variables $\pi$ and modified cost $C^{\pi}$
to meet constraint (5), as the new basis as $C^{\pi}$ will soon be queried in
the next simplex iteration.

The selecting step performs a query on the data structure on $C^{\pi}$ about
the minimum element and the dual updating performs an update. Though the
primal updating step can be done within $O(\log\absolutevalue{V})$ amount of
time [31], conventional network simplex maintain $C^{\pi}$ through brute
force. Upon an query, conventional network simplex brutally traverse through
all adjusted costs and selects the minimum, while upon an update, conventional
network simplex updates the adjusted cost one by one after the adjusted cost
is computed. This indicates that the time complexity of each simplex iteration
is $O(\absolutevalue{V}^{2})$.

Our goal is to maintain $C^{\pi}$ so that it can answer the global minimum
query and perform update when the primal basis changes. In a simplex method,
when we decide to pivot in the variable $x_{i_{\rm in}j_{\rm in}}$, we update
the dual variables as equation (7) describes where $V_{1}$ is the set of nodes
connected to $i_{\rm in}$ and $V_{2}$ is the set of nodes connected to $j_{\rm
in}$ after the edge $x_{i_{\rm out}j_{\rm out}}$ is cut [8]. $\pi$ are the
dual variables before pivoting and $\pi^{\prime}$ are the dual variables after
pivoting.

$\pi_{i}^{\prime}\leftarrow\pi_{i}+\begin{cases}C_{i_{\rm in}j_{\rm
in}}^{\pi}&i\in V_{1}\\\ 0&i\in V_{2}\end{cases}$ (7)

From the definition of the adjusted cost function $C^{\pi}$, we obtain that
the update objective can be formulated as (8) where $C^{\pi}$ is the adjusted
cost function with regard to the old dual variables $\pi$ while
$C^{\pi^{\prime}}$ regards the new dual variables $\pi^{\prime}$. Therefore,
Algorithm 1 can help us with the adjusted cost function update. Our Skip
Orthogonal List in Section 3 is capable of performing cut, add and link
operation within $O(\absolutevalue{V})$ amount of time. Therefore Theorem 4.1
holds.

$C^{\pi^{\prime}}_{ij}=C_{ij}^{\pi}+\pi^{\prime}_{i}-\pi^{\prime}_{j}=C^{\pi}_{ij}+\begin{cases}0&(i,j)\in
V_{1}\times V_{1}\\\ -C^{\pi}_{i_{\rm in}j_{\rm in}}&(i,j)\in V_{1}\times
V_{2}\\\ C^{\pi}_{i_{\rm in}j_{\rm in}}&(i,j)\in V_{2}\times V_{1}\\\
0&(i,j)\in V_{2}\times V_{2}\end{cases}$ (8)

###### Theorem 4.1

Each simplex iteration in conventional network simplex can be completed within
expected $O(\absolutevalue{V})$ amount of time.

Input: Adjusted cost function $C^{\pi}$, entering variable $x_{i_{\rm
in}j_{\rm in}}$, leaving variable $x_{i_{\rm out}j_{\rm out}}$

Output: Updated adjusted cost function $C^{\pi^{\prime}}$

1 $t\leftarrow C^{\pi}_{a_{\rm in}b_{\rm in}}$

2
$C^{\pi}_{11},C^{\pi}_{12},C^{\pi}_{21},C^{\pi}_{22}\leftarrow\textit{cut}(C^{\pi},a_{\rm
in},b_{\rm out})$

3 $\textit{add}(C^{\pi}_{12},-t)$

4 $\textit{add}(C^{\pi}_{21},t)$

5
$C^{\pi^{\prime}}\leftarrow\textit{link}(C^{\pi}_{11},C^{\pi}_{12},C^{\pi}_{21},C^{\pi}_{22},a_{\rm
in},b_{\rm in})$

Algorithm 1 Adjusted Cost Matrix Update

##  5 Experiments

All the results were obtained on a server equipped with AMD Ryzen Threadripper
PRO 5975WX 32-Cores CPU and 512GB main memory with frequency 3200 MHz; the
data structures are implemented in C++20 compiled by G++ 13.1.0 on Ubuntu
22.04.3. The data structures are compiled to shared objects to be called by
Python 3.11.5. Our code uses Network Simplex library from [5] to get an
initial feasible flow and NumPy [16] to represent points in the space.

In our experiment, we used Network Simplex algorithm [23] and Sinkhorn
algorithm [9] from the Python Optimal Transport (POT) library [13]. We run our
algorithm in both Space Position Modification and Point Insertion scenario as
described in Section 4.1. We use the running time to measure the performance
of algorithms.

Datasets. We study the performance of our algorithm on both synthetic datasets
and real datasets. For synthetic datasets, we construct a mixture of 2
Gaussian distribution in $\mathbb{R}^{2}$, where points in the same Gaussian
distribution share the same label [4]. We also used MNIST dataset as a popular
real world dataset [19]. We partition the labels into 2 groups and run optimal
transport between these groups, as we wish to run optimal transport on large
datasets, which indicates that a single label cannot meet the number of nodes
we demand.

Setup. We set the Sinkhorn regularization parameter $\lambda$ as $0.1$ and
scale the median of the cost matrix $C$ to be 1. We test on sampled datasets
where we vary the the node size $\absolutevalue{V}$ up to $4\times 10^{4}$.
For each dataset, we test the static running time of POT on our machine and
executed each dynamic operations 100 times to calculate the means and standard
errors of our algorithm. For Space Position Modification, we randomly choose
some point and add a random Gaussian noise with a variance of 0.5 in each
dimension to it. For Point Insertion, we randomly select a point in the
dataset that is not currently in the OT system to insert. We perform a query
after these updates and compare with static algorithms implemented in POT to
ensure our correctness.

Result and Analysis. We illustrate our experiment results in Figure 10 with
the help of Matplotlib [17]. In dynamic models, our algorithm is about 1000
times faster than Static Network Simplex algorithm and 10 times faster than
Sinkhorn Algorithm when $\absolutevalue{V}$ reaches 40000 depending on the
amplitude of modification. As the number of nodes grow larger, our algorithm
will be even faster. This indicates that our algorithm is fast in the case
where simplex iterations are not much, as discussed in Section 4.1. Also, the
running time of our algorithm is in a super linear trend. As each simplex
iteration in our method is strictly linear in expectation theoretically, this
could be caused by increment in number of simplex iterations or decrement in
cache hit rate as the node size grows larger.

Figure 10: Experiment Result

##  6 Conclusion and Future Work

Our work provides a data structure which is useful in traditional network
simplex method. With the help of our data structure, the time complexity of
the whole pivoting process is $O(\absolutevalue{V})$ in expectation as Theorem
4.1 describes. Based on the fact that conventional simplex method performs
well in average cases, our algorithm has a better theoretical performance than
most existing algorithms.

However, our algorithm lead to several performance issues in practice. First,
as our algorithm stores the entire 2D Skip Orthogonal List data structure, our
algorithm require additional space. Second, as our algorithm is based on
linked data structures, the cache hit rate is not high. Therefore, improving
our implementation and developing algorithms and data structures with similar
complexity but more memory friendly is our next goal.

## References

  * [1] Alaya, M.Z., Berar, M., Gasso, G., Rakotomamonjy, A.: Screening sinkhorn algorithm for regularized optimal transport. Advances in Neural Information Processing Systems 32 (2019)
  * [2] Alvarez-Melis, D., Fusi, N.: Geometric dataset distances via optimal transport. In: NeurIPS 2020. ACM (February 2020), https://www.microsoft.com/en-us/research/publication/geometric-dataset-distances-via-optimal-transport/
  * [3] Bentley, J.L.: Multidimensional binary search trees used for associative searching. Communications of the ACM 18(9), 509–517 (1975)
  * [4] Blum, A., Hopcroft, J., Kannan, R.: Foundations of data science. Cambridge University Press (2020)
  * [5] Bonneel, N., van de Panne, M., Paris, S., Heidrich, W.: Displacement Interpolation Using Lagrangian Mass Transport. ACM Transactions on Graphics (SIGGRAPH ASIA 2011) 30(6) (2011)
  * [6] Butterfield, A., Ngondi, G.E., Kerr, A.: A dictionary of computer science. Oxford University Press (2016)
  * [7] Cheng, K., Aeron, S., Hughes, M.C., Miller, E.L.: Dynamical wasserstein barycenters for time-series modeling. Advances in Neural Information Processing Systems 34, 27991–28003 (2021)
  * [8] Cunningham, W.H.: A network simplex method. Mathematical Programming 11, 105–116 (1976)
  * [9] Cuturi, M.: Sinkhorn distances: Lightspeed computation of optimal transport. Advances in neural information processing systems 26 (2013)
  * [10] Dantzig, G.B., Orden, A., Wolfe, P., et al.: The generalized simplex method for minimizing a linear form under linear inequality restraints. Pacific Journal of Mathematics 5(2), 183–195 (1955)
  * [11] Edmonds, J., Karp, R.M.: Theoretical improvements in algorithmic efficiency for network flow problems. Journal of the ACM (JACM) 19(2), 248–264 (1972)
  * [12] Flamary, R., Courty, N., Tuia, D., Rakotomamonjy, A.: Optimal transport for domain adaptation. IEEE Trans. Pattern Anal. Mach. Intell 1, 1–40 (2016)
  * [13] Flamary, R., Courty, N., Gramfort, A., Alaya, M.Z., Boisbunon, A., Chambon, S., Chapel, L., Corenflos, A., Fatras, K., Fournier, N., et al.: Pot: Python optimal transport. The Journal of Machine Learning Research 22(1), 3571–3578 (2021)
  * [14] Gramfort, A., Peyré, G., Cuturi, M.: Fast optimal transport averaging of neuroimaging data. In: Information Processing in Medical Imaging: 24th International Conference, IPMI 2015, Sabhal Mor Ostaig, Isle of Skye, UK, June 28-July 3, 2015, Proceedings 24. pp. 261–272. Springer (2015)
  * [15] Haker, S., Zhu, L., Tannenbaum, A., Angenent, S.: Optimal mass transport for registration and warping. International Journal of computer vision 60, 225–240 (2004)
  * [16] Harris, C.R., Millman, K.J., van der Walt, S.J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N.J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M.H., Brett, M., Haldane, A., del Río, J.F., Wiebe, M., Peterson, P., Gérard-Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., Oliphant, T.E.: Array programming with NumPy. Nature 585(7825), 357–362 (Sep 2020). https://doi.org/10.1038/s41586-020-2649-2, https://doi.org/10.1038/s41586-020-2649-2
  * [17] Hunter, J.D.: Matplotlib: A 2d graphics environment. Computing in Science & Engineering 9(3), 90–95 (2007). https://doi.org/10.1109/MCSE.2007.55
  * [18] Janati, H., Bazeille, T., Thirion, B., Cuturi, M., Gramfort, A.: Group level meg/eeg source imaging via optimal transport: minimum wasserstein estimates. In: Information Processing in Medical Imaging: 26th International Conference, IPMI 2019, Hong Kong, China, June 2–7, 2019, Proceedings 26. pp. 743–754. Springer (2019)
  * [19] LeCun, Y., Cortes, C., Burges, C.: Mnist handwritten digit database. ATT Labs [Online]. Available: http://yann.lecun.com/exdb/mnist 2 (2010)
  * [20] Lee, Y.T., Sidford, A.: Path finding methods for linear programming: Solving linear programs in o (vrank) iterations and faster algorithms for maximum flow. In: 2014 IEEE 55th Annual Symposium on Foundations of Computer Science. pp. 424–433. IEEE (2014)
  * [21] Li, W., Ryu, E.K., Osher, S., Yin, W., Gangbo, W.: A parallel method for earth mover’s distance. Journal of Scientific Computing 75(1), 182–197 (2018)
  * [22] Métivier, L., Brossier, R., Mérigot, Q., Oudet, E., Virieux, J.: Measuring the misfit between seismograms using an optimal transport distance: Application to full waveform inversion. Geophysical Supplements to the Monthly Notices of the Royal Astronomical Society 205(1), 345–377 (2016)
  * [23] Orlin, J.B.: A polynomial time primal network simplex algorithm for minimum cost flows. Mathematical Programming 78, 109–129 (1997)
  * [24] Pele, O., Werman, M.: Fast and robust earth mover’s distances. In: 2009 IEEE 12th international conference on computer vision. pp. 460–467. IEEE (2009)
  * [25] Peyré, G., Cuturi, M., et al.: Computational optimal transport: With applications to data science. Foundations and Trends® in Machine Learning 11(5-6), 355–607 (2019)
  * [26] Pugh, W.: Skip lists: a probabilistic alternative to balanced trees. Communications of the ACM 33(6), 668–676 (1990)
  * [27] Sherman, J.: Generalized preconditioning and undirected minimum-cost flow. In: Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms. pp. 772–780. SIAM (2017)
  * [28] Sleator, D.D., Tarjan, R.E.: A data structure for dynamic trees. In: Proceedings of the thirteenth annual ACM symposium on Theory of computing. pp. 114–122 (1981)
  * [29] Sleator, D.D., Tarjan, R.E.: Self-adjusting binary search trees. Journal of the ACM (JACM) 32(3), 652–686 (1985)
  * [30] Spielman, D.A., Teng, S.H.: Smoothed analysis of algorithms: Why the simplex algorithm usually takes polynomial time. Journal of the ACM (JACM) 51(3), 385–463 (2004)
  * [31] Tarjan, R.E.: Dynamic trees as search trees via euler tours, applied to the network simplex algorithm. Mathematical Programming 78(2), 169–177 (1997)
  * [32] Tarjan, R.E., Vishkin, U.: Finding biconnected componemts and computing tree functions in logarithmic parallel time. In: 25th Annual Symposium onFoundations of Computer Science, 1984. pp. 12–20. IEEE (1984)
  * [33] Thibault, A., Chizat, L., Dossal, C., Papadakis, N.: Overrelaxed sinkhorn–knopp algorithm for regularized optimal transport. Algorithms 14(5), 143 (2021)
  * [34] Torres, L.C., Pereira, L.M., Amini, M.H.: A survey on optimal transport for machine learning: Theory and applications. arXiv preprint arXiv:2106.01963 (2021)