# Multimodal Shape Completion via Implicit Maximum Likelihood Estimation

First Author
Institution1
Institution1 address
<EMAIL_ADDRESS>Second Author
Institution2
First line of institution2 address
<EMAIL_ADDRESS>
###### Abstract

Shape completion is the problem of completing partial input shapes such as
partial scans. This problem finds important applications in computer vision
and robotics due to issues such as occlusion or sparsity in real-world data.
However, most of the existing research related to shape completion has been
focused on completing shapes by learning a one-to-one mapping which limits the
_diversity_ and _creativity_ of the produced results. We propose a novel
multimodal shape completion technique that is effectively able to learn a one-
to-many mapping and generates diverse complete shapes. Our approach is based
on the conditional Implicit Maximum Likelihood Estimation (IMLE) technique
wherein we condition our inputs on partial 3D point clouds. We extensively
evaluate our approach by comparing it to various baselines both quantitatively
and qualitatively. We show that our method is superior to alternatives in
terms of completeness and diversity of shapes.

##  1 Introduction

As 3D scanning devices are now widely available, scenes and objects in the
form of point clouds are rapidly becoming more and more prevalent. However,
many such objects are not fully scanned due to occlusion; as a result, the
scanned shapes may be missing substantial portions of the complete shape.
Therefore, the ability to complete such shapes by generating the missing
portions is desired. More precisely, the objective of shape completion is to
meaningfully complete a partial shape and obtain a complete shape without
violating the geometry of the partial shape.

Early work has focused on generating a single complete shape; the mapping from
the input and output is one-to-one. However, shape completion is ill-posed by
nature as the target shape is not unique – given a partial shape, there are
multiple ways to generate the missing portions, especially when the input
point cloud is noisy and coarse, which is the case in real-world scanned data.
Therefore, given a partial shape, it is desired to generate a _set_ of many
possible complete shapes, all of which respecting the given input partial
shape. Not only is multimodality a way to resolve ill-posedness, but it also
adds more _diversity_ and _creativity_ to the final outputs.

Figure 1: Our network receives a partial point cloud $\mathcal{X}_{i}$ and a
complete point cloud (in training time) $\mathcal{Y}_{i}$. The latent code
generated from the incomplete shape, $\mathbf{x}_{i}$ along with a noise,
$\mathbf{z}_{j}$ is given to generator $G$ to be mapped to
$\mathbf{\tilde{y}}_{i,j}$. Closest $\mathbf{\tilde{y}}_{i,j}$ to the latent
code $\mathbf{y}_{i}$ of the complete point cloud is found and passed to
decoder $D$ to produce a novel and complete shape $\mathcal{\tilde{X}}_{i,j}$.

Recently, Wu et al. [wu2020multimodal] proposed the task of multimodal shape
completion, where the goal is to generate various complete shapes for a single
partial shape. This is challenging, because the training data only contains
one ground truth complete shape for each partial shape. In [wu2020multimodal],
a conditional GAN-based model was proposed to learn a one-to-many mapping from
partial shapes to complete shapes.

While [wu2020multimodal] took a step towards multimodal shape completion, the
different shapes it generates for a given partial point cloud tend to be
similar in structure. For example, for an incomplete chair, the different
complete planes it generates are of a similar style and vary mostly in width
and height, as shown in Figure 4. In this paper, our goal is to generate
_structurally diverse_ complete shapes for the same incomplete shape. In
particular, consider the case where a part of the object is entirely missing,
the different generated complete shapes should ideally vary in the style of
the missing part while remaining compatible with the other parts that are
present.

Our key insight is that structural similarity between the generated complete
shapes comes from the tendency of GANs to collapse modes. To achieve
structural diversity, we must therefore get around mode collapse. To this end,
we propose a novel method by devising a model formulation that allows for
training with the method of Implicit Maximum Likelihood Estimation (IMLE).
Prior to our work, IMLE has never been successfully used for point cloud
generation. In this paper, we show how to do so, which results in a simple and
elegant model with fewer parameters than [wu2020multimodal]. Yet, despite this
simplicity, we demonstrate that our method is capable of generating
structurally diverse results. Quantitatively, our method achieves superior
performance than prior methods in terms of Total Mutual Difference, and
comparable or superior performance than prior methods in terms of
Unidirectional Haussdorf Distance (UHD). Moreover, we demonstrate the
capability of our method to generalize to new kinds of incomplete shapes.
Despite not having part-level supervision, our method can produce complete
shapes for incomplete shapes with an entire part missing.

##  2 Related Work

Since the missing regions of an incomplete point cloud are unknown and might
be geometrically complex, simple surface repair algorithms are not fully
effective in reproducing parts and geometric features. In fact, missing
regions should be better _generated_ after observing many relevant shapes to
learn their geometric properties. Therefore, here, we mostly provide a review
of deep generative models that are reasonable solutions for shape completion.
In addition, since our input is a point cloud, we provide a review on deep
learning models designed for point clouds, and we particularly discuss point
cloud completion methods. Point cloud completion may have similarities with
other problems, including point upsampling, 3D part assembly, 3D
reconstruction, and single view reconstruction. However, since the focus of
our paper is multimodal point cloud completion, we refrain from discussing
works related to those areas and refer the interested reader to the available
survey papers such as [egstar2020_struct, berger2014state].

###  2.1 Generative Models

Generating new images or 3D models is an important area in computer vision and
computer graphics. One of the most well-known techniques to generate new
shapes is _variational autoencoder_ [kingma2013auto] or VAE that tries to fit
a known distribution (e.g., Gaussian) to the latent space of an autoencoder so
that by sampling the known distribution, the decoder can produce novel data.
VAEs have been used in many applications, including generating new 3D datasets
and shape completion [stutz2020learning, egstar2020_struct].

_Generative Adversarial Networks_ (GANs) [goodfellow2014generative] are
alternative generative models composed of a generator and a discriminator,
where the generator tries to produce novel data that can fool the
discriminator, and the discriminator tries to distinguish the produced results
(i.e., fake) from real samples. There are various forms of GAN that have been
designed to tackle different problems, including conditional GANs that try to
generate data based on a condition such as an incomplete image for the purpose
of inpainting [yu2018generative, isola2017image] or a partial shape for shape
completion [chen2019unpaired]. In fact, GAN-based networks have been
successful in many 3D tasks and applications [egstar2020_struct].

Despite the popularity and success of GANs, they have limitations including
mode collapse. This problem is more evident when variety in the generated data
is required (such as the case of multimodal shape completion). This has been
the primary motivation behind IMLE, which is to avoid mode collapse and it has
been applied to various multimodal image generation problems including image
super-resolution [Li2020MultimodalIS], image synthesis from semantic layout
[Li2019DiverseIS] and image decompression [Peng2020GeneratingUA].

###  2.2 Deep Learning on Point Clouds

One of the possibilities is to feed a neural network by an array of points
(i.e., point clouds) capturing the geometry of shapes. One of the earliest
networks to deal with point clouds is PointNet [qi2017pointnet] that tries to
learn feature points to perform different tasks including segmentation or
classification. As PointNet only learns global features, PointNet++
[qi2017pointnet++] is introduced to learn multi-scale features and better
abstract local patterns. Due to the effectiveness of PointNet and PointNet++,
they have been used as the main component of many networks with various
applications including point to point translation [yin-sig18, yin-siga19],
edge detection on point clouds [wang2020], learning implicit surfaces
[erler2020], denoising [rakotosaona2020] and many more [guo2020].

###  2.3 Shape Completion

Classical methods. Shape completion has been studied in the area of mesh
processing even before the popularity of deep learning. For example, the
missing region of a mesh can be repaired by minimizing an energy function via
Laplacian or Poisson equation [sorkine2004least, zhao2007robust,
nealen2006laplacian, kazhdan2013screened, kazhdan2006poisson], benefiting from
shapes’ skeletal structures [Dpoints15, yin2014morfit], or retrieving a
template from a database and deforming it into the incomplete input point
cloud [pauly2005example, rock2015completing]. Interested readers can refer to
[berger2014state] for a thorough study on shape reconstruction and completion
for point clouds.

Deep shape completion. The missing region of a 3D shape, however, might be
geometrically complex that can be only _generated_ through a deep network
trained on a dataset with variety. As a result, deep learning techniques have
been proposed to leverage available 3D shape datasets to learn how to
reproduce the missing regions successfully.

In [dai2017shape], an autoencoder structure is used to first predict a
complete distance field of the voxelized incomplete input shape. The distance
field is then contrasted against a shape database to obtain the final shape
through a multi-resolution patch-based 3D shape synthesis. Stutz and Geiger
[stutz2020learning] provide a weakly supervised VAE-based method for shape
completion. They first train a denoising variational autoencoder on ShapeNet
in a supervised manner. They then fix and reuse the decoder to train a new
encoder on incomplete shapes where alignment of the predicted result with the
incomplete input shape is respected by applying a maximum likelihood (ML)
loss. ASFM-Net [xia2021asfm] maps the partial and complete input point clouds
into a common latent space to capture detailed shape priors and better respect
the input.

A series of methods have utilized a multi-stage approach to complete a shape.
Inspired from patch-match technique in image processing, Han et al.
[han2017high] employed two networks to produce missing parts of an incomplete
shape, where one network tries to respect the overall structure of the shape
and the other one produces local patches to repair the given incomplete shape.
PCN [yuan2018pcn] is a multi-stage encoder-decoder network where the decoder
first produces a coarse complete shape, and then for each point of the coarse
shape, a local patch is generated and deformed to reproduce details. Liu et
al. [liu2020morphing] also introduced a multi-stage method where a coarse
shape is first approximated by a set of parametric surfaces and then it is
refined by learning point-wise residuals. Similarly, Wang et al.
[wang2020cascaded] designed a network that performs the shape completion task
in two stages (i.e., coarse and dense reconstruction). They have also designed
a patch discriminator to guarantee that local patches resemble the same
pattern as those in the ground truth. In addition to the coarse/fine point
completion techniques, TopNet [tchapmi2019topnet] has been designed as a tree
structure where each branch represents a subset of the point cloud.

Detailed shape completion. There are recent works that have focused on
completing partial input shapes that attain a great amount of detail. These
works try to replicate existing features in the partial shape in the proper
location of the generated complete shape. Wen et al. [wen2020point] have
defined a skip-attention mechanism to avoid using a global representation that
may suffer from altering the initial geometry of the input shape. PF-Net
[huang2020pf] employs a multi-resolution encoding for the shape and tries to
only generate the missing parts of the partial shape to respect the input’s
features. Most recently, VRCNet [pan2021variational] benefits from the multi-
stage completion (i.e., coarse and fine stages) and is capable of replicating
fine details when they are available. Although these works produce nice
results when the inputs are detailed, our work is naturally different from
these methods. Here, we focus on completing input shapes that are relatively
coarse and therefore, multiple interpretations are possible. Only under this
setting, the multimodality is meaningful since when the partial input shape is
dense and detailed, the complete shape will be limited to only one ground
truth shape with the same features as the input. In Section 5 and
supplementary material, we show that the uni-modal version of our network
performs comparably to VRCNet for partial shapes for which many complete
shapes can be imagined.

GANs and multimodality. GAN-based techniques have also been employed for shape
completion. Wang et al. [wang2017shape] use a GAN in the form of a
convolutional autoencoder to produce a coarse shape which is then sliced into
images and passed to a recurrent convolutional network for refinement. To
setup an unpaired shape completion, two autoencoders are initially trained on
complete and incomplete shapes in [chen2019unpaired]. The latent code of an
incomplete shape is then mapped to a complete shape’s latent code via a
generator, and two discriminators are used to evaluate complete shapes and
their latent codes. The other generative method is proposed by Wu et al.
[wu2020multimodal] where they combine VAE and conditional GAN to achieve
multimodality. Later in this paper, we show that although [wu2020multimodal]
is successful in point cloud completion due to mode collapse, our IMLE-based
approach is more successful in producing shapes with more variety. A series of
experiments have been performed to examine the effectiveness of our method,
which is discussed in Section 5. In addition to GAN-based approaches that are
capable of providing multimodality, a voxel-based generative model (GCA) is
introduced by Zhang et al. [zhang2021learning] that produce results by
gradually completing voxels. Although we have provided a qualitative
comparison with GCA in the supplementary material, it might not be directly
relevant to our work as GCA is designed for voxels instead of point clouds,
and also it does not use GAN in its structure while we try to provide an
alternative to GAN.

##  3 Overview

In shape completion problem, given a partial point cloud, we want to predict a
complete point cloud based on the partially observed portion. Formally, given
an input $\mathbf{x}$ in an one-to-many problem setting, the goal is to learn
$p(\mathbf{y}|\mathbf{x})$ where $\mathbf{y}$ is a prediction. Here,
$\mathbf{x}$ is the partial point cloud and $\mathbf{y}$ is a predicted
complete point cloud. Since there are many possible complete shapes that are
consistent with $\mathbf{x}$, we want to generate non-deterministic
predictions $p(\mathbf{y}|\mathbf{x})$ instead of a single prediction. The
simplest way to generate non-deterministic predictions is to add a random code
input $\mathbf{z}\sim\mathcal{N}(0,\mathbf{I})$ to a deterministic function
$T_{\theta}$ where $\mathbf{y}\coloneqq T_{\theta}(\mathbf{x},\mathbf{z})$.
This is also known as an implicit generative model [mohamed2016learning]. One
way to train this implicit generative model is to use conditional GAN, in this
context $T_{\theta}(\mathbf{x},\mathbf{z})$ is often known as the generator.
However, the generator may fail to generate some plausible outputs due to mode
collapse. Therefore, we use Implicit Maximum Likelihood Estimation (IMLE) to
train the generative model which overcomes the issue of mode collapse.

Unlike GANs, IMLE does not suffer from mode collapse as it reverses the
direction of how generated samples are matched to the target data: instead of
making the generated samples similar to some target data, IMLE ensures that
every target data has a similar sample being generated. In addition, IMLE does
not require an adversarial objective which makes it easy to train. If we
denote the generator as $T_{\theta}(·,·)$, an input point cloud as
$\mathbf{x}_{i}$ and the corresponding ground truth as $\mathbf{y}_{i}$, a
latent code $\mathbf{z}_{i,j}$, then IMLE would optimize the following
equation:

$\min_{\theta}\mathbb{E}_{\mathbf{z}_{1,1},\ldots,\mathbf{z}_{n,m}\sim\mathcal{N}(0,\mathbf{I})}\left[\sum_{i=1}^{n}\min_{j\in\\{1,\ldots,m\\}}d(T_{\theta}(\mathbf{x}_{i},\mathbf{z}_{i,j}),\mathbf{y}_{i})\right]$
(1)

where $d(·,·)$ is a distance metric and $m$ is a hyperparameter. Algorithm 4
shows the conditional IMLE training procedure.

Conditional IMLE Training Procedure The set of inputs
$\left\\{\mathbf{x}_{i}\right\\}_{i=1}^{n}$ and the set of corresponding
observed outputs $\left\\{\mathbf{y}_{i}\right\\}_{i=1}^{n}$ Initialize the
parameters $\theta$ of the generator $T_{\theta}$ $p=1$ to $N$ Pick a random
batch $S\subseteq\\{1,\ldots,n\\}$ $i\in S$ Randomly generate i.i.d. $m$
latent codes $\mathbf{z}_{1},\ldots,\mathbf{z}_{m}$
$\tilde{\mathbf{y}}_{i,j}\leftarrow
T_{\theta}(\mathbf{x}_{i},\mathbf{z}_{j})\;\forall j\in[m]$
$\sigma(i)\leftarrow\arg\min_{j}d(\mathbf{y}_{i},\tilde{\mathbf{y}}_{i,j})\;\forall
j\in[m]$ $q=1$ to $M$ Pick a random mini-batch $\tilde{S}\subseteq S$
$\theta\leftarrow\theta-\eta\nabla_{\theta}\left(\sum_{i\in\tilde{S}}d(\mathbf{y}_{i},\tilde{\mathbf{y}}_{i,\sigma(i)})\right)/|\widetilde{S}|$
$\theta$

##  4 Method

Overview. In the shape completion problem, we want to predict a complete point
cloud based on a given partial point cloud with missing regions or parts. In
the multimodal version of this problem, a set of complete shapes is desired
under only one given partial point cloud. Formally, given an input
$\mathbf{x}$ in a one-to-many problem setting, the goal is to learn
$p(\mathbf{y}|\mathbf{x})$ where $\mathbf{y}$ is a prediction. Here,
$\mathbf{x}$ is our partial point cloud and $\mathbf{y}$ is a predicted
complete point cloud. Since there are many possible complete shapes that are
consistent with $\mathbf{x}$, we want to generate non-deterministic
predictions $p(\mathbf{y}|\mathbf{x})$ instead of a single prediction. The
simplest way to generate non-deterministic predictions is to add a random code
input $\mathbf{z}\sim\mathcal{N}(0,\mathbf{I})$ to a deterministic function
$T_{\theta}$ where $\mathbf{y}\coloneqq T_{\theta}(\mathbf{x},\mathbf{z})$.
This is also known as an implicit generative model [mohamed2016learning]. One
way to train this implicit generative model is to use a conditional GAN, in
this context, where $T_{\theta}(\mathbf{x},\mathbf{z})$ is often known as
_generator_. However, the generator may fail to produce diverse outputs due to
mode collapse. Therefore, we use conditional Implicit Maximum Likelihood
Estimation (IMLE) to train the generative model to overcome the mode collapse
issue.

Figure 2: Autoencoder is trained on both partial (blue) and complete (purple)
shapes using Earth mover’s distance (a). Two losses are used on shapes and
latent codes in our network (b). Losses are highlighted in red.

IMLE reverses the direction of how generated samples are matched to the target
data: instead of making the generated samples similar to some target data,
IMLE ensures that every target data has a similar sample being generated. In
addition, IMLE does not require an adversarial objective which simplifies
training. If we denote the generator as $T_{\theta}(\cdot,\cdot)$, an input
point cloud as $\mathbf{x}_{i}$, the corresponding ground truth as
$\mathbf{y}_{i}$, and a latent code as $\mathbf{z}_{i,k}$, then IMLE would
optimize the following equation:

$\min_{\theta}\mathbb{E}_{\mathbf{Z}}\left[\sum_{i=1}^{n}\min_{k\in\\{1,\ldots,m\\}}d(T_{\theta}(\mathbf{x}_{i},\mathbf{z}_{i,k}),\mathbf{y}_{i})\right]$
(2)

where
$\mathbf{Z}={\mathbf{z}_{1,1},\ldots,\mathbf{z}_{n,m}\sim\mathcal{N}(0,\mathbf{I})}$,
$d(\cdot,\cdot)$ is a distance metric and, $m$ is a hyperparameter. Algorithm
4 shows the conditional IMLE training procedure. In our case, IMLE is
conditioned on the given partial shape.

[h] Conditional IMLE Training Procedure The set of inputs
$\left\\{\mathbf{x}_{i}\right\\}_{i=1}^{n}$ and the set of corresponding
observed outputs $\left\\{\mathbf{y}_{i}\right\\}_{i=1}^{n}$ Initialize the
parameters $\theta$ of the generator $T_{\theta}$ $p=1$ to $N$ Pick a random
batch $S\subseteq\\{1,\ldots,n\\}$ $i\in S$ Randomly generate i.i.d. $m$
latent codes $\mathbf{z}_{1},\ldots,\mathbf{z}_{m}$
$\tilde{\mathbf{y}}_{i,k}\leftarrow
T_{\theta}(\mathbf{x}_{i},\mathbf{z}_{k})\;\forall k\in[m]$
$j\leftarrow\arg\min_{k}d(\mathbf{y}_{i},\tilde{\mathbf{y}}_{i,k})\;\forall
k\in[m]$ $q=1$ to $M$ Pick a random mini-batch $\tilde{S}\subseteq S$
$\theta\leftarrow\theta-\eta\nabla_{\theta}\left(\sum_{i\in\tilde{S}}d(\mathbf{y}_{i},\tilde{\mathbf{y}}_{i,j})\right)/|\widetilde{S}|$
$\theta$

Shape completion network’s overall structure. The structure of our network is
composed of an auto-encoder (i.e., encoder $E$ and decoder $D$) along with
generator $G$ (see Figure 1). The auto-encoder encodes partial shape
$\mathcal{X}_{i}$ and complete shape $\mathcal{Y}_{i}$ via its encoder and
generates a novel and complete shape $\mathcal{\tilde{X}}_{i}$ via its
decoder. Generator $G$ plays the role of $T_{\theta}$ in Algorithm 4 by
mapping the latent code $\mathbf{x}_{i}$ of the partial shape, aggregated by a
noise vector $\mathbf{z}_{j}$ sampled from normal distribution, to a new code
$\mathbf{\tilde{y}}_{i,j}$. $\mathbf{\tilde{y}}_{i,j}$ is then passed to
decoder $D$ to produce a new shape $\mathcal{\tilde{X}}_{i,j}$. Therefore,
each partial shape $\mathcal{X}_{i}$ can be mapped to many shapes
$\mathcal{\tilde{X}}_{i,j}$.

Multi-modality and training. To train the network on the entire data, for each
input $\mathbf{x}_{i}$, $m$ forward passes with different random vectors
($\mathbf{z}_{1}$,$\mathbf{z}_{2}$,…,$\mathbf{z}_{m}$) are processed.
Therefore, the output from the generator is $m$ latent vectors
($\tilde{\mathbf{y}}_{i,1}$,$\tilde{\mathbf{y}}_{i,2}$,…,$\tilde{\mathbf{y}}_{i,m}$)
where each latent vector corresponds to a different complete shape or _mode_.
As also described in Algorithm 4, we then use a nearest neighbour algorithm to
select $\tilde{\mathbf{y}}_{i,j}$ which is closest to the latent vector
encoded from the corresponding complete shape $\mathbf{y}_{i}$. We compute an
$\text{L}_{1}$ loss between $\tilde{\mathbf{y}}_{i,j}$ and $\mathbf{y}_{i}$ to
bring $\tilde{\mathbf{y}}_{i,j}$ closer to the selected mode (see Figure 2):

$L_{\text{L}_{1}}=d(\mathbf{y}_{i},\tilde{\mathbf{y}}_{i,j})$ (3)

Enforcing that every complete shape in our training data has a _similar_ shape
in generated samples ensures that mode collapse does not occur and we have
diversity in generated shapes. Finally, we decode $\tilde{\mathbf{y}}_{i,j}$
through decoder $D$ and produce a novel complete shape. We use Unidirectional
Hausdorff Distance (UHD) loss between the given partial shape and the
generated complete shape to ensure that the generated complete shape encloses
the partial shape effectively. This loss can be defined as:

$L_{UHD}=\max_{x\in\mathcal{X}_{i}}\min_{y\in\mathcal{\tilde{X}}_{i,j}}||x-y||,$
(4)

where $x$ and $y$ are single points belonging to partial and the generated
complete point clouds.

Note that our auto-encoder is trained for reconstruction task on a point cloud
dataset which includes both partial and complete shapes. We employed Earth
Mover’s Distance (EMD) loss between input and output point clouds
$\mathbf{p}_{in}$ and $\mathbf{p}_{out}$ to train our auto-encoder:

$L_{EMD}(\mathbf{p}_{in},\mathbf{p}_{out})=\min_{\xi:\mathbf{p}_{in}\rightarrow\mathbf{p}_{out}}\sum_{\alpha\in\mathbf{p}_{in}}||\alpha-\xi(\alpha)||$
(5)

Testing. During testing, our model does not require any accompanied complete
shape $\mathcal{Y}_{i}$ because of which there is no nearest neighbour
selection involved. We simply pass the partial shapes along with different
randomly sampled noise vectors to the network and decode all the produced
latent codes into complete generated shapes.

##  5 Experiments and Results

In this section, we provide various experiments and compare our work both
qualitatively and quantitatively with baseline methods. In the end, we are
providing some statistics about timing along with qualitative results to show
how much our method respects the partial input shape. We have also provided
some complementary experiments on baselines and real-world scan data in the
supplementary material.

Baselines. We compare the result of our point completion tasks against two
multimodal methods: MSC [wu2020multimodal], and kNN-latent (a baseline
introduced by [wu2020multimodal]). kNN-latent is a simple baseline to look for
most similar complete shapes in dataset using encoded latent vector of partial
shape. MSC uses a conditional GAN architecture to generate multiple shapes. We
perform quantitative and qualitative comparisons with these baselines in the
following. In supplementary material, we provide quantitative and qualitative
comparisons with other methods including uni-modal shape completion techniques
such as pcl2pcl [chen2019unpaired] and VRCNet [pan2021variational]. We show
that our method performs comparably.

Datasets. We use three datasets to evaluate our work in different settings.
First, we use 3Depn dataset [dai2017complete] which is generated from ShapeNet
using virtual scans. The dataset has varied proportions of incompleteness in
the samples. We evaluate our method on three classes of this dataset - Chair,
Plane, and Table. Second, we evaluate Partnet dataset [chang2015shapenet] in
which parts are randomly removed from a complete shape while training/testing.
In Partnet, we evaluate on Chair, Table, and Lamp because of high number of
samples for these classes. Finally, we use the recently released
RobustPointSet [taghanaki2020robustpointset] dataset for evaluating
transferablity of our network. Shapes in this dataset have missing regions
synthetically created on various objects.

Evaluation. Here we provide quantitative and qualitative results to evaluate
our method against baselines. In Table 1 we compare our result using Total
Mutual Difference on 3Depn (left) and Partnet (right) datasets. Total Mutual
Difference (TMD) quantifies the diversity among the generated samples. This is
calculated using Chamfer distance between every pair among generated samples.
We show that we have significant gains in terms of diversity in generated
samples in comparison with multimodal baselines. The difference can be
explained by the fact that our work tends to use all modes in the dataset.
Similarly, in Table 2 we compare our result using Unidirectional Hausdorff
Distance (UHD) on 3Depn (left) and Partnet (right) datasets. UHD is a measure
of how well the partial shape is enclosed by the generated complete shape. We
either outperform the baselines or perform relatively similarly. To provide
qualitative results, we show our generated samples for various shape
categories in Figure 3. It can be observed that in the generated shapes, our
method tends to produce various geometric and shape attributes.

Methods | Chair | Plane | Table
---|---|---|---
kNN-latent | 2.84 | 1.13 | 3.25
Multimodal-GAN | 2.56 | 2.03 | 4.49
Ours | 2.93 | 2.31 | 4.92

Methods | Chair | Lamp | Table
---|---|---|---
kNN-latent | 2.28 | 4.18 | 2.85
Multimodal-GAN | 2.56 | 3.31 | 3.88
Ours | 2.76 | 5.49 | 4.45

Table 1: Comparison between our methods and baselines in terms of Total Mutual
Difference (TMD) or diversity on two datasets - 3Depn (left) and PartNet
(right). Higher TMD values indicate more diversity. TMD values in the table
are multiplied by a factor of $10^{2}$.

Methods | Chair | Plane | Table
---|---|---|---
kNN-latent | 8.94 | 9.54 | 10.04
Multimodal-GAN | 9.14 | 9.59 | 8.74
Ours | 8.51 | 9.55 | 8.52

Methods | Chair | Plane | Table
---|---|---|---
kNN-latent | 8.58 | 8.47 | 7.61
Multimodal-GAN | 6.65 | 5.40 | 5.38
Ours | 6.17 | 5.58 | 5.16

Table 2: Comparison between our method and baselines in terms of
Unidirectional Haussdorf Distance (UHD) or completeness on two datasets -
3Depn (left) and PartNet (right). Lower UHD values indicate better completion.
UHD values in the table are multiplied by a factor of $10^{2}$. Figure 3: Our
completion results for multiple object categories (e.g., Chair, Lamp, Plane
and Bottle). It can be observed that the generated complete shapes are diverse
for every object and respects partial shapes.

Figure 4 contrasts our results to conditional GAN [wu2020multimodal] and shows
that our model has more diversity. Note that our method respects the partial
input very well as apparent in Figure 5. To make a better contrast, we have
overlapped the partial input highlighted in green over the generated complete
shape highlighted in blue.

Figure 4: Qualitative Comparison of our work (left) and conditional GAN
[wu2020multimodal] (right). As it is apparent, our method produces more
diverse objects. Figure 5: Partial shapes (green) and generated complete
shapes (blue) are overlapped to illustrate how well our complete shape
respects the partial shape.

Part-based completion. In Figure 6, we demonstrate the ability of our method
to generalize to unseen types of partial shapes as input. We remove a part
entirely and visualize different completions generated by our method. Even
though no part labels were used during training and so the model has never
seen an entire part taken out, the model can generate reasonable and diverse
completions. Moreover, the model can generate a diverse range of parts in
place of the missing part that is functionally similar but structurally
different, which suggest several interesting properties: (1) the model has
learned which parts have the same function and are interchangeable, suggesting
it has acquired an implicit understanding of the semantics of the part in
question even though it has never seen part labels during training. (2) The
model does not necessarily generate the particular shape of the part that was
taken out in all completions, suggesting it did not collapse to a single mode.
(3) The model is able to generate structurally diverse shapes for missing
parts in different completions, suggesting that it has learned the space of
structural variations for a part. Moreover, the variations can be localized to
a part even though global noise is injected into the architecture.

Figure 6: Three partial portions are taken from complete shapes for two
categories (Chair and Lamp) and completed by our technique.

Transfer Learning. To examine whether our model can be used on other datasets,
possibly with more non-uniform missing parts and noises, we perform an
experiment on Robustpointset dataset [taghanaki2020robustpointset].
Robustpointset contains various artifacts including missing regions,
rotations, noise, etc. We use our model trained on 3Depn dataset and fine-tune
it on this dataset. Our model successfully fills in or completes the missing
regions and parts. In Figure 7, we show that our model can be successfully
transferred to other datasets. It is worth noting that these samples do not
have much scope for a diverse generation, so our model naturally converges to
the actual shape.

Figure 7: Generator trained on 3Depn dataset is fine-tuned on RobustPointSet
dataset and results are evaluated on Chair class.

Timing Comparison. Since our network does not have a discriminator, its
training and inference times are lower than conditional GAN. Our method needs
13 ms per shape for training and 2 ms for inference, while conditional GAN
needs 20ms per second and 2ms for inference. [wu2020multimodal]. On a Geforce
RTX 2080 GPU, the overall training time for conditional GAN work for the main
class Chair (3Depn) is 10 hours for the autoencoder, 10 hours for the VAE and
5 hours for the generator/discriminator network. In comparison, for the same
case, our work takes total 13.5 hours (10 hours for the autoencoder and 3.5
hours for the generator training).

Noise Tolerance. We perform noise tolerance analysis of our work by adding
Gaussian noise to input point clouds. We perform this on 3Depn dataset for
Chair category. As it can be seen in Figure 8 , we are able to achieve
reasonable completion despite having noise in partial shapes. We achieve
$2.71$ x $10^{-}2$ as Total Mutual Difference value and $8.84$ x $10^{-}2$ as
Unidirectional Hausdorff Distance value with noisy input point clouds.

Figure 8: Noisy partial shapes are fed through the network for shape
completion through our technique

Implementation details. We tune our hyperparameters through experiments to
have a standard set of parameters for all classes. For the encoder and decoder
model, we use 1D convolutional layers coupled with batch-normalization layers.
We train this auto-encoder network for 1000 epochs with initial learning rate
of $5e-4$ . Similarly, we use linear layers for the generator and train it for
500 epochs with the same initial learning rate. During testing, we set $m=10$
for the number of samples we generate for every partial shape.

Limitations. Our current network cannot generate very detailed features in a
multimodal manner. This is primarily due to the network’s simple architecture
and can be resolved by adding refinement networks.

##  6 Conclusion and Future Work

In this paper, we tackled the problem of multimodal shape completion problem,
in which a given partial shape can be completed in several different forms. In
fact, when a partial shape is missing a region or part, there is not a single
and fixed complete shape and the result can be subject to interpretation. As a
result, a method capable of producing several plausible results is desired to
expand the creativity and diversity of shape completion techniques. Here, we
used IMLE that tends to avoid mode collapse as opposed to conditional GAN to
make sure that we can deliver diverse completed shapes. Through several
experiments, we have shown that our method is superior in respecting the input
partial shapes’ geometry while producing more diverse results. We have also
shown that without direct supervision, the network learns to respect geometric
properties of shapes such as symmetry and it learns to generate missing parts
in the right location with accurate geometry. The results of our shape
completion are also comparable with the state-of-the-art uni-modal shape
completion techniques (see supplementary material), which shows that IMLE is
capable of respecting features when they are present.

Several avenues can be explored related to shape completion and also IMLE. Our
shape completion only adds parts and missing pieces. However, it would be
interesting to study how we can add styles and fine details in a multimodal
manner. IMLE can also be used in other tasks where a one-to-many mapping is
needed. An example can be single view reconstruction, where several objects
can be produced according to a single input image.


# Scalable Scheduling Policies for Quantum Satellite Networks

Albert Williams, Nitish K. Panigrahy, Andrew McGregor, and Don Towsley
University of Massachusetts Amherst, USA.
Email: {abwillia, nitish, mcgregor<EMAIL_ADDRESS>
###### Abstract

As Low Earth Orbit (LEO) satellite mega constellations continue to be deployed
for satellite internet and recent successful experiments in satellite-based
quantum entanglement distribution emerge, a natural question arises: _How
should we coordinate transmissions and design scalable scheduling policies for
a quantum satellite internet?_ In this work, we consider the problem of
transmission scheduling in quantum satellite networks subject to resource
constraints at the satellites and ground stations. We show that the most
general problem of assigning satellites to ground station pairs for
entanglement distribution is NP-hard. We then propose four heuristic
algorithms and evaluate their performance for Starlink mega constellation
under various amount of resources and placements of the ground stations. We
find that the maximum number of receivers necessary per ground station grows
very slowly with the total number of deployed ground stations. Our proposed
algorithms, leveraging optimal weighted b-matching and the global greedy
heuristic, outperform others in entanglement distribution rate, entanglement
fidelity, and handover cost metrics. While we develop these scheduling
algorithms, we have also designed a software system to simulate, visualize,
and evaluate satellite mega-constellations for entanglement distribution.

##  I Introduction

Long distant entanglement is a key resource for a variety of distributed
quantum applications, such as quantum key distribution (QKD) [1], quantum
sensing [2], and quantum teleportation [3]. However, superposition and
entanglements are susceptible to environmental noise and losses. In order to
distribute entanglements between distant users, ongoing efforts focus on
mitigating noise and increase robustness. Terrestrial entanglement
distribution invariably suffers loss that decays exponentially as a function
of distance [4] and must make use of quantum repeaters [5] for long distance
distribution. Locally created entangled particles can, instead, travel through
free space and incur loss that is only quadratic in the distance travelled [6,
7]. This makes satellite based entanglement distribution [8, 7, 9, 10] an
attractive alternative to long distance terrestrial links with quantum
repeaters.

Satellite based entanglement distribution, however, has its own set of
challenges. The rate and fidelity of entanglement distribution through free
space depend heavily on atmospheric conditions such as cloud cover, weather,
and time of day. Additionally, one needs to carefully select the optimal
altitude for satellite deployment. A higher orbit incurs greater loss, but
offers longer contact periods. Conversely, a lower orbit experiences lower
loss, but is more dynamic with shorter contact periods. For example,
Geostationary (GEO) satellites orbit at altitudes over $35,000$ kilometers, a
range where very few entanglements ($\sim 10^{-4}$) propagate through free
space, even though the decay rates are only quadratic. For Low Earth Orbit
(LEO) satellites, even in the best case, a single satellite rarely stays
within line of sight of a given ground station for more than five minutes.
However, the shorter contact periods can be offset by deploying multiple
satellites in constellations, ensuring a few of them remain visible to the
ground stations most of the time.

Recently, there has been a renewed interest in the deployment of LEO satellite
constellations for global scale broad band connectivity [11, 12, 13]. These
constellations, also known as mega constellations, consist of hundreds to
thousands of low cost communication satellites orbiting around the earth at
low-orbits. For example, SpaceX has already deployed $6000$ satellites into
its Starlink constellation [11]. While these satellites aim to provide
conventional classical data communication using RF based satellite-to-ground
and Laser based inter satellite links, their potential for enabling photonic
entanglement distribution has yet to be explored. Such mega constellations
hold promise for significantly improving the rate of satellite based
entanglement distribution while providing this service on a global scale.

In a satellite based entanglement distribution, entanglements are distributed
to ground station pairs which may be located far from each other. First, an
entangled pair of qubits is generated at a source located on the satellite.
Subsequently, these qubits are transmitted to each of the ground stations in
the pair through free space. Unlike classical satellite communication, for
successful entanglement distribution, the satellite must be simultaneously
visible to both the ground stations. We refer to this as the _pair visibility
requirement_.

A key aspect of mega constellations is the capability of multiple satellites
to maintain visibility with ground stations. Similarly, multiple ground
station pairs can be visible to a single satellite at any point of time.
Therefore to optimize the operation of such a set-up, i.e. a constellation of
satellites serving entanglements to a set of ground station pairs, one needs
to solve a scheduling problem - which satellite to distribute entanglements to
which ground station pairs and when to do it. This scheduling process must
align with the resource constraints present at both satellites and ground
stations. We refer to the scheduling problem as the _Quantum Satellite
Scheduling Problem_ (QSSP). One of the contributions of this work is to devise
scalable solutions to QSSP and subsequently evaluate these solutions in the
context of satellite mega constellations.

QSSP has been investigated by Panigrahy et. al. [6] with an emphasis on
formulating an optimal schedule to maximize the aggregate entanglement
distribution rate across all ground station pairs. The authors modeled QSSP as
an Integer Linear Program (ILP). While the ILP solution works well for a
modest number of ground station pairs and satellites, it does not scale when
confronted with a large number of ground station pairs and mega-
constellations. In our work, we prove that the general version of QSSP is NP-
hard, by reducing it from a known NP-hard problem. Our focus in this work is
on a broader setup where we identify essential features of a scalable and
efficient solution to QSSP. In particular, we summarize the major
contributions of our work as follows.

  * •

We show that the most general version of QSSP is NP-hard by reducing it from a
known NP-hard problem, the 3-dimensional matching problem.

  * •

We introduce four heuristic algorithms for QSSP and assess their performance
for Starlink mega-constellation across diverse scenarios, including various
placements and resource levels at ground stations, while factoring in
practical considerations such as time-of-day/year and realistic entanglement
generation sources. Our results demonstrate that algorithms based on optimal
weighted b-matching and global greedy heuristics outperform others in
different performance metrics.

  * •

We build a software system for simulation, visualization, and evaluation of
satellite mega-constellations for entanglement distribution.

##  II Preliminaries

In this Section, we present the system model used in the remainder of the
paper.

###  II-A Satellites and Ground Stations

Figure 1: The Dual Downlink architecture for entanglement distribution. Symbol |  Description
---|---
$S$ |  set of satellites
$T_{i}$ |  ground station pairs satellite $i$ can connect
$G$ |  set of ground stations
$R_{g}$ |  receivers on ground station $g\in G$
$F$ |  set of ground station pairs $\left(F\subset G\times G\right)$
$L_{j}$ |  maximum connections between $g_{1}$ and $g_{2}$,
$j=(g_{1},g_{2})$
$w_{ij}(t)$ |  value of connecting $g_{1}$, $g_{2}$ with $i$, $j=(g_{1},g_{2})$
$x_{ij}(t)$ |  decision variable: number of connections $i$ makes between $g_{1}$ and $g_{2}$, $j=(g_{1},g_{2})$, $x_{ij}\in\mathbb{Z}^{\geq 0}$
TABLE I: Summary of notation

Let $S$ denote the set of satellites in the constellation. For each satellite
$i\in S$, let $T_{i}$ denote the number of ground station pairs that satellite
$i$ can connect under ideal circumstances. Because each ground station pair
comprises two ground stations, $T_{i}$ is less than the number of transmitters
on satellite $i$. The exact relationship between $T_{i}$ and the number of
transmitters depends on the satellite’s design. Let $G$ be the set of ground
stations, and $R_{g}$ be the number of receivers at ground station $g$ for
each $g\in G$. Let $F\subset G\times G$ be the set of ground station pairs
that can be and desire to be connected, and let $L_{j}$ be the maximum number
of connections ground station pair $j\in F$ will accept between them. We
assume $L_{j}\leq\min\\{R_{g_{l}},R_{g_{m}}\\}$ for pair $j=(g_{l},g_{m}).$

###  II-B Dual Downlink Architecture

In this work, we consider a dual downlink architecture as shown in Figure 1.
Here, entangled pairs of qubits are generated by satellites and transmitted to
two receiving stations on the ground via free space. We assume satellites have
no onboard quantum memories, so the transmission of each qubit in an entangled
pair must occur immediately. This architecture, while pragmatic for its
limited use of quantum memory, is limited in that only pairs of ground
stations that are in view of the same satellite can be connected. Therefore
connections are limited by geographic distance111In our simulations, which use
low Earth orbit satellites and require that satellites be a minimum of 20
degrees above the horizon, the average satellite can directly connect ground
stations at most 2,400 kilometers apart., depending upon satellite location.
This limitation can be overcome with inter-satellite links which is out of the
scope of this work. In this work we focus on the problem of maximizing
entanglement distribution rate under the constraint that there are no inter-
satellite links. This is a useful metric in its own right, and serves as a
preliminary metric for evaluating the distribution rate of protocols using
entanglement swapping and inter-satellite links.

###  II-C Entanglement Generation Source

In our setup, bipartite entanglements are generated using a Spontaneous
Parametric Down-Conversion (SPDC) based dual-rail polarization entangled pair
source located on the satellite [6, 14]. We assume the generated entanglements
have the following form.

$\displaystyle\ket{\psi^{\pm}}$
$\displaystyle\\!=N_{0}\bigg{[}\\!\sqrt{p(0)}\ket{0,0;0,0}\\!$
$\displaystyle+\\!\sqrt{\frac{p(1)}{2}}\left(\ket{1,0;0,1}\pm\ket{0,1;1,0}\right)$
$\displaystyle+\sqrt{\frac{p(2)}{3}}\left(\ket{2,0;0,2}\pm\ket{1,1;1,1}+\ket{0,2;2,0}\right)\bigg{]},$
(1)

where $N_{0}$ is a normalization constant and

$p(n)=(n+1)\frac{{N_{s}}^{n}}{(N_{s}+1)^{n+2}}.$ (2)

Here, $N_{s}$ denotes the mean photon number per mode, also known as the pump
power. The source outputs a superposition state comprised of the following
three elements: (i) an entangled pair
($[\ket{1,0;0,1}\pm\ket{0,1;1,0}]/\sqrt{2}$), (ii) the vacuum
($\ket{0,0;0,0}$), and (iii) the unwanted two-pair photon terms
($[\ket{2,0;0,2}\pm\ket{1,1;1,1}+\ket{0,2;2,0}]/\sqrt{3}$). Note that,
choosing the right amount of pump power is important as higher pump power not
only enhances the rate of entanglement production but also increases the rate
of noisy two photon pair emissions.

###  II-D Loss and Noise

The entanglements generated at the source are split and transmitted to their
respective ground stations through free space. While undergoing transmission,
entanglements experience different forms of losses, with two significant ones
being - free space loss and atmospheric loss. Free space loss primarily occurs
due to diffraction and scales quadratically, whereas absorption and scattering
contribute to atmospheric loss, which grows as an exponential function of
distance. Unlike ground based transmissions where losses are exponential, the
majority of the losses in free space are the lower quadratic losses.

In this work. we consider two sources of noise - (i) Two pair photon
emissions, and (ii) Background photon flux. Due to two pair photon emissions,
losses in the channel can result in the loss of a single photon from two pair
photons, leading them to resemble entangled pairs. As these pairs are not
inherently entangled, they cause a decrease in fidelity of a distributed pair.
As stated earlier, two photon pair emissions can be controlled by choosing
appropriate pump power.

Background photons are being continuously generated from different light
sources (e.g. sun, moon, stars). They introduce noise, since it is difficult
to distinguish random and uncorrelated background photons from source-
generated entangled photons. We model the noise introduced by background
photons as an increase in the number of dark click counts at the detector at
the receiver. We treat the arrival of any background photon as a spurious
detection event registered by the receiver even though no logical photon is
present in the system. Typically, the background photon flux is orders of
magnitude higher during daytime compared to night time.

###  II-E Scheduling in a time snapshot

In this work, we assume that time is divided into slots, and the scheduling
problem is solved for each slot. This means that the satellite-to-ground
station pair-assignments, which are solutions to the scheduling problems
presented in this paper, are intended to be used only for a brief period. In
the next timeslot, the pair-assignments must be recalculated with new inputs.
This assumption simplifies the problem, and reflects likely network deployment
limitations in a few ways. In practice, the changes to the input parameters of
the problem (ambient conditions, such as weather) may occur at arbitrary times
and may involve only a small set of satellites and ground stations. In this
work, we do not model the continuous time version of the problem and instead
treat it as a discrete time system where scheduling changes occur at the
beginning of each time unit (slot).

Also, there may be a cost associated with such changes (e.g. delay due to
realignment of satellite transmitters or ground station receivers), which
might influence scheduling decisions. Any scheduling algorithm in deployment
should attempt to reduce the frequency with which communications equipment
will need to be aligned, a consideration that requires coordination across
multiple slots in time, something we do not address in this work. However, it
is worth noting that emerging entanglement generation sources, like the _zero-
added-loss multiplexing_ (ZALM) sources [15], can make the changes with very
low overhead. Designing scheduling policies that account for switching cost is
a topic of our future work.

##  III Quantum Satellite Scheduling Problem

We now formally introduce the Quantum Satellite Scheduling Problem (QSSP).

For each $i\in S$ and $j\in F$, let the decision variables $x_{ij}(t)$
represent the number of connections to be made between ground stations $g_{1}$
and $g_{2}$ via satellite $i$ at time slot $t$, for $j=(g_{1},g_{2})$. As the
number of connections is discrete and non-negative,
$x_{ij}(t)\in\mathbb{Z}^{\geq 0}$ for all $i\in S,\;j\in F$. Let $w_{ij}(t)$
be the acquired reward of connecting the ground station pairs in $j$ using
satellite $i$ at time $t$. This can be variously defined, but in our
simulations (described in section IV) we use the entanglement distribution
rate as the weight associated with a connection, which is calculated with
respect to atmospheric and free space distance between satellite $i$ and each
ground station in $j$. Note that if either ground station in $j$ is out of
range for satellite $i$, $w_{ij}(t)=0$.

In its most general form, the problem of scheduling transmissions from
satellites to ground station pairs with arbitrary resource/utility limits on
satellites, ground station pairs, and individual ground stations, can be
formalized as:

QSSP [6]: Given $S,F,$ and $G$, and $T_{i},L_{j},R_{g},$ and $w_{ij}(t)$ for
all $i\in S,j\in F,g\in G$ (i.e. values for the elements of Table I), find
values for $x_{ij}(t)$ at each time slot $t$ that

maximize $\displaystyle\sum_{i\in S,j\in F}w_{ij}(t)x_{ij}(t)$ (3) subject to
$\displaystyle\forall i\in S:$ $\displaystyle\sum_{j\in F}x_{ij}(t)\leq
T_{i},$ (4) $\displaystyle\forall j\in F:$ $\displaystyle\sum_{i\in
S}x_{ij}(t)\leq L_{j},$ (5) $\displaystyle\forall g\in G:$
$\displaystyle\sum_{i\in S}\sum_{j\in F|g\in j}x_{ij}(t)\leq R_{g},$ (6)
$\displaystyle\forall i\in S,\forall j\in F:$ $\displaystyle\quad
x_{ij}(t)\in\\{0,1,\cdots\\}.$ (7)

###  III-A NP-hardness of QSSP

The main result of this Section is as follows.

######  Theorem III.1.

The Quantum Satellite Scheduling Problem is NP-hard.

To show that QSSP is NP-hard, we give a reduction from a known NP-hard
problem: 3-Dimensional Matching.

3D Matching [16]: Let $G^{3D}=\left(V_{1}\cup V_{2}\cup V_{3},E\right)$ be an
undirected tripartite hypergraph such that $E\subseteq V_{1}\times V_{2}\times
V_{3}$. Find the maximum cardinality set of hyperedges $E^{*}\subseteq E$ such
that $(v_{1},v_{2},v_{3})\neq(v_{1}^{\prime},v_{2}^{\prime},v_{3}^{\prime})$
and $(v_{1},v_{2},v_{3}),(v_{1}^{\prime},v_{2}^{\prime},v_{3}^{\prime})\in
E^{*}\implies v_{1}\neq v_{1}^{\prime},v_{2}\neq v_{2}^{\prime},$ and
$v_{3}\neq v_{3}^{\prime}$.

For any arbitrary instance $G^{3D}\left(V_{1}\cup V_{2}\cup V_{3},E\right)$ of
the 3D Matching problem, we build an instance of QSSP222We remove the time
dependency, i.e., use $w_{ij},x_{ij}$ instead of $w_{ij}(t),x_{ij}(t)$ for
brevity., i.e., we construct an input:

$\displaystyle S,G,F,\left\\{T_{i}:i\in S\right\\},\left\\{R_{g}:g\in
G\right\\},\left\\{L_{j}:j\in F\right\\},$ $\displaystyle\left\\{w_{ij}:i\in
S,j\in F\right\\}$

to QSSP. Initially, let $S=F=G=\phi.$ For every hyperedge
$(v_{k},v_{l},v_{m})\in E$, update $S\leftarrow S\cup\\{i_{k}\\}$,
$F\leftarrow F\cup\\{j=(g_{l},g_{m})\\}$, and $G\leftarrow
G\cup\\{g_{l},g_{m}\\}$. Lastly, let $w_{ij}=T_{i}=L_{j}=R_{g}=1$ for all
$i\in S,j\in F,g\in G$.

Solving QSSP on this input gives a solution $x_{ij}\in\\{0,1\\}$ ($x_{ij}\leq
1$ since $\forall i\in S:T_{i}=1$, so $\forall i\in S,j\in
F:x_{ij}\leq\sum_{j^{\prime}\in F}x_{ij^{\prime}}\leq 1$). Let $E^{\prime}$ be
the set of hyperedges such that $(v_{k},v_{l},v_{m})\in E^{\prime}$ if and
only if $x_{i_{k}(g_{l},g_{m})}=1$. We aim to show that $E^{\prime}$ is a
solution to the 3D Matching problem on input $V_{1},V_{2},V_{3},E$.

First we show that $E^{\prime}$ is a matching of $G^{3D}$. Note that, because
$T_{i}=L_{j}=R_{g}=1$ for all $i\in S,j\in F,g\in G$ and by construction, we
have the following.

$\displaystyle(v_{k},v_{l},v_{m})\neq(v_{k}^{\prime},v_{l}^{\prime},v_{m}^{\prime})\text{
and }$
$\displaystyle(v_{k},v_{l},v_{m}),(v_{k}^{\prime},v_{l}^{\prime},v_{m}^{\prime})\in
E^{\prime}$ $\displaystyle\implies$
$\displaystyle(i_{k},j=(g_{l},g_{m}))\neq(i_{k}^{\prime},j^{\prime}=(g_{l}^{\prime},g_{m}^{\prime}))\text{
and }$ $\displaystyle x_{i_{k}j}=x_{i_{k}^{\prime}j^{\prime}}=1$
$\displaystyle\implies$ $\displaystyle i_{k}\neq i_{k}^{\prime},g_{l}\neq
g_{l}^{\prime},\text{ and }g_{m}\neq g_{m}^{\prime}$ $\displaystyle\implies$
$\displaystyle v_{k}\neq v_{k}^{\prime},v_{l}\neq v_{l}^{\prime},\text{ and
}v_{m}\neq v_{m}^{\prime}$

Now it remains to show that $E^{\prime}$ is a maximum cardinality matching.

Suppose that there exists some $E^{*}\subseteq E$ such that
$(v_{k},v_{l},v_{m})\neq(v_{k}^{\prime},v_{l}^{\prime},v_{m}^{\prime})$ and
$(v_{k},v_{l},v_{m}),(v_{k}^{\prime},v_{l}^{\prime},v_{m}^{\prime})\in
E^{*}\implies v_{k}\neq v_{k}^{\prime},v_{l}\neq v_{l}^{\prime},$ and
$v_{m}\neq v_{m}^{\prime}$, and $|E^{*}|>|E^{\prime}|$. Let $x^{*}_{i_{k}j}$
be defined as $x^{*}_{i_{k}j}=1$ when $(v_{k},v_{l},v_{m})\in E^{*}$, and
$x^{*}_{i_{k}j}=0$ when $(v_{k},v_{l},v_{m})\in E\setminus E^{*}$ with
$j=(g_{l},g_{m})$. Also, $\forall i_{k}\in S,\;j\in F,\text{ and }g\in G,$ we
have,

$\displaystyle\sum_{j^{\prime}\in
F}x^{*}_{i_{k}j^{\prime}}=|\\{j^{\prime}=(g_{l},g_{m})|(v_{k},v_{l},v_{m})\in
E^{*}\\}|\leq 1,$ $\displaystyle\sum_{i^{\prime}_{k}\in
S}x^{*}_{i^{\prime}_{k}j}(t)=|\\{i^{\prime}_{k}|(v_{k},v_{l},v_{m})\in
E^{*}\\}|\leq 1,$ $\displaystyle\sum_{i^{\prime}_{k}\in S}\sum_{j^{\prime}\in
F|g\in j^{\prime}}x^{*}_{i^{\prime}_{k}j^{\prime}}$
$\displaystyle=|\\{i^{\prime}_{k},j^{\prime}=(g_{l},g_{m})|(v_{k},v_{l},v_{m})\in
E^{*}\\}|\leq 1.$

Thus the vector: $[x^{*}_{ij},\forall i\in S,j\in F]$ satisfies the
constraints defined for the above instance of QSSP. But, we also have,

$\sum_{i\in S,j\in F}x^{*}_{ij}=|E^{*}|>|E^{\prime}|=\sum_{i\in S,j\in
F}x_{ij},$

which is a contradiction, since $[x_{ij},\forall i\in S,j\in F]$ is the
optimal solution of the defined instance of QSSP. Therefore, $\nexists
E^{*}:|E^{*}|>|E^{\prime}|$ and $E^{\prime}$ is a solution to the 3D Matching
problem. Because any arbitrary instance of the 3D Matching problem can be
transformed into an instance of QSSP in polynomial time, and 3D Matching is
NP-hard, QSSP is also NP-hard.

######  Remark 1.

It’s important to highlight that QSSP is more difficult than its classical
counterpart. The difficulty of this problem arises from the concept of
assigning satellites to ground station pairs, as opposed to individual ground
stations, while still facing resource constraints on individual ground
stations. Pair assignment, an uncommon feature in classical satellite
networks, adds complexity in the context of quantum satellite networks.

###  III-B Special Cases

While the most general version of QSSP is NP-hard, below we present a few
special cases for which polynomial time solutions exist.

####  III-B1 Satellites outnumbering ground stations

We assume a scenario where the satellites significantly outnumber the ground
stations. This scenario ensures the availability of a satellite for connecting
any two ground stations within each other’s range. Under this assumption, we
can solve the problem of pairing ground stations together before assigning
satellites to them. Therefore, removing the satellite transmitter constraints,
we get the following optimization problem.

$\displaystyle\text{maximize}\sum_{j\in F}w_{j}x_{j}$ (8)
$\displaystyle\text{subject to: }\forall j\in F:$ $\displaystyle\quad
x_{j}\leq L_{j},$ $\displaystyle\forall g\in G:$ $\displaystyle\sum_{j\in
F|g\in j}x_{j}\leq R_{g},$ $\displaystyle\forall j\in F:$ $\displaystyle\quad
x_{j}\in\\{0,1,\cdots\\}.$

Here, the weights $w_{j}$, represent increasing affine functions of the
distance between ground stations in pair $j.$ For this scenario, we can
construct a graph $U(V,E)$ where $V=G$ and $E=F$. We set $w_{j}$ and $L_{j}$
to be the weight and capacity of edge $j\in E,$ respectively. Define
$b(g)=R_{g}$ as the capacity of vertex $g\in V.$ Then solving problem (8) is
equivalent to solving the weighted capacitated b-matching problem for $U$. The
separation algorithm by Letchford et al. [17] can be used to solve Problem (8)
with a running time complexity of $O(|G|^{2}|F|\log(|G|^{2}/|F|)).$

####  III-B2 No resource constraints at ground stations

Such a scenario corresponds to removing the constraints on the number of
receivers (Constraint (6) in QSSP) at each ground station, or, equivalently,
setting $R_{g}=\infty$ for all $g\in G$. This can be accomplished by including
many receivers at each ground station which makes QSSP tractable. For this
case, we can construct a bipartite graph: $U(V,E)$ with bi-partitions $S$ and
$F$, such that $V=S\cup F$ and $E=C=\\{(i,j),\forall i\in S,\forall j\in
F,\text{ and }w_{ij}>0\\}.$ We set $w_{ij}$ as the weight of edge $(i,j)\in
E.$ Solving QSSP, when $R_{g}=\infty,$ is analogous to solving the
uncapacitated weighted b-matching problem for $U.$ QSSP, in this case, can be
solved using the modified minimum cost flow algorithm by Anstee et al. [18]
with a running time complexity of $O((|S|+|F|)^{2}|C|).$ This scenario
($R_{g}=\infty$) suggests a heuristic that we present in the next subsection.

Algorithm 1 Heuristic Algorithms for QSSP

1:procedure UpdateState($i,j,C$)

2: $x_{ij}\leftarrow x_{ij}+1$

3: $(g_{l},g_{m})=j$

4: $T_{i}\leftarrow T_{i}-1$; $L_{j}\leftarrow L_{j}-1$

5: $R_{g_{l}}\leftarrow R_{g_{l}}-1$; $R_{g_{m}}\leftarrow R_{g_{m}}-1$

6: if $T_{i}=0$ then $C\leftarrow C\setminus\\{(i,j^{\prime}),\forall
j^{\prime}\in F\\}$

7: if $L_{j}=0$ then $C\leftarrow C\setminus\\{(i^{\prime},j),\forall
i^{\prime}\in S\\}$

8: if $R_{g_{l}}=0$ then

9: $C\leftarrow C\setminus\\{(i^{\prime},j^{\prime}),\forall i^{\prime}\in
S,\forall j^{\prime}|g_{l}\in j^{\prime}\\}$

10: if $R_{g_{m}}=0$ then

11: $C\leftarrow C\setminus\\{(i^{\prime},j^{\prime}),\forall i^{\prime}\in
S,\forall j^{\prime}|g_{m}\in j^{\prime}\\}$

12:procedure RANDOM(C)

13: while $C\neq\phi$ do

14: Pick $(i,j)$ uniformly at random from $C$.

15: Call Procedure UpdateState($i,j,C$).

16:procedure LOCAL_GREEDY(C)

17: while $C\neq\phi$ do

18: $F^{\prime}=F\cap\\{j^{\prime}|(i^{\prime},j^{\prime})\in C,\forall
i^{\prime}\in S\\}$

19: Pick $j$ uniformly at random from $F^{\prime}$.

20: Pick $i\in S$ s.t. $(i,j)\in C$, $i$ maximizes $w_{ij}$.

21: Call Procedure UpdateState($i,j,C$).

22:procedure GLOBAL_GREEDY(C)

23: while $C\neq\phi$ do

24: Pick $(i,j)$ from $C$ that maximizes $w_{ij}$.

25: Call Procedure UpdateState($i,j,C$).

26:Initialize $C=\\{(i,j),\forall i\in S,\forall j\in F,\text{ and
}w_{ij}>0\\}$, $x_{ij}=0,\forall i\in S,j\in F$.

###  III-C Heuristic Algorithms

Since the most general version of QSSP is NP-hard and lacks a polynomial time
solution, we shift our attention to developing heuristic algorithms. We
consider three heuristic algorithms corresponding to different orderings of
the connections - (i) RANDOM, (ii) LOCAL_GREEDY, and (iii) GLOBAL_GREEDY. The
techniques are presented as procedures in Algorithm 1.

Each technique starts with an instance of all available connections $C$ (See
line 2 in Algorithm 1). Subsequently, a specific connection is selected based
on various criteria, each corresponding to a distinct heuristic rule. For
example, in the RANDOM algorithm, a connection is chosen uniformly at random
from $C$. In LOCAL_GREEDY, a ground station pair is first chosen uniformly at
random, and then a satellite is selected to establish the connection in a
manner that maximizes the entanglement distribution rate of the pair. In
GLOBAL_GREEDY, the connection that maximizes the entanglement distribution
rate among all connections in $C$ is selected. Upon selecting a connection,
the corresponding resources utilized to serve this connection are removed from
the available resource pool for all heuristic methods. If a satellite or a
ground station depletes all its resources, all of its connections are removed
from $C$. This procedure continues until $C$ becomes empty.

Algorithm 2 GREEDY_BACKOFF algorithm for QSSP

1:procedure GREEDY_BACKOFF(C)

2: $\mathtt{done\leftarrow False}$

3: while $\mathtt{done}$ is $\mathtt{False}$ do

4: Use modified minimum cost flow algorithm by Anstee et al [18] to find an
uncapacitated maximum weight b-matching with edge weights $w_{ij}$ to assign
$x_{ij}$ for all $i\in S,j=(g_{l},g_{m})\in F$ such that $i,j$ are not already
at capacity for $T_{i}$, $L_{j}$, $R_{g_{l}}$, and $R_{g_{m}}$.

5: $\mathtt{done\leftarrow True}$

6: $\mathtt{compliant\leftarrow False}$

7: while $\mathtt{compliant}$ is $\mathtt{False}$ do

8: $\mathtt{compliant\leftarrow True}$

9: if any constraints are violated then

10: $\mathtt{done\leftarrow False,compliant\leftarrow False}$

11: Remove the non-compliant connection with the lowest $w_{ij}$ value.

12:Initialize $C=\\{(i,j),\forall i\in S,\forall j\in F,\text{ and
}w_{ij}>0\\}$, $x_{ij}=0,\forall i\in S,j\in F$.

Additionally, we introduce a fourth heuristic method, GREEDY_BACKOFF, which
relies on an optimal algorithm for uncapacitated maximum weight b-matching
(Anstee et al algorithm [18]) between satellites and ground station pairs.
This is outlined in Algorithm 2. The algorithm involves initially establishing
a b-matching between satellites and ground stations, with edge weights
corresponding to entanglement distribution rate of the connections. Any
violation of the ground station receiver constraint (Constraint (6) in QSSP)
is addressed by removing connections with the lowest entanglement distribution
rate. The process is repeated until all resources are exhausted, and no more
connections remain to be served. Note that GREEDY_BACKOFF is optimal when
$R_{g}=\infty,\forall g\in G.$

##  IV Performance Evaluation

(a)

(b)

(c)

Figure 2: Satellites and ground station distributions, $|S|=3967,|G|=100.$ (a)
Satellites in the Starlink LEO mega-constellation. (b) Ground stations are
placed uniformly at random on the land. (c) Ground stations are chosen based
on highest population density.

In this section, we evaluate the solutions generated by various heuristics and
investigate the impact of factors such as satellite and ground station density
and ground station placement on their performance. We test the performance of
these heuristics using the Starlink LEO mega-constellation as a case study.
The Starlink constellation is a leading candidate among all mega-
constellations intended for classical satellite internet deployment. It yields
very large scale instances of QSSP, and the resulting instances cannot achieve
optimal solutions due to the computational hardness of QSSP. We detail the
Starlink case in the following section.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 3: (a) Satellite altitudes and inclination angles in the Starlink LEO
mega-constellation. (b) Maximum number of satellites visible per ground
station pair over time. (c) Maximum number of ground station pairs visible per
satellite over time. (d) Entanglement distribution rates achieved by different
heuristics. (e) Average number of connections per ground station for different
heuristics. (f)-(g) Effect of number of receivers on performance. (h)-(i)
Longevity of the assigned connections. (j)-(k) Average fidelity vs number of
receivers across day time and night time. (l) Effect of different time of the
year on entanglement fidelity.

###  IV-A The Starlink constellation

The Starlink constellation considered in this work consists of $3967$
satellites, distributed across four different orbital inclinations. Figure
3(a) illustrates the altitudes and orbital inclination angles of these
satellites. While the majority orbit at an altitude of $550$ km, there are
instances where satellites are positioned at lower altitudes, with a minimum
altitude of approximately $250$ km. The constellation includes four distinct
orbital inclination angles: $43,53,70,$ and $97$ degrees. At $550$ km
altitude, the satellites take around $100$ minutes to complete a full
revolution around the earth. To track the positions and trajectories of these
satellites, we utilize a two-line element set (TLE) file (available in our
code repository). Figure 2(a) displays the spatial distribution of the
satellites above Earth’s surface. Note that, the density of the satellites is
higher around altitudes of approximately $52$ degrees North and South due to
orbital characteristics, whereas density decreases when not between those
latitudes.

###  IV-B Ground Stations

Param | Value | Param | Value | Param | Value
---|---|---|---|---|---
$\lambda_{s}$ |  $737$ nm | $N_{S}$ | $0.078$ | $\tau$ | $10^{9}$
$\eta_{s}$ | $0.707$ | $\eta_{g}$ | $0.707$ | $r_{s}$ |  $0.1$ m
$r_{g}$ |  $1$ m | $t_{A}$ |  $5$ km | $\theta_{e}$ |  $20$ Deg
TABLE II: Simulation Parameters: $\lambda_{s},N_{S}:$ Wavelength and pump
power of SPDC source; $\tau:$ Repetition rate of SPDC source;
$\eta_{s},\eta_{g}:$ Inefficiencies at satellite transmitters and ground
station receivers; $r_{s},r_{g}:$ Radii of circular apertures of the
transmitter and receiver telescopes; $t_{A}:$ Thickness of the atmospheric
shell; $\theta_{e}:$ Elevation angle limit.

We examine two primary placement configurations for $100$ ground stations
defined as follows.

  * •

RANDOM_ON_LAND: In this configuration, ground stations are uniformly placed at
random on land. We first choose a latitude and a longitude uniformly at
random. We then resample if the chosen coordinates do not fall on land. Figure
2(b) illustrates the spatial distribution of the ground stations in this
configuration.

  * •

POPULATION_CENTERS: Here, we select ground stations corresponding to the $100$
most populous population centers. These centers are chosen from a list
maintained at [19]. A portion of the ground station distribution is depicted
in Figure 2(c).

We assume that for a satellite to be visible by a particular ground station,
its elevation angle must be above $\theta_{e}=20$ degrees above the horizon.
Consequently, each satellite at an altitude of $550$ km has a footprint with
an approximate radius of $1125$ km. This also suggests that ground stations
exceeding $2250$ km apart will not be able to establish direct connections via
a Starlink satellite positioned at $550$ km above the earth’s surface.

###  IV-C Simulation Setup

The simulation parameters are outlined in Table II. To calculate the weights,
$w_{ij}$, we use the quadratic free space and exponential atmospheric loss
model as detailed in Section II. From this, we compute transmissivities, which
in turn, are utilized to determine entanglement distribution rates based on
the distances between satellite $i$ and the ground stations in pair $j$, along
with an atmospheric thickness of $5$ km above ground level. We consider three
different $24$ hour periods: the two solstices (June 21 and December 22) and
one equinox (September 23). We divide each $24$ hour interval into $1$ minute
intervals. We divide each day into a night period and a day period. We set the
values of detector dark click probability ($P_{d}$) to $3\times 10^{-3}$ and
$3\times 10^{-7}$ during day and night time respectively [20]. At the
beginning of each time interval, we use the four heuristic algorithms to solve
QSSP. Our simulation is coded in Python 3 [21], and leverages third party
packages such as the PyEphem [22], Matplotlib [23], Networkx [24], NumPy [25],
and global-land-mask [26].

###  IV-D Satellite and ground station pair visibility

We investigate the visibility of satellites and ground station pairs in
Figures 3(b) and (c). The temporal progression of the peak number of
satellites visible to a pair of ground stations is presented in Figure 3(b).
It is interesting to observe that when ground stations are randomly
distributed on land, ground station pairs observe fewer satellites compared to
when ground stations are placed in highly populous centers, due to the latter
resulting in more geographic clusters of ground stations. We plot the maximum
number of ground station pairs visible to a satellite in Figure 3(c). The
results demonstrate that the Starlink constellation configuration suits more
to population centers even in the quantum case with dual ground station pair
visibility requirement. We observe that up to $50$ satellites are visible to a
single ground station pair during some of the time slots, while there are
times that $40$ ground station pairs are visible to a satellite. This
motivates the importance of scheduling, given the large number of potential
combinations of satellite to ground station pair assignments.

Figure 4: Ground stations in Asia connected with LOCAL_GREEDY, and
GREEDY_BACKOFF. Red links are the connections assigned only by GREEDY_BACKOFF,
blue links are assigned only by LOCAL_GREEDY, and black links are assigned by
both.

###  IV-E Performance of different heuristics

In Figure 3(d), we plot the average entanglement distribution rates achieved
across different heuristics and ground station placements. We find that rates
are lower when ground stations are randomly positioned compared to placements
based on population across all heuristic policies. The reason is evident from
the visibility plots (Figure 3 (b) and (c)). There are more satellites visible
per pair and more ground station pairs visible per satellite in the population
center case. So it easier to resolve conflicts, particularly when $R_{g}=1.$
Additionally, we observe that GLOBAL_GREEDY, and GREEDY_BACKOFF perform the
best, achieving nearly identical entanglement distribution rates regardless of
ground station placement.

We observe that LOCAL_GREEDY performs the worst compared to all other
algorithms, even falling below RANDOM. Upon further investigation, we observe
that LOCAL_GREEDY tends to prioritize either nearby connections or distant
connections, unlike other algorithms (e.g. GREEDY_BACKOFF) which tend to
choose connections of relatively medium distance. This is evident from Figure
4 where we show connections established by LOCAL_GREEDY and GREEDY_BACKOFF for
a single time interval between ground stations in Asia, when positioned at
population centers. LOCAL_GREEDY tends to choose the best and the worst links,
while GREEDY_BACKOFF opts for average distance connections. This relatively
egalitarian assignment in GREEDY_BACKOFF leads to substantial improvements in
entanglement distribution rate over LOCAL_GREEDY, as entanglement distribution
rate is a decreasing convex function of distance.

###  IV-F Effect of number of receivers

We study the sensitivity of the performance of different heuristics algorithms
to the number of receivers for population center based placement in Figure
3(f). We observe up to $70\%$ improvement in performance when transitioning
from $R_{g}=1$ to $R_{g}=5.$ When $R_{g}=5,$ the performance is very similar
to that of a scenario with an infinite number of receivers ($R_{g}=\infty$).
This is evident from the number of connections per ground station plot with
$R_{g}=\infty$ shown in Figure 3(e). When ground stations are placed at
population centers, they are form clusters, such as in East Asia and Europe.
With more ground station pairs competing for the same limited number of
satellites within range of that cluster, most ground stations do not need more
than five receivers. Note that, in Figure 3(e), a good portion of the ground
stations in LOCAL_GREEDY are assigned to more than five receivers. This
explains the performance gap for LOCAL_GREEDY between $R_{g}=5$ and
$R_{g}=\infty$ in Figure 3(f).

We also explore the impact of the number of receivers on performance of
RANDOM_ON_LAND placement, as shown in Figure 3(g). We observe that the
performance gap between $R_{g}=1$ and $R_{g}=5$ is further increased.
Moreover, the performance gap between $R_{g}=5$ and $R_{g}=\infty$ is also
increased.

###  IV-G Longevity of connections

We track the duration of connections between given ground station pair and a
given satellite for each heuristic algorithm in the $R_{g}=1$ setting,
illustrated in Figures 3 (h) and (i). We note that the best-performing
heuristic algorithm, GREEDY_BACKOFF, maintains more connections for multiple
timeslots than other algorithms do, with the worst-performing algorithm
LOCAL_GREEDY, maintaining fewer connections for multiple timeslots than any
other algorithm. This suggests that GREEDY_BACKOFF’s performance may improve
further when accounting for any potential overhead incurred when re-aligning
satellite transmitters and receivers to form new connections.

###  IV-H Effect of time of day and year

We plot the average fidelity of the distributed entanglements as a function of
number of receivers, during day and night times as shown in Figures 3(j) and
(k). We observe that increasing the number of receivers results in a slight
decrease in average fidelity across most heuristics. LOCAL_GREEDY has the
lowest fidelity, while GREEDY_BACKOFF and GLOBAL_GREEDY perform the best. The
day time fidelities are considerably lower compared to nighttime fidelities
across all heuristics and ground station placements. Infidelities caused
during the daytime are primarily due to the presence of background photons
emitted from sun, in addition to infidelities caused by two-photon pairs
emitted at the source during both day and night.

We also observe that daytime fidelities are lower when ground stations are
randomly placed on land versus when they are placed in population centers.
This is due to the shorter average distances between ground station pairs in
the latter case. We plot the achieved entanglement fidelity across different
times of the year for different heuristics in Figure 3(l). When ground
stations are placed in populous centers, the majority of them are in the
northern hemisphere. So the average fidelity is smaller for the June solstice
due to the extended daylight hours. Maximum fidelity is observed during the
December solstice, which has the longest night in northern hemisphere.

##  V Conclusion and Future directions

In this work, we analyzed QSSP for a LEO satellite mega-constellation with a
large set of ground stations. Via NP-hardness reduction, we proved that the
general version of QSSP is intractable. We ran extensive simulations to
explore the effectiveness of four heuristic algorithms to find a schedule for
the Starlink mega-constellation. Going forward, we would like to simulate more
diverse hardware, such as satellites with more transmitters, satellites
capable of establishing inter-satellite links, or even satellites with onboard
quantum memories. We would also like to investigate the constellation design
problem for quantum satellites which would involve determining the optimal
arrangement and configuration of satellites in space for entanglement
distribution.

## References

  * [1] A. K. Ekert, “Quantum cryptography based on bell’s theorem,” _Phys. Rev. Lett._ , vol. 67, pp. 661–663, Aug 1991. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett.67.661
  * [2] C. Degen, F. Reinhard, and P. Cappellaro, “Quantum sensing,” _Reviews of Modern Physics_ , vol. 89, no. 3, jul 2017. [Online]. Available: https://doi.org/10.1103%2Frevmodphys.89.035002
  * [3] X.-S. Ma, T. Herbst, T. Scheidl, D. Wang, S. Kropatschek, W. Naylor, B. Wittmann, A. Mech, J. Kofler, E. Anisimova, V. Makarov, T. Jennewein, R. Ursin, and A. Zeilinger, “Quantum teleportation over 143 kilometres using active feed-forward,” _Nature_ , vol. 489, no. 7415, pp. 269–273, sep 2012. [Online]. Available: https://doi.org/10.1038%2Fnature11472
  * [4] S. Pirandola, R. Laurenza, C. Ottaviani, and L. Banchi, “Fundamental limits of repeaterless quantum communications,” _Nature Communications_ , vol. 8, no. 1, apr 2017. [Online]. Available: https://doi.org/10.1038%2Fncomms15043
  * [5] W. Dür, H.-J. Briegel, J. I. Cirac, and P. Zoller, “Quantum repeaters based on entanglement purification,” _Physical Review A_ , vol. 59, no. 1, pp. 169–181, jan 1999. [Online]. Available: https://doi.org/10.1103%2Fphysreva.59.169
  * [6] N. K. Panigrahy, P. Dhara, D. Towsley, S. Guha, and L. Tassiulas, “Optimal entanglement distribution using satellite based quantum networks,” 2022.
  * [7] C. Y. Lu, Y. Cao, C. Z. Peng, and J. W. Pan, “Micius quantum experiments in space,” _Reviews of Modern Physics_ , vol. 94, no. 3, 2022.
  * [8] M. Gündoğan, J. S. Sidhu, V. Henderson, L. Mazzarella, J. Wolters, D. K. Oi, and M. Krutzik, “Proposal for space-borne quantum memories for global quantum networking,” _npj Quantum Information_ , vol. 7, no. 1, pp. 1–11, 2021.
  * [9] L. de Forges de Parny, O. Alibart, J. Debaud, S. Gressani, A. Lagarrigue, A. Martin, A. Metrat, M. Schiavon, T. Troisi, E. Diamanti, P. Gélard, E. Kerstel, S. Tanzilli, and M. Van Den Bossche, “Satellite-based quantum information networks: use cases, architecture, and roadmap,” _Communications Physics_ , vol. 6, no. 1, pp. 1–17, 2023.
  * [10] S. Khatri, A. J. Brady, R. A. Desporte, M. P. Bart, and J. P. Dowling, “Spooky action at a global distance: analysis of space-based entanglement distribution for the quantum internet,” _npj Quantum Inf_ , vol. 7, no. 4, 2021.
  * [11] N. K. Panigrahy, P. Dhara, D. Towsley, S. Guha, and L. Tassiulas, “Spacex starlink,” 2017.
  * [12] A. F. M. O. M. A. AUTHORIZATION, “Spacex starlink,” 2018.
  * [13] A. Boyle, “Amazon to offer broadband access from orbit with 3,236 satellite ‘project kuiper’ constellation,” 2019.
  * [14] P. Dhara, S. J. Johnson, C. N. Gagatsos, P. G. Kwiat, and S. Guha, “Heralded multiplexed high-efficiency cascaded source of dual-rail entangled photon pairs using spontaneous parametric down-conversion,” _Physical Review Applied_ , vol. 17, no. 3, p. 034071, 2022.
  * [15] K. C. Chen, P. Dhara, M. Heuck, Y. Lee, W. Dai, S. Guha, and D. Englund, “Zero-added-loss entangled-photon multiplexing for ground- and space-based quantum networks,” _Phys. Rev. Appl._ , vol. 19, p. 054029, May 2023. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevApplied.19.054029
  * [16] R. M. Karp, _Reducibility among Combinatorial Problems_. Boston, MA: Springer US, 1972, pp. 85–103. [Online]. Available: https://doi.org/10.1007/978-1-4684-2001-2_9
  * [17] A. N. Letchford, G. Reinelt, and D. O. Theis, “A faster exact separation algorithm for blossom inequalities,” _Proceedings of IPCO, LNCS 3064_ , pp. 196–205, 2004.
  * [18] R. P. Anstee, “A polynomial algorithm for b-matchings: An alternative approach,” _InformationProcessing Letters_ , vol. 24, no. 3, pp. 153–157, 1987.
  * [19] C. Youderain, “World cities database,” Jun 2021. [Online]. Available: https://simplemaps.com/data/world-cities
  * [20] C. Harney, A. I. Fletcher, and S. Pirandola, “End-To-End Capacities of Hybrid Quantum Networks,” _Physical Review Applied_ , vol. 18, no. 1, p. 1, 2022. [Online]. Available: https://doi.org/10.1103/PhysRevApplied.18.014012
  * [21] Python Core Team, _Python: A dynamic, open source programming language_ , Python Software Foundation, 2019. [Online]. Available: https://www.python.org/
  * [22] B. Rhodes, Jun 2021. [Online]. Available: https://rhodesmill.org/pyephem/index.html
  * [23] J. D. Hunter, “Matplotlib: A 2d graphics environment,” _Computing in Science & Engineering_, vol. 9, no. 3, pp. 90–95, 2007.
  * [24] A. Hagberg, P. Swart, and D. S Chult, “Exploring network structure, dynamics, and function using networkx,” Los Alamos National Lab.(LANL), Los Alamos, NM (United States), Tech. Rep., 2008.
  * [25] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy,” _Nature_ , vol. 585, no. 7825, pp. 357–362, Sep. 2020. [Online]. Available: https://doi.org/10.1038/s41586-020-2649-2
  * [26] T. Karin, “toddkarin/global-land-mask: Release of version 1.0.0,” Oct 2020.