# DCoflow: Deadline-Aware Scheduling Algorithm for Coflows in Datacenter
Networks

Quang-Trung Luu1, Olivier Brun1, Rachid El-Azouzi2, Francesco De Pellegrini2,
Balakrishna J. Prabhu1, Cédric Richier2

{qtluu, brun<EMAIL_ADDRESS>{rachid.elazouzi, francesco.de-pellegrini,
<EMAIL_ADDRESS>1LAAS-CNRS, University of Toulouse, CNRS,
31400 Toulouse, France  2CERI/LIA, University of Avignon, 84029 Avignon,
France

###### Abstract

Datacenter networks routinely support the data transfers of distributed
computing frameworks in the form of coflows, i.e., sets of concurrent flows
related to a common task. The vast majority of the literature has focused on
the problem of scheduling coflows for completion time minimization, i.e., to
maximize the average rate at which coflows are dispatched in the network
fabric. Modern applications, though, may generate coflows dedicated to online
services and mission-critical computing tasks which have to comply with
specific completion deadlines. In this paper, we introduce $\mathtt{DCoflow}$,
a lightweight deadline-aware scheduler for time-critical coflows in datacenter
networks. The algorithm combines an online joint admission control and
scheduling logic and returns a $\sigma$-order schedule which maximizes the
number of coflows that attain their deadlines. Extensive numerical results
demonstrate that the proposed solution outperforms existing ones.

###### Index Terms:

Time-sensitive coflow scheduling, coflow admission control, $\sigma$-order,
deadline.

##  I Introduction

Modern traffic engineering in datacenter networks is based on the notion of
coflow originally defined in [1]. The interest for this traffic abstraction
has originally been motivated by the need to capture the structure of the data
exchanges occurring in distributed computing frameworks such as MapReduce or
Spark [2, 3]. Such software frameworks rely on the so called _dataflow_
computing model for large-scale data processing, i.e., a distributed computing
paradigm, where each intermediate computation stage is distributed over a set
of nodes and its output is transferred to nodes hosting the next stage. In
between two computation stages, such dataflows are producing a set of network
flows traversing the datacenter fabric and are abstracted as a _coflow_. A
popular example of such data transfer is the _shuffle phase_ of Hadoop
MapReduce [2].

The customary performance metric for the data transfer phase is the makespan
or the weighted coflow completion time (CCT). Minimizing the average CCT is an
appropriate goal in order to increase the number of computing jobs dispatched
per hour in a datacenter. The weighted CCT minimization has thus been
addressed in several works, e.g., [1, 4, 5, 6]. A decade’s research on the
problem has shed light on its complexity and several algorithmic solutions
have been devised. The problem was proven NP-hard and inapproximable below a
factor of $2$ by reduction to the job scheduling problem on multiple
correlated machines. Near-optimal algorithms with $4$-approximation
performance bounds have been proposed in the literature [6, 7, 5]. However,
the context changes radically in the case of time-sensitive jobs, where the
data transfer phase may be subject to strict coflow deadlines.

Here, coflow scheduling is typically combined with admission control in order
to reduce the number of violations, i.e., the number of coflows to complete
after their deadlines. The resulting problem is the coflow deadline
satisfaction (CDS) problem introduced in [8]. Each coflow is subject to a
completion deadline and the target is to operate joint coflow admission
control and scheduling in such a way to maximize the number of admitted
coflows which respect their deadlines. This problem is NP-hard as well and it
is proved inapproximable within any constant factor from the optimum [8].

Even though the problem has been identified quite early in the literature [9],
with a few exceptions, the vast majority of works on coflow scheduling have
not dealt with the problem of time-sensitive coflows. On the other hand, as
confirmed later in our performance analysis, even near-optimal algorithms for
CCT minimization may fail to respect the coflow deadline. In reality, the
notion of time-sensitive coflows has become pervasive in the way how modern
datacenters operate as distributed networks. In fact, not only computing
frameworks are often tasked with time-sensitive jobs: modern web and mobile
applications are implemented using microservice architectures, so that the
users requests issued to an application may activate hundreds or thousands of
services from as many servers to retrieve the users data. The last incoming
bit of data, i.e., the CCT of this batch of flows, determines the lag to the
service response, and large delays degrade the quality of experience.

In this paper, we address the problem of maximizing the Coflow Acceptance Rate
(CAR), in which each coflow is subject to a completion time deadline. In
principle, one can solve this problem by formulating a suitable Mixed Integer
Linear Program (MILP). However, in datacenters with tens of thousands of
coflows [4], techniques based on MILPs or their relaxations may be not viable.
For designing scalable algorithms, the main idea appearing in many research
works is to schedule coflows using a priority order of coflows. Once an
ordering, denoted as $\sigma$, is determined, it is enough to adopt a work-
conserving transmission policy. We focus on $\sigma$-order schedulers since
they offer a key implementation advantage: at the level of rate control, any
work conserving preemptive dynamic rate allocation is allowed as long as it is
compatible with the input coflow priority (the maximal performance loss within
said rate allocation policies is bounded by a factor of $2$ [6]). For
instance, using fixed coflow priorities under DiffServ satisfies the
definition of a $\sigma$-order scheduler. On the other hand, in order to avoid
per-flow rate control, commercial switches have built-in priority queues and
per-flow tagging can be used to prioritize active coflows. In principle, this
permits to perform a greedy rate allocation which is compatible with a target
$\sigma$-order. The exact mapping from a coflow $\sigma$-order to the switch
priority queuing mechanism (and the inevitable limitations of legacy hardware
therein) while an interesting subject, is out of the scope of the present
paper.

Contributions. This paper proposes a lightweight method to perform coflow
scheduling under deadlines. The proposed solution provably outperforms
existing ones in the literature and does not rely on the solution of a linear
program. In particular, we propose first a baseline offline admission control
policy which is combined with a scheduler drawn in the class of $\sigma$-order
coflow schedulers [6]. The output of the algorithm is an order of priority
restricted to the set of admitted coflows. The algorithm is hence extended to
perform joint admission control and scheduling in the online scenario, when
coflows are generated at runtime at unknown release times. Through extensive
numerical experiments on a wide variety of scenarios, we show that our
algorithm systematically outperforms the ones currently available in the
literature. These experiments are performed on both offline and online setting
using both synthetic traces and real traces obtained from the Facebook data
[9]. The main observation is that, under higher workloads (i.e., when the
acceptance ratio is lower), our algorithm outperforms significantly the
existing ones. Thus, it proves robust to workload variations as well to the
type of data set used to generate the coflows.

The rest of the paper is organized as follows. Sec. II describes the general
problem tackled in the paper and the coflow ordering models, whereas Sec. III
describes the proposed algorithms. Numerical results are then provided in Sec.
IV. In Sec. V, we describe the literature on deadline-aware coflow scheduling.
Concluding remarks and future research directions are given in Sec. VI.

##  II Problem Statement and Theoretical Analysis

In this section, we present the system model and formulate the acceptance rate
maximization problem for a given input set of coflows. Table I summarizes the
main notations used throughout the paper. The datacenter network (or
datacenter fabric) is represented as a Big-Switch model [9], a non-blocking
switch whose ingress (egress) port $\ell$ has capacity $B_{\ell}$ equal to the
corresponding bandwidth inbound (outbound) capacity to connect servers to the
top-of-rack (ToR) switch. Due to large bisection capacity and customary usage
of load balancing, in fact, traffic congestion is typically observed only at
the rack access ports leading to the ToR switches.

We consider a batch of $N$ coflows ${\cal C}=\\{1,2,...,N\\}$. A coflow is a
collection of flows, in which each flow represents a shuffle connection over a
pair of fabric ingress-egress ports. Denote ${\cal F}_{k}$ as the set of flows
of coflow $k\in{\cal C}$ and assume that the volume $v_{k,j}$ of each flow
$j\in{\cal F}_{k}$ is known. For the sake of clarity, we suppose that all
coflows arrive at the same time, i.e., their release time is zero. Also, each
coflow $k$ is subject to a completion deadline $T_{k}$. Similarly, we let
${\cal F}_{\ell,k}$ be the set of flows in ${\cal F}_{k}$ which uses port
$\ell\in{\cal L}$ either as ingress port or as egress port. The total volume
of data sent by coflow $k$ on port $\ell$ is then given by
$\hat{v}_{\ell,k}=\sum_{j\in{\cal F}_{\ell,k}}v_{k,j}$. Let $p_{\ell,k}$ be
the transfer completion time in isolation, i.e., at full rate, of coflow $k$
at port $\ell$, which is given by $p_{\ell,k}=\hat{v}_{\ell,k}/B_{\ell}$. Let
$r_{k,j}(t)\in\mathbb{R}_{+}$ be the rate allocated to flow $j\in{\cal F}_{k}$
at time $t$. The CCT of coflow $k$, denoted as $c_{k}$, thus writes

$c_{k}=\max_{j\in{\cal F}_{k}}\frac{v_{k,j}}{\bar{r}_{k,j}},\mbox{
where}~{}\bar{r}_{k,j}=\frac{1}{CT_{k,j}}\int_{0}^{CT_{k,j}}r_{k,j}(t)\mathop{}\\!\mathrm{d}{t},$
(1)

and $\bar{r}_{k,j}$ is the average rate of flow $j$ through its lifetime and
$CT_{k,j}$ is its completion time. This quantity rules the dispatching time
for the volume traversing the so-called coflow bottleneck which determines the
CCT of the said coflow. Hence a coflow $k$ satisfies the deadline when the
last flow of the coflow finishes before $T_{k}$, i.e., $c_{k}\leq T_{k}$.

Let $z_{k}\in\\{0,1\\}$ be an indicator of whether coflow $k$ finishes before
$T_{k}$. Maximizing the number of coflows meeting their deadline corresponds
to maximizing the sum of $z_{k}$ subject to the constraint of bandwidth
capacity at ingress and egress ports. The CAR maximization problem can be
formulated as

$\displaystyle\underset{r}{\mathrm{max}}\enskip$
$\displaystyle\sum_{k\in\mathcal{C}}z_{k}$ (P1)
$\displaystyle\mathrm{s.t.}\enskip$
$\displaystyle\sum_{k\in\mathcal{C}}\sum_{j\in\mathcal{F}_{k,l}}r_{k,j}(t)\leq
B_{\ell},\quad\forall\ell\in\mathcal{L},\forall t\in\mathcal{T},$ (2)
$\displaystyle\int_{0}^{T_{k}}r_{k,j}(t)\mathop{}\\!\mathrm{d}{t}\geq
v_{k,j}z_{k},\quad\forall j\in\mathcal{F}_{k},\forall k\in\mathcal{C},$ (3)

Constraint (2) expresses that, at any instant $t$ within the time horizon
${\cal T}$, the total rate that port $\ell$ assigns to flows cannot exceed its
capacity $B_{\ell}$. Constraint (3) ensures that the data of flows of each
accepted coflow $k$ should be completely transmitted before the deadline
$T_{k}$.

######  Lemma 1 (Proposition 1 in [8]).

There exists a polynomial time reduction of the CAR problem (P1) to the
problem of minimizing the number of late jobs in a concurrent open shop [10].
Hence, the CAR problem is NP-hard.

TABLE I: Main notations.  Symbol | Description
---|---
${\cal L}$ | set of fabric ports
$M$ |  number of machines, $M=|{\cal L}|/2$
$B_{\ell}$ |  available bandwidth of port $\ell\in\mathcal{L}$
${\cal C}$ |  set of coflows. ${\cal C}$ has cardinality of $N$
$\sigma$ |  scheduling order of coflows, $\sigma=\left\\{\sigma_{1},\cdots,\sigma_{N-1},\sigma_{N}\right\\}$
$T_{k}$ |  deadline of coflow $k$
${\cal F}_{k}$ (${\cal F}_{\ell,k}$) |  set of flows of coflow $k$ (that use port $\ell$)
$v_{k,j}$ |  volume of flow $j\in{\cal F}_{k}$ of coflow $k$
$\hat{v}_{\ell,k}$ |  total volume transmitted by coflow $k$ on port $\ell$
$p_{\ell,k}$ |  processing time of coflow $k$ on port $\ell$
$c_{k}$ ($c_{\ell,k}$) |  completion time of coflow $k$ on port $\ell$
${CT}_{k,j}$ |  completion time of flow $j$ of coflow $k$

###  II-A Motivating Example

We will take $\mathtt{CS\text{-}MHA}$ [11] as our starting point to address
the CAR problem. $\mathtt{CS\text{-}MHA}$ has introduced a new direction to
solve the scheduling problem by means of a static coflow prioritisation. The
prioritisation is used in order to approximate the solution of the coflow
scheduling problem that maximizes the CAR. $\mathtt{CS\text{-}MHA}$ computes
the scheduling orders at each port using the Moore-Hodgon’s algorithm which
also determines the set of admitted coflows at each port
$\ell\in{\mathcal{L}}$. Since different ports could have different sets of
admitted coflows, a coflow is admitted if it is admitted at all ports. For all
rejected coflows, a second round is applied to check if some rejected coflows
can actually satisfy their deadline. The algorithm selects a coflow with the
minimum bandwidth required at the bottleneck port since it is more likely to
catch up with its deadline.

The simple example depicted in Fig. 1 illustrates the limitations of
$\mathtt{CS\text{-}MHA}$. It uses the standard _Big-Switch_ model to abstract
a datacenter fabric: the example will be used as a running example throughout
the paper. The instance contains $5$ coflows: $C_{1}$ has $4$ flows and
$C_{2}$, $C_{3}$, $C_{4}$ and $C_{5}$ each have one flow. To ease the
presentation, the flows are organised in virtual output queues at the ingress
ports. The virtual queue index represents the flow output port, modulo the
number of machines. The numbers on the flows’ representations correspond to
their normalized volumes. All fabric ports have the same normalized bandwidth
of $1$.

Figure 1: Motivating example with a Big-Switch fabric composed of $4$
ingress/egress ports connecting to $4$ machines. Flows in ingress ports are
organized by destinations and are color-coded by coflows. The example has $5$
coflows. Coflow $C_{1}$ (blue) has $4$ flows, with each ingress port sending
$1$ units of data to one egress port: its deadline is $1$; coflows $C_{2}$
(green), $C_{3}$ (red), $C_{4}$ (orange) and $C_{5}$ (purple) have a single
flow, each sending $(1+\varepsilon)$ unit of data. The deadline of these
coflows is $2$.

At the first iteration, $\mathtt{CS\text{-}MHA}$ computes the scheduling order
on each port using Moore-Hodgson algorithm [12], which is based on the
Earliest Due Date rule and is known to minimize the number of missed deadlines
on a single machine (or port in the coflow context). Since coflow $C_{1}$ uses
all ports and has the smallest deadline ($T_{1}=1$), then at each port,
$C_{1}$ will be scheduled first. As consequence, all other coflows are
rejected since they cannot satisfy their deadline when scheduled after
$C_{1}$. Given that coflow scheduling, the CAR is $\frac{1}{5}$. However, an
optimal scheduling solution would be $\\{C_{2},C_{3},C_{4},C_{5},C_{1}\\}$ or
any combination that has coflow $C_{1}$ as the last one to be scheduled. The
latter scheduler attains a CAR of $\frac{4}{5}$. From this example, we can see
that $\mathtt{CS\text{-}MHA}$ can be made arbitrarily worse compared to the
optimal solution.

The above example can easily be extended to show that in the worst-case
scenario, the CAR of $\mathtt{CS\text{-}MHA}$ can be made arbitrarily small
(as close to $0$ as needed). Consider $M$ machines, coflow $C_{1}$ that uses
all ports, and coflows $C_{2},\ldots,C_{M}$ with one flow each. The other
parameters remain the same as above. It can be seen that the CAR obtained
using $\mathtt{CS\text{-}MHA}$ and $\mathtt{DCoflow}$ are $\frac{1}{M}$,
$\frac{M-1}{M}$, respectively. In this setting, the CAR obtained using
$\mathtt{CS\text{-}MHA}$ is close to zero when $M$ is high, while using
$\mathtt{DCoflow}$, the CAR is close to one.

The key observation is that $\mathtt{CS\text{-}MHA}$ neglects the impact that
a coflow may have on other coflows on multiple ports. Indeed, a coflow causing
multiple deadlines to be missed should have lesser priority, even when its
deadline is the earliest. When this is neglected, the coflow ordering is
misjudged thus degrading the CAR. In this work, we start from this observation
and propose a new $\sigma$-order scheduler, called $\mathtt{DCoflow}$.

###  II-B Coflow Ordering

Problem (P1) implicitly depends on the CCT. An alternative approach to
maximize the CAR is to order the coflows in some appropriate way, and then to
leverage the priority forwarding mechanisms of the underlying transport
network [13].

In this approach, once such an ordering $\sigma$ is determined, it is enough
to adopt work-conserving transmission policies and use forwarding priorities
in such a way that a flow $j\in{\cal F}_{\sigma_{k}}$ is blocked if and only
if either its ingress or egress port is busy serving a flow
$j^{\prime}\in{\cal F}_{\sigma_{k^{\prime}}}$ for some $k^{\prime}<k$, i.e., a
flow of a coflow with higher priority according to $\sigma$. Such a flow
scheduling is called _$\sigma$ -order-preserving_.

The formulation of Problem (P1) can be transformed into an Integer Linear
Program (ILP) as follows. We define the binary variable $\delta_{i,k}$ as $1$
if coflow $i$ has higher priority than coflow $k\neq i$, and $0$ otherwise.
The ordering of coflows can then be modeled using the following standard
disjunctive constraints

$\displaystyle\delta_{k,k^{\prime}}+\delta_{k^{\prime},k}$ $\displaystyle=$
$\displaystyle 1,\quad\forall k,k^{\prime}\in{\cal C},$ (4)
$\displaystyle\delta_{k,k^{\prime}}+\delta_{k^{\prime},k"}+\delta_{k",k}$
$\displaystyle\leq$ $\displaystyle 2,\quad\forall k,k^{\prime},k"\in{\cal C}.$
(5)

It should be clear that the ordering $\sigma$ can easily be derived from the
variables $\\{\delta_{k,k^{\prime}}\\}_{k,k^{\prime}\in{\cal C}}$. This
ordering should be such that as many coflows as possible are accepted, that
is, $\sum_{i\in{\cal C}}z_{k}$ is maximized. A central difficulty here is to
compute the completion time of a coflow, which stems from the fact that data
transmissions on the various ports are not independent: the transmission of a
flow may be blocked until the ingress port becomes available even if the
egress port is idle, and vice versa. However, assuming that the ports are
independent, we can obtain a lower bound on the completion time of coflow $k$
on port $\ell$ as follows

$c_{\ell,k}\geq\sum_{k^{\prime}\neq
k}p_{\ell,k^{\prime}}\delta_{k^{\prime},k}z_{k^{\prime}}+p_{\ell,k},\quad\forall\ell\in{\cal
L},k\in{\cal C}.$ (6)

This lower bound assumes that the transmission of coflow $k$ on port $\ell$
can start as soon as all flows of all coflows $k^{\prime}$ scheduled before
$k$ have been transmitted on port $\ell$. Note that the above constraint is
not linear due to the product $\delta_{k^{\prime},k}z_{k^{\prime}}$. However,
it can easily be linearized by introducing binary variables $y_{k^{\prime},k}$
satisfying the constraints

$y_{k^{\prime},k}\leq z_{k^{\prime}};\quad
y_{k^{\prime},k}\leq\delta_{k^{\prime},k};\quad y_{k^{\prime},k}\geq
z_{k^{\prime}}+\delta_{k^{\prime},k}-1,$ (7)

(6) thus becomes

$c_{\ell,k}\geq\sum_{k^{\prime}\neq
k}p_{\ell,k^{\prime}}y_{k^{\prime},k}+p_{\ell,k},\quad\forall\ell\in{\cal
L},k\in{\cal C}.$ (8)

A lower bound on the completion time of coflow $k$ is then obtained as
$c_{k}=\max_{\ell\in{\cal L}}c_{\ell,k}$, and this coflow can meet its
deadline only if $c_{k}\leq T_{k}$. The deadline constraint can then be
described as

$c_{\ell,k}\leq T_{k}z_{k},\quad\forall\ell\in{\cal L},k\in{\cal C}.$ (9)

We can thus obtain an upper bound on the number of accepted coflows by solving
the following ILP,

$\mathrm{max}\enskip\sum_{k\in\mathcal{C}}z_{k},\quad\mathrm{s.t.}\enskip(\ref{eq:order-
disjunct-1},\ref{eq:order-disjunct-2},\ref{eq:y1},\ref{eq:lower-bound-CCT-
new},\ref{eq:deadline-constraint}).$

It is easy to verify that Problem (II-B) is NP-hard and so that it may become
unfeasible to find an optimal solution in large-scale datacenter networks. In
the next section, we hence propose an efficient heuristic algorithm to
determine the order in which coflows must finish in order to maximize the
coflow acceptance rate.

##  III $\sigma$-Order Scheduling with $\mathtt{DCoflow}$

In this section, we present $\mathtt{DCoflow}$, an algorithm to solve the
problem of joint coflow admission control and scheduling. Given a list of $N$
coflows, it provides a permutation
$\sigma=(\sigma_{1},\sigma_{2},..,\sigma_{N})$ of these coflows, with the aim
of maximizing the coflow acceptance rate. We observe that, when a coflow is
served according to the corresponding $\sigma$-order-preserving schedule,
lower-priority coflows will be impacted. It is hence possible that some of
these lower-priority coflows will not respect their deadline. Hence our
objective is to evaluate the impact of a scheduled coflow with respect to the
original CAR problem. The problem is combinatorial in nature, since there are
$N!$ possible coflow orderings. Hence, we propose a new approach that uses the
formulation of Problem (II-B) by deriving a necessary condition satisfied by
all feasible coflow orders solving Problem (II-B).

Consider a feasible solution to Problem (II-B) and a subset ${\cal
S}\subseteq{\cal C}$ of coflows. Given a coflow $k\in{\cal C}$, let ${\cal
S}_{k}^{-}=\left\\{k^{\prime}\in{\cal S}\ :\ \delta_{k^{\prime},k}=1\right\\}$
be the set of coflows in ${\cal S}$ which are scheduled before $k$. Condition
(6) may be rewritten as

$c_{\ell,k}\geq p_{\ell,k}+\sum_{k^{\prime}\in{\cal
S}_{k}^{-}}p_{\ell,k^{\prime}}z_{k^{\prime}},$ (10)

which implies

$c_{\ell,k}p_{\ell,k}z_{k}\geq
p_{\ell,k}^{2}z_{k}^{2}+p_{\ell,k}z_{k}\sum_{k^{\prime}\in{\cal
S}_{k}^{-}}p_{\ell,k^{\prime}}z_{k^{\prime}}.$ (11)

Using the inequality $T_{k}\geq c_{\ell,k}z_{k}$ and summing over all coflows
$k\in{\cal S}$, we obtain

$\displaystyle\sum_{k\in\mathcal{S}}p_{\ell,k}T_{k}\geq$
$\displaystyle\sum_{k\in\mathcal{S}}c_{\ell,k}p_{\ell,k}z_{k}\geq\sum_{k\in\mathcal{S}}\left(p_{\ell,k}z_{k}\right)^{2}+\hskip-11.38109pt\sum_{k\in\mathcal{S},k^{\prime}\in\mathcal{S}_{k}^{-}}\hskip-11.38109ptp_{\ell,k}p_{\ell,k^{\prime}}z_{k}z_{k^{\prime}}$
$\displaystyle\hskip-54.06023pt=\frac{1}{2}\sum_{k\in\mathcal{S}}\left(p_{\ell,k}z_{k}\right)^{2}+\frac{1}{2}\left[\sum_{k\in\mathcal{S}}\left(p_{\ell,k}z_{k}\right)^{2}+2\sum_{k\in\mathcal{S}}p_{\ell,k}z_{k}\hskip-5.69054pt\sum_{k^{\prime}\in\mathcal{S}_{k}^{-}}p_{\ell,k^{\prime}}z_{k^{\prime}}\\!\right]$
$\displaystyle\hskip-54.06023pt=\frac{1}{2}\sum_{k\in\mathcal{S}}\left(p_{\ell,k}z_{k}\right)^{2}+\frac{1}{2}\left(\sum_{k\in\mathcal{S}}p_{\ell,k}z_{k}\right)^{2}.$
(12)

We thus conclude that any feasible solution to Problem (II-B) satisfies the
condition $\sum_{k\in{\cal S}}p_{\ell,k}T_{k}\geq f_{\ell}({\cal S})$ for any
subset ${\cal S}\subseteq{\cal C}$ of _accepted_ coflows, where
$f_{\ell}({\cal S})=\frac{1}{2}\sum_{k\in{\cal
S}}p_{\ell,k}^{2}+\frac{1}{2}\left(\sum_{k\in{\cal S}}p_{\ell,k}\right)^{2}$.
These conditions are the so-called parallel inequalities and provide valid
inequalities for the concurrent open shop problem [14]. Note that they do not
depend on the ordering of the coflows.

We shall use the parallel inequalities to determine the coflows that should
not be admitted. Consider a solution to Problem (II-B) and assume that, in
this solution, there exists a subset ${\cal S}$ of accepted coflows (i.e.,
$z_{k}=1$ for all $k\in{\cal S}$) such that $\sum_{k\in{\cal
S}}p_{\ell,k}T_{k}<f_{\ell}({\cal S})$ for at least one port $\ell\in{\cal
L}$. This implies that this solution is not feasible and can only become
feasible by rejecting a coflow $k^{\prime}$ among the accepted coflows using
port $\ell$. We choose this coflow $k^{\prime}$ so as to minimize the quantity
$f_{\ell}({\cal S}\setminus\\{k^{\prime}\\})-\sum_{k\in{\cal
S}\setminus\\{k^{\prime}\\}}p_{\ell,k}T_{k}$, in the hope that it becomes
negative. Observe that

$\displaystyle f_{\ell}\left({\cal S}\right)$
$\displaystyle\\!=\\!\left[\frac{1}{2}p_{\ell,k^{\prime}}^{2}\\!+\\!\frac{1}{2}\\!\sum_{k\in{\cal
S}\setminus\\{k^{\prime}\\}}p_{\ell,k}^{2}\right]\\!+\\!\frac{1}{2}\\!\left(p_{\ell
k^{\prime}}+\\!\sum_{k\in{\cal S}\setminus\\{k^{\prime}\\}}p_{\ell
k}\right)^{2}$ $\displaystyle=f_{\ell}\left({\cal
S}\setminus\\{k^{\prime}\\}\right)+p_{\ell,k^{\prime}}\sum_{k\in{\cal
S}}p_{\ell,k},$ (13)

from which it follows that

$\displaystyle f_{\ell}\left({\cal
S}\setminus\\{k^{\prime}\\}\right)-\sum_{k\in{\cal
S}\setminus\\{k^{\prime}\\}}p_{\ell,k}T_{k}$
$\displaystyle=f_{\ell}\left({\cal S}\right)-\sum_{k\in{\cal
S}}p_{\ell,k}T_{k}$
$\displaystyle+p_{\ell,k^{\prime}}\left(T_{k^{\prime}}-\sum_{k\in{\cal
S}}p_{\ell,k}\right)$ (14)

The above relation will be used to define our coflow admission control
algorithm. In general, it suggests possible heuristics to remove a coflow
$k^{\prime}$ from ${\cal S}$ in order to satisfy the parallel inequalities.
For instance, one such heuristics is to minimize the quantity

$\Psi_{\ell,k^{\prime}}:=p_{\ell,k^{\prime}}\left(T_{k^{\prime}}-\sum_{k\in{\cal
S}}p_{\ell,k}\right).$ (15)

In fact, the term $\sum_{k\in{\cal S}}p_{\ell,k}-T_{k^{\prime}}$ represents
how much the completion time of coflow $k^{\prime}$ on port $\ell$ exceeds its
deadline, assuming that it is scheduled latest. Hence, we could reject coflow
$k^{\prime}$ with both large processing time $p_{\ell,k^{\prime}}$ on port
$\ell$ and large deadline violation. More admission rationales based on (14)
will be detailed in describing $\mathtt{DCoflow}$, whose pseudocode is
reported in Algorithm 1. Generally, $\mathtt{DCoflow}$ takes a list of
unsorted coflows and provides as output the scheduling order of accepted
coflows. It works in rounds and at each round, it either accepts a coflow or
it rejects one.

1

2${\cal S}=\left\\{1,2,\ldots,N\right\\}$; $\triangleright$ initial set of
unscheduled coflows

3$\sigma=\varnothing;$ $\triangleright$ initial scheduling order

4$\sigma^{\star}=\varnothing;$ $\triangleright$ initial set of pre-rejected
coflows

5$n=N$; $\triangleright$ round counter

6while  _ ${\cal S}\neq\varnothing$ _ do

7

8 $\ell_{b}=\underset{\ell\in{{\cal L}}}{\mathrm{arg\,max}}\,\sum_{k\in{\cal
S}}p_{\ell,k}$; $\triangleright$ bottleneck port

9 ${\cal S}_{\ell_{b}}=\left\\{k\in{{\cal S}}:p_{\ell_{b},k}>0\right\\}$;
$\triangleright$ set of coflows using $\ell_{b}$

10 # Coflows that can finish in time when scheduled last

11 ${\cal S}^{d}_{\ell_{b}}\leftarrow\\{j\in{\cal
S}_{\ell_{b}}|\sum_{k\in{\cal S}_{\ell_{b}}}p_{\ell_{b},k}\leq T_{j}\\}$;

12 if  _$\displaystyle{\cal S}^{d}_{\ell_{b}}\not=\varnothing$_ then

13

14 $k_{n}=\underset{k\in{{\cal C}}}{\mathrm{arg\,max}}\,T_{k}$;
$\triangleright$ admit coflow with largest deadline

15 $\sigma_{n}=k_{n}$; $\triangleright$ append $k_{n}$ to $\sigma$

16 ${\cal S}={\cal S}\setminus\left\\{k_{n}\right\\}$; $\triangleright$ remove
$k_{n}$ from ${\cal S}$

17 else

18

19 $k^{\star}=\mathtt{RejectedCoflow}({\cal S})$; $\triangleright$ select a
coflow to reject

20 $\sigma_{n}=k^{\star}$; $\triangleright$ append $k^{\star}$ to $\sigma$;

21 $\sigma^{\star}=\sigma^{\star}\cup\left\\{k^{\star}\right\\}$;
$\triangleright$ append $k^{\star}$ to $\sigma^{\star}$

22 ${\cal S}={\cal S}\backslash\left\\{k^{\star}\right\\}$; $\triangleright$
remove $k^{\star}$ from ${\cal S}$

23

24 $n=n-1$; $\triangleright$ update the round index

25$\sigma=\mathtt{RemoveLateCoflows}\left(\sigma,\sigma^{\star}\right)$;

26return $\sigma$; $\triangleright$ final scheduling order

27

28

29

30Function  _RemoveLateCoflows(_$\sigma$ , $\sigma^{\star}$_)_:

31

32 while  _ $\sigma^{\star}\neq\varnothing$ _ do

33

34
$k^{\star}\leftarrow\underset{k}{\mathrm{arg\,min}}\left\\{\sigma_{k}\in\sigma^{\star}\right\\}$;
$\triangleright$ first coflow in $\sigma^{\star}$

35
$c_{\sigma_{k^{\star}}}\leftarrow\mathtt{evalCCT}\left(\left\\{\sigma_{k}\right\\}_{k\in\left[1,k^{\star}\right]}\right)$;
$\triangleright$ CCT of coflow $\sigma_{k^{\star}}$

36 if  _ $c_{\sigma_{k^{\star}}}>T_{\sigma_{k^{\star}}}$_ then

37 $\sigma\leftarrow\sigma\backslash\left\\{\sigma_{k^{\star}}\right\\}$;

38
$\sigma^{\star}\leftarrow\sigma^{\star}\backslash\left\\{\sigma_{k^{\star}}\right\\}$;

39 return $\sigma$;

Algorithm 1 $\mathtt{DCoflow}$

$\mathtt{DCoflow}$ starts by computing the total completion time at each port
and finds the bottleneck $\ell_{b}$, i.e., the port with the largest
completion time. Having this, it determines $\mathcal{S}^{d}_{\ell_{b}}$, the
set of coflows active on the bottleneck port $\ell_{b}$ (line $8$) which
complete before their deadline if scheduled last on $\ell_{b}$. If
$\mathcal{S}^{d}_{\ell_{b}}$ is not empty (Lines 10–13), it selects a coflow
$k_{n}\in\mathcal{S}^{d}_{\ell_{b}}$ that has the largest deadline in
$\mathcal{S}^{d}_{\ell_{b}}$. If $\mathcal{S}^{d}_{\ell_{b}}$ is empty (Lines
14–18), i.e., no coflow respects its deadline when scheduled last on the
bottleneck (Line $13$), the algorithm selects a coflow to be removed using
function $\mathtt{RejectedCoflow}$, whose aim is to comply with parallel
inequalities according to (14). We consider two variants, each of which
corresponds to a slightly different criterion to select the candidate coflow
$k^{\star}$ by taking into account its weight.

The first variant, namely $\mathtt{DCoflow\\_v1}$, finds the candidate
$k^{\star}\in\mathcal{S}_{b}$ and every port $\ell$ used by $k$ where it fails
to meet the deadline (that is, $\Psi_{\ell,k}<0$),

$k^{\star}=\underset{k\in\mathcal{S}_{b}}{\mathrm{arg\,min}}\,{\displaystyle\Big{(}\sum_{\ell:\Psi_{\ell,k}<0}}\Psi_{\ell,k}\Big{)}.$
(16)

The second variant, namely $\mathtt{DCoflow\\_v2}$, finds the candidate
$k^{\star}\in\mathcal{S}_{b}$ and every port $\ell$ used by $k$ that has at
least $\gamma$ times the congested level of the bottleneck, i.e.,
$\ell:\sum_{j}p_{\ell j}\geq\gamma\sum_{j}p_{bj}$,

$k^{\star}=\underset{k\in\mathcal{S}_{b}}{\mathrm{arg\,min}}\,{\displaystyle\Big{(}\sum_{\ell:\sum_{j}p_{\ell
j}\geq\gamma\sum_{j}p_{bj}}}\Psi_{\ell k}\Big{)}.$ (17)

Once the initial scheduling order $\sigma$ is obtained, $\mathtt{DCoflow}$
uses the function $\mathtt{RemoveLateCoflows}$ to estimate the CCT and remove
from $\sigma$ the coflows that do not satisfy the deadline constraint (Line
$20$). Briefly, $\mathtt{RemoveLateCoflows}$ considers each pre-rejected
coflow $k^{\star}$ in $\sigma^{\star}$ (Line $24$). It calls the function
$\mathtt{evalCCT}$ to evaluate the CCT of $k^{\star}$, given all the coflows
in $\sigma$ that are scheduled before $k^{\star}$ (Line $25$). If $k^{\star}$
cannot meet its deadline, i.e.,
$c_{\sigma_{k^{\star}}}>T_{\sigma_{k^{\star}}}$, it will be removed
permanently from $\sigma^{\star}$ and $\sigma$. $\mathtt{RemoveLateCoflows}$
performs a swipe on the admitted coflows (as in [11]) and iteratively removes
coflows in $\sigma$ that belong to $\sigma^{\star}$ until the estimated CCT
$c_{\sigma_{k}}$ of each coflow $\sigma_{k}\in\sigma$ is of at most
$T_{\sigma_{k}}$.

It is important to note that the final solution yielded by $\mathtt{DCoflow}$
does not guarantee all coflows in $\sigma$ to satisfy their deadlines. In Sec.
IV, the prediction error of $\mathtt{DCoflow}$ indicates the gap between the
estimated CAR and the actual CAR after applying resource allocation.

Example. To illustrate the difference between $\mathtt{DCoflow}$ and
$\mathtt{CS\text{-}MHA}$, we consider again the example illustrated in Fig. 1.
Table II shows the execution of $\mathtt{DCoflow\\_v1}$ on that example. At
the first step, $\mathtt{DCoflow\\_v1}$ chooses bottleneck ingress port $1$,
which is used by coflows $C_{1}$ and $C_{2}$. It computes
$\overline{\Psi}_{k}=\sum_{\ell:\Psi_{\ell,k}<0}\Psi_{\ell,k}$ for both
coflows and selects the coflow with smallest $\bar{\Psi}_{k}$ (in this case,
$C_{1}$) to be scheduled last. Since all unscheduled coflows do not share any
port in the fabric, any ordering of remaining coflows gives the same average
CAR. Given the final schedule, $\mathtt{DCoflow}$ obtains $\frac{4}{5}$ as the
average CAR, which is the optimal solution.

TABLE II: Execution of $\mathtt{DCoflow}$ on the example of Fig. 1.  Unscheduled coflows (set $\mathcal{S}$) |  $\ell_{b}$ |  $\left\\{\overline{\Psi}_{1},\overline{\Psi}_{2},\overline{\Psi}_{3},\overline{\Psi}_{4},\overline{\Psi}_{5}\right\\}$
---|---|---
$\mathcal{S}=\left\\{\boldsymbol{C_{1}},\boldsymbol{C_{2}},C_{3},C_{4},C_{5}\right\\}$ |  1 |  $\left\\{-4\left(1+\varepsilon\right),-\varepsilon,\,\cdot\,,\,\cdot\,,\,\cdot\,\right\\}$
$\mathcal{S}=\left\\{\boldsymbol{C_{2}},C_{3},C_{4},C_{5}\right\\}$ |  1 |  $\left\\{\,\cdot\,,0,\,\cdot\,,\,\cdot\,,\,\cdot\,\right\\}$
$\mathcal{S}=\left\\{\boldsymbol{C_{3}},C_{4},C_{5}\right\\}$ |  2 |  $\left\\{\,\cdot\,,\,\cdot\,,0,\,\cdot\,,\,\cdot\,\right\\}$
$\mathcal{S}=\left\\{\boldsymbol{C_{4}},C_{5}\right\\}$ |  3 |  $\left\\{\,\cdot\,,\,\cdot\,,\,\cdot\,,0,\,\cdot\,\right\\}$
$\mathcal{S}=\left\\{\boldsymbol{C_{5}}\right\\}$ |  4 |  $\left\\{\,\cdot\,,\,\cdot\,,\,\cdot\,,\,\cdot\,,0\right\\}$

Online Implementation. $\mathtt{DCoflow}$ can be also be run online, when
coflows arrive sequentially and possibly in batches. For this, define $f$ to
be the frequency of updates, i.e., instants at which $\mathtt{DCoflow}$
recomputes a schedule. The updates can be performed either at arrival instants
of coflows (in which case we set $f=\infty$) or periodically with period
$1/f$. We assume that the scheduler knows the volumes of the flows of each
arrived coflow. However, it neither knows the volumes nor the release times of
future coflows.

At each update instant, $\mathtt{DCoflow}$ recomputes the scheduling order for
coflows currently available in the network. These include the ones that were
scheduled in the previous scheduling instants and have not yet finished; the
ones that were rejected in the previous scheduling instants but whose deadline
has not yet expired and the arrivals during the update interval.
$\mathtt{DCoflow}$ calculates the new order for this set of coflows based on
the remaining volumes of the flows and not on the original volumes. Note that
it is assumed that coflows can be preempted [4]. This process is repeated at
each update instant.

Complexity Analysis. The complexity of $\mathtt{DCoflow}$ is ${\cal
O}(N^{2})$. Specifically, in $\mathtt{DCoflow}$, the values $\sum_{i\in{\cal
S}}p_{\ell,i}$ and $\Psi_{\ell,x}$ of each remaining coflow $x$ can be pre-
stored by calculating them at a cost ${\cal O}(NL)$, where $L=|{\cal L}|$ is
the number of ports. Then these values can be updated at a cost ${\cal O}(L)$
per coflow at each iteration. The number of operations required at Lines
$11$–$13$ is ${\cal O}(N)$ and at Lines $15$–$16$ is ${\cal O}(N)$, so that
finally across iterations it adds to ${\cal O}(N^{2})$. On the other hand, it
is easy to verify the complexity of $\mathtt{RemoveLateCoflows}$ is ${\cal
O}(NL)$.

##  IV Performance Evaluation

We evaluate via simulations our proposed heuristics (two variants of
$\mathtt{DCoflow}$111 The flow-level simulator and the implementation of all
algorithms tested in this paper are available at
https://github.com/luuquangtrung/CoflowSimulator.) along with some existing
algorithms such as $\mathtt{CS\text{-}MHA}$222 We only reimplemented the
centralized algorithm ($\mathtt{CS\text{-}MHA}$) presented in [11], which has
been reported to be better than the decentralized version
($\mathtt{D\textsuperscript{2}\text{-}CAS}$) in terms of CAR.  and the
solution provided by the optimization method $\mathtt{CDS\text{-}LP}$ proposed
in [8]. The relaxed version of $\mathtt{CDS\text{-}LP}$, named
$\mathtt{CDS\text{-}LPA}$, is also implemented333 It is worth noting that the
formulation of $\mathtt{CDS\text{-}LP}$ and $\mathtt{CDS\text{-}LPA}$ use the
same decision variables $\\{z_{k}\\}_{k\in{\cal C}}$ as those introduced in
Problem (P1). In $\mathtt{CDS\text{-}LP}$, $z_{k}$ are binaries, whereas in
$\mathtt{CDS\text{-}LPA}$, $z_{k}$ are continuous variables taking values in
$[0,1]$. For any solution yielded by $\mathtt{CDS\text{-}LPA}$, only coflows
$k$ whose ${z_{k}}$ strictly equals $1$ are considered as accepted ones. . By
using the solution derived from $\mathtt{CDS\text{-}LP}$ as an upper bound, we
would get the sense of how close the algorithms are to the optimum. A brief
description of the reference algorithms has been given in Sec. I. We also
compare the performance of our schedulers against $\mathtt{Sincronia}$ [13]
and $\mathtt{Varys}$ [9] that aim to minimize the average CCT.

Once we obtained the $\sigma$-order, the actual coflow resource allocation for
our solution is implemented by the greedy rate allocation algorithm
$\mathtt{GreedyFlowScheduling}$ [13]. At any given point of time,
$\mathtt{GreedyFlowScheduling}$ reserves the full port bandwidth to one flow
at the time. It does so by complying to the scheduling order in $\sigma$ of
the coflow to which the flow belongs [13]. We note that, in the case of
$\mathtt{CDS\text{-}LP}$, $\mathtt{CDS\text{-}LPA}$, and $\mathtt{Varys}$,
instead, the rate allocation is part of the algorithmic solution.

The network fabric is represented by $M$ machines or end-hosts connected to a
non-blocking Big-Switch fabric, of which each access port has a normalized
capacity of $1$. The algorithms will be evaluated on both small-scale and
large-scale networks, where a network is denoted by $[M,N]$ to indicate
different fabric size and number of coflows ($N$) used in the simulations.
Small-scale networks have a fabric of size $M=10$, whereas large-scale
networks have a fabric with either $50$ or $100$ machines. The coflows in
these networks are generated from both synthetic and real traffic traces.

The MILP solver $\mathtt{gurobi}$ is used to solved $\mathtt{CDS\text{-}LP}$
and $\mathtt{CDS\text{-}LPA}$. Due to the high complexity,
$\mathtt{CDS\text{-}LP}$ and $\mathtt{CDS\text{-}LPA}$ are only evaluated on
small-scale networks. In what follows, the detailed setup, comparison metrics,
and simulation results are presented.

###  IV-A Simulation Setup

Synthetic Traffic. The synthetic traffic comprises two types of coflows.
Type-$1$ coflows have only one flow, whereas the number of flows of Type-$2$
coflows follows a uniform distribution in $[2M/3,M]$. Each generated coflow is
randomly assigned to either Class $1$ or Class $2$ with probability of
respectively $0.6$ and $0.4$. Moreover, each coflow $k$ is assigned a random
deadline within $[\text{CCT}_{k}^{0},2\text{CCT}_{k}^{0}]$), where
$\text{CCT}_{k}^{0}$ is the completion time of coflow $k$ _in isolation_.
Flows of Class-$1$ coflows are assigned a random volume of mean of $1$ and
standard deviation of $0.2$. The volume ratio for the flows of Class-$1$ and
Class-$2$ coflows is $0.8$.

Real Traffic. Real traffic datasets are obtained by the Facebook traces
dataset [9], based on a MapReduce shuffle trace collected from one of
Facebook’s $3000$-machine cluster with $150$ racks. The data traces consist of
$526$ coflows. It has a skewed coflow width distribution, ranging from coflows
with a single flow to very large ones (the largest coflow has $21170$ flows).
For detailed statistics of the Facebook traces we refer the reader to [15].
For each configuration $[M,N]$, $N$ coflows are randomly sampled from the
Facebook dataset. Coflows are only selected from the ones that have at most
$M$ flows. The volume of each flow is given by the dataset.

Metric. We evaluate the algorithms based on the average CAR. We also present
the gains in percentiles of each algorithm with respect to the solution
yielded by $\mathtt{CDS\text{-}LP}$ in terms of CAR. These gains are
calculated using the formula: $\text{average gain in CAR}=\frac{\text{compared
CAR}}{\text{CAR under }\mathtt{CDS\text{-}LP}}-1.$

###  IV-B Results with Offline Setting

In the offline setting, we consider that all coflows arrive at the same time,
i.e., their release time is zero. For each simulation with a specific scale of
the network and either synthetic or real traffic traces, we randomly generate
$100$ different instances and compute the average performance of algorithms
over $100$ runs.

####  IV-B1 Average CAR Under Synthetic Traffic

Figs. 2a–2b show the average CAR with respectively small-scale networks and
large-scale networks. The percentile gains of each algorithm with respect to
$\mathtt{CDS\text{-}LP}$ are shown in Fig. 4a, in terms of average CAR for the
configuration $[10,60]$.

(a) Synthetic traffic traces on a small-scale network.  (b) Synthetic traffic
traces on a large-scale network.

Figure 2: Average CAR with synthetic traffic traces using (a) small-scale and
(b) large-scale networks. Each point in the x-axis represents the network
$[M,N]$.

It is observed that our proposed heuristics are closest in terms of CAR to the
optimum ($\mathtt{CDS\text{-}LP}$) than all other algorithms, with both small-
and large-scale networks. Among two variants of $\mathtt{DCoflow}$,
$\mathtt{DCoflow\\_v1}$ yields the best performance and even outperforms
$\mathtt{CDS\text{-}LPA}$, the approximation version of
$\mathtt{CDS\text{-}LP}$. For instance, on the network $[10,10]$,
$\mathtt{DCoflow\\_v1}$ improves the CAR on average by $6.5\%$, $11.5\%$,
$15.1\%$, and $26.6\%$, compared respectively to $\mathtt{CDS\text{-}LPA}$,
$\mathtt{CS\text{-}MHA}$, $\mathtt{Sincronia}$, and $\mathtt{Varys}$.
Interestingly, the improvement becomes higher when the load is increased. For
example, the corresponding improvement on average CAR on a $[10,60]$ network
are $67.2\%$, $98.3\%$, $59.9\%$, and $36.8\%$ (see Fig. 2a). The improvement
is even higher when performed on a large-scale netork. For example, compared
to $\mathtt{CS\text{-}MHA}$, $\mathtt{Sincronia}$, and $\mathtt{Varys}$, on
the network $[100,400]$, the improvement in terms of average CAR are
respectively $648.1\%$, $32.3\%$, and $17.9\%$. (see Fig. 2b). It is worth
noting how the performance of $\mathtt{CS\text{-}MHA}$ falls drastically when
dealing with large-scale networks. This is expected since
$\mathtt{CS\text{-}MHA}$ computes prioritizes coflows that use a large number
of ports over those that use a few. In instances with a large number of
coflows of the latter type, the CAR of $\mathtt{CS\text{-}MHA}$ goes to $0$
(see detailed explanation with the motivating example in Sec. II-A)

The result in Fig. 4a shows that the two variants of $\mathtt{DCoflow}$
achieve a smaller gap to the optimal in almost all values of percentile
compared to other algorithms. For instance, compared to $\mathtt{Sincronia}$,
$\mathtt{DCoflow\\_v1}$ improves the CAR in $50\%$ of $100$ instances by
$50\%$ and it achieves around $43\%$ at $99$th percentile.

####  IV-B2 Average CAR Under Real Traffic Traces

This section presents the results obtained with the Facebook traffic traces,
using the same configurations as those used in Sec. IV-B1. Figs. 3a–3b show
the average CAR with respectively small- and large-scale networks. The gains
in percentiles of each algorithm with respect to $\mathtt{CDS\text{-}LP}$, in
terms of average CAR when using a $[10,60]$ fabric are shown in Fig. 4b.
Similar to what observed in the results with the synthetic traces (see Sec.
IV-B1), $\mathtt{DCoflow\\_v1}$ and $\mathtt{DCoflow\\_v2}$ yield a
significant improvement in terms of average CAR compared to other heuristics.
For instance, with a $[10,60]$ configuration, the two variants of
$\mathtt{DCoflow}$ improve the average CAR on average by $24.4\%$, $25\%$,
$52.2\%$, $93.1\%$, compared respectively to $\mathtt{CDS\text{-}LPA}$,
$\mathtt{CS\text{-}MHA}$, $\mathtt{Sincronia}$, and $\mathtt{Varys}$. (see
Fig. 3a). The improvement is even higher when performed on a large-scale
network. For example, compared to $\mathtt{CS\text{-}MHA}$,
$\mathtt{Sincronia}$, and $\mathtt{Varys}$, on $[100,400]$ network, the
improvement in terms of average CAR are respectively $36.6\%$, $55.3\%$, and
$147.5\%$.

Moreover, the results in Fig. 4b show that the two variants of
$\mathtt{DCoflow}$ achieve a smaller gap to the optimal in almost all values
of percentile compared to the other algorithms. For instance, compared to
$\mathtt{Sincronia}$, $\mathtt{DCoflow\\_v1}$ improves the CAR in $57\%$ of
$100$ instances by $50\%$ and it achieves around $35\%$ at $99$th percentile.

(a) Facebook traffic traces on a small-scale network.  (b) Facebook traffic
traces on a large-scale network.

Figure 3: Average CAR with Facebook traces using (a) small-scale network and
(b) large-scale network. Each point in the x-axis represents network $[M,N]$.

(a) Synthetic traces.  (b) Facebook traces.

Figure 4: The $1$st-$10$th -$50$th-$90$th-$99$th percentiles of the average
gain in CAR with small-scale network $[10,60]$ using (a) synthetic and (b)
Facebook traces.

####  IV-B3 Prediction Error of $\mathtt{DCoflow}$

As mentioned at the end of Sec. III, we also evaluate the prediction error
$(|\sigma|-|\hat{\sigma}|)/|\sigma|$ of our heuristics, where
$\hat{\sigma}\subseteq\sigma$ is the set of coflows in $\sigma$ that satisfy
the deadline constraint after performing the actual resource allocation using
$\mathtt{GreedyFlowScheduling}$. Both variants of $\mathtt{DCoflow}$ provide a
prediction of CAR with an average error below $3.6\%$ for both traffic traces.

###  IV-C Online Setting

We now present a series of numerical results on the performance of the online
version of $\mathtt{DCoflow}$. The metric used for the performance evaluation
is the average CAR obtained over $40$ instances. In each instance, coflows
arrive sequentially according to a Poisson process of rate $\lambda$, i.e.,
the inter-arrival time of coflows is exponentially distributed with rate
$\lambda$. Unless stated otherwise, coflow priorities are computed upon
arrival of a new coflow ($f=\infty$).

As both versions of $\mathtt{DCoflow}$ provide similar results, we only
present the results obtained with $\mathtt{DCoflow\\_v1}$. The average CAR
obtained with $\mathtt{DCoflow\\_v1}$ is compared against those obtained with
the online version of Varys with deadline [16], $\mathtt{CS\text{-}MHA}$ and
$\mathtt{Sincronia}$. We investigate the effect of two main parameters: (i)
the coflow arrival rate $\lambda$ and (ii) the frequency $f$ at which coflow
priorities are updated.

####  IV-C1 Impact of Arrival Rate

We first study the impact of the arrival rate $\lambda$ on the CAR obtained
with the various algorithms. The CAR is averaged over $40$ instances, each one
with $4000$ coflow arrivals. The deadline of a coflow $k$ is drawn from a
uniform distribution in $[\text{CCT}_{k}^{0},4\text{CCT}_{k}^{0}]$. We also
consider two scenarios: a small fabric with $M=10$ machines, and a large
fabric with $M=50$ machines. For each scenario, results are presented for the
following values of $\lambda$: $\lambda=8$, $\lambda=12$, $\lambda=16$, and
$\lambda=20$.

Our results are shown in Figs. 5a and 5b, for the small fabric scenario and
the large one, respectively. We observe that $\mathtt{DCoflow\\_v1}$ achieves
a higher average CAR for all values of $\lambda$, and that the gain with
respect to the other scheduling algorithms increases with the value of
$\lambda$. If all algorithms achieve more or less the same CAR for a lightly
loaded fabric, $\mathtt{DCoflow\\_v1}$ clearly outperforms the other
algorithms when the fabric is highly congested.

(a) Small fabric.  (b) Large fabric.

Figure 5: Average CAR using synthetic traffic with varying $\lambda$ and (a)
$M=10$ and (b) $M=50$.

Figs. 6a and 6b show respectively the average CAR obtained with $M=10$ and
$M=100$, both with $4000$ coflows, using the Facebook dataset. Similar to what
obtained with the synthetic traffic traces, $\mathtt{DCoflow\\_v1}$
outperforms all other methods with significant gains. When the fabric is
highly congested (i.e., with $M=10$), again $\mathtt{DCoflow\\_v1}$ yields a
higher gain compared to other algorithms. For instance,
$\mathtt{DCoflow\\_v1}$ achieves $9.3\%$ higher CAR than $\mathtt{Sincronia}$
when $M=100$ (see Fig. 6b), while with $M=10$, the gap becomes $16.4\%$ (see
Fig. 6a).

(a) Small fabric.  (b) Large fabric.

Figure 6: Average CAR using Facebook traffic with varying $\lambda$ and (a)
$M=10$ and (b) $M=100$.

####  IV-C2 Impact of Update Frequency

We now evaluate the impact of the update frequency $f$ on the average CAR. We
consider the following values of $f$: $f=\frac{\lambda}{2}$, $f=\lambda$,
$f=2\lambda$, and $f=\infty$. Recall that $f=\infty$ means that priorities are
updated upon each coflow arrival. We assume that $M=10$ and compute the CAR by
averaging over $40$ instances. For each instance, we simulate $8,000$ coflow
arrivals according to a Poisson process at rate $\lambda$, assuming that the
deadline of a coflow $k$ is uniformly distributed in
$[\text{CCT}_{k}^{0},2\text{CCT}_{k}^{0}]$. We present the average CAR
obtained for different values of $f$
($f\in\\{\frac{\lambda}{2},\lambda,2\lambda,\infty\\}$) and for different
values of the arrival rate $\lambda$ (from $2$ to $10$).

(a) Without batch.  (b) With batch.

Figure 7: Average CAR of $\mathtt{DCoflow\\_v1}$ and $\mathtt{CS\text{-}MHA}$
using synthetic traffic with $[10,8000]$ and varying $\lambda$, when obtaining
(a) one single coflow per arrival; and (b) a random batch of coflow per
arrival.

The results in Fig. 7a are obtained from a simulation, in which each arrival
corresponds to one single coflow. We note, as before, that for a low arrival
rate, both algorithms provide a similar average CAR (for $\lambda=2$,
$\mathtt{CS\text{-}MHA}$ achieves a slightly higher CAR than
$\mathtt{DCoflow\\_v1}$), but that $\mathtt{DCoflow\\_v1}$ clearly outperforms
$\mathtt{CS\text{-}MHA}$ when the fabric is highly congested. We also note
that a higher frequency $f$ significantly improves the CAR for both
algorithms. For instance, for $\lambda=2$ (resp. $\lambda=10$), the average
CAR is increased by $52\%$ (resp. $46\%$) if we update coflow priorities upon
arrival of each new coflow instead of using the periodic scheme with
$f=\frac{\lambda}{2}$. These results suggest that there is a need for a trade-
off between the computational complexity of updating coflow priorities at a
high frequency and the CAR achieved. In Fig. 7b, we present a similar result,
but assuming that coflows arrive in batches. The batch size follows a uniform
distribution $\mathcal{U}([5,15])$. Since the average number of coflows in
each batch is $10$, in this setting, we divide the batch arrival rate by $10$
to obtain the same coflow arrival rates as in Fig. 7a. The results obtained
for batch arrivals are similar to those obtained previously, but we note that
the gains of $\mathtt{DCoflow\\_v1}$ with respect to $\mathtt{CS\text{-}MHA}$
are significantly higher in this case. Moreover, we note that the benefit of
using a higher update frequency is lower in this case (e.g., for $\lambda=10$,
the average CAR is increased only by $17\%$ if we use $f=\infty$ instead of
$f=\frac{\lambda}{2}$).

##  V Related Work

As discussed in Sec. I, most works in the literature focus on CCT
minimization, and deadline scheduling has received comparatively less
attention. $\mathtt{Varys}$ [9] was one of the first algorithms for deadline-
sensitive coflow scheduling. It uses a cascade of coflow admission control and
scheduling. The scheduler strives for CCT minimization combining (i) a coflow
ordering heuristic based on the per coflows bottleneck’s completion time; and
(ii) an allocation algorithm to assign bandwidth to individual flows of each
coflow. Rate allocation in $\mathtt{Varys}$ is performed to approximately
align the completion time of all coflows to the bottleneck one.

$\mathtt{Chronos}$ [16] is a heuristic for deadline scheduling which avoids
starvation for flows that do not meet their deadlines by granting them the
residual bandwidth. A priority order is determined first, and coflows are
hence allocated the minimum necessary bandwidth to meet their individual
deadlines. Once all the flows that meet their deadlines have been allocated
bandwidth, the residual bandwidth is shared by the remaining coflows in
proportion to their demands.

In [11], a connection between deadline scheduling of coflows and the well-
known problem of minimizing late jobs in a concurrent open shop—a known NP-
hard problem, is made. A heuristic based on the Moore-Hodgson’s algorithm [12]
for a single link is proposed. Both centralized as well as decentralized
heuristics are introduced (namely $\mathtt{CS\text{-}MHA}$ and
$\mathtt{D\textsuperscript{2}\text{-}CAS}$, respectively).

A formal description of the deadline scheduling problem including bandwidth
allocation of flows was given in [8]. The CDS maximization problem is
formulated as an MILP (called $\mathtt{CDS\text{-}LP}$). Time is divided into
intervals whose boundaries are the coflow deadlines arranged in increasing
order. The program determines which coflows to accept and the amount of
bandwidth to allocate in each interval. $\mathtt{CDS\text{-}LP}$ is shown to
be NP-hard, and an approximation based on LP relaxation (called
$\mathtt{CDS\text{-}LPA}$) of the binary variables is also proposed.
$\mathtt{CDS\text{-}LPA}$ only accepts coflows for which the relaxed variable
is strictly equal to $1$, i.e., only coflows that are completely accepted by
the LP relaxation are retained.

For completeness, we also cite a few works considering the minimization of CCT
[4, 6, 17, 5, 18] as well as the survey article [19]. Popular among CCT
minimization algorithms is $\mathtt{Sincronia}$ proposed in [13]. It considers
scheduling on the network bottlenecks and returns a scheduling coflow order
achieving a $4$–approximation factor.

##  VI Conclusion and Future Work

In this paper, we introduced a new joint coflow admission control and
scheduling algorithm for a batch of coflows with deadlines. The proposed
schemes leverage results from open-shop scheduling to determine a subset of
coflows to schedule and a corresponding $\sigma$-order which is then employed
in order to schedule coflows in priority.

Numerical results show that on small-scale networks, our algorithms perform
similar to or better than other deadline-sensitive algorithms proposed in the
literature. On large-scale networks, however, it shows significant
improvements with respect to existing algorithms, e.g., $98\%$ higher CAR than
$\mathtt{CS\text{-}MHA}$ in an offline setting. Our scheme also has a low
prediction error: even though the admission control is performed using a
bottleneck approximation for the CCT, almost all accepted coflows actually
finish within their deadline when being actually scheduled.

This behaviour is observed in both offline and online settings, both with
synthetic traces and for real traces from the Facebook data set. This shows
that the algorithm is robust with respect to the coflow size distribution and
performs very well across a wide range of network sizes.

Several extensions of this research line are possible. In future works, we
shall study the performance of our online solution in the case when coflows
tend to be released in batches, e.g., when a distributed framework polls
worker nodes with given period. Furthermore, a relevant case is that of
incomplete information on flow volumes, i.e., because the volume of a flow is
not directly available to the scheduler but only, for instance, via some a
priori distribution. Finally, issues of starvation and fairness issue among
coflows represent interesting issues we have not addressed yet.

## References

  * [1] M. Chowdhury and I. Stoica, “Coflow: A networking abstraction for cluster applications,” in _Proc. ACM HotNets_ , Redmond, Washington, 2012, p. 31–36.
  * [2] J. Dean and S. Ghemawat, “MapReduce: Simplified data processing on large clusters,” _Commun. ACM_ , vol. 51, no. 1, pp. 107–113, 2008.
  * [3] M. Zaharia, M. Chowdhury, M. J. Franklin, S. Shenker, I. Stoica _et al._ , “Spark: Cluster computing with working sets.” _HotCloud_ , vol. 10, no. 10-10, p. 95, 2010.
  * [4] N. M. K. Chowdhury, “Coflow: A networking abstraction for distributed data-parallel applications,” Ph.D. dissertation, University of California, Berkeley, 2015.
  * [5] M. Shafiee and J. Ghaderi, “An improved bound for minimizing the total weighted completion time of coflows in datacenters,” _IEEE/ACM Trans. Netw._ , vol. 26, no. 4, pp. 1674–1687, 2018.
  * [6] S. Agarwal, S. Rajakrishnan, A. Narayan, R. Agarwal, D. Shmoys, and A. Vahdat, “Sincronia: Near-optimal network design for coflows,” in _Proc. ACM SIGCOMM_ , 2018, pp. 16–29.
  * [7] M. Chowdhury _et al._ , “Near optimal coflow scheduling in networks,” in _Proc. ACM SPAA_ , Phoenix, AZ, USA, June 22-24 2019, pp. 123–134.
  * [8] S.-H. Tseng and A. Tang, “Coflow deadline scheduling via network-aware optimization,” in _Proc. Annu. Allert. Conf. Commun. Control Comput._ , 2018, pp. 829–833.
  * [9] M. Chowdhury, Y. Zhong, and I. Stoica, “Efficient Coflow Scheduling with Varys,” in _Proc. ACM SIGCOMM_ , 2014, pp. 443–454.
  * [10] B. Lin and A. Kononov, “Customer order scheduling to minimize the number of late jobs,” _Eur. J. Oper. Res._ , vol. 183, no. 2, pp. 944–948, 2007.
  * [11] S. Luo, H. Yu, and L. Li, “Decentralized deadline-aware coflow scheduling for datacenter networks,” in _Proc. IEEE ICC_ , 2016, pp. 1–6.
  * [12] J. M. Moore, “An n job, one machine sequencing algorithm for minimizing the number of late jobs,” _Manag. Sci._ , vol. 15, no. 1, pp. 102–109, 1968\.
  * [13] S. Agarwal, R. Agarwal, S. Rajakrishnan, D. Shmoys, A. Narayan, and A. Vahdat, “Sincronia: Near-Optimal Network Design for Coflows,” in _Proc. ACM SIGCOMM_ , 2018, pp. 16–29.
  * [14] M. Mastrolilli, M. Queyranne, A. S. Schulz, O. Svensson, and N. A. Uhan, “Minimizing the sum of weighted completion times in a concurrent open shop,” _Oper. Res. Lett._ , vol. 38, no. 5, pp. 390–395, 2010.
  * [15] M. Chowdhury, “Coflow: A networking abstraction for distributed data-parallel applications,” Ph.D. dissertation, University of California, Berkeley, Nov. 2015\.
  * [16] S. Ma, J. Jiang, B. Li, and B. Li, “Chronos: Meeting Coflow Deadlines in Data Center Networks,” in _Proc. IEEE ICC_ , 2016.
  * [17] L. Chen, W. Cui, B. Li, and B. Li, “Optimizing coflow completion times with utility max-min fairness,” in _Proc. IEEE INFOCOM_ , 2016, pp. 1–9.
  * [18] L. Shi, Y. Liu, J. Zhang, and T. Robertazzi, “Coflow scheduling in data centers: routing and bandwidth allocation,” _IEEE Trans. Parallel Distrib. Syst._ , vol. 32, no. 11, pp. 2661–2675, 2021.
  * [19] S. Wang, J. Zhang, T. Huang, J. Liu, T. Pan, and Y. Liu, “A survey of coflow scheduling schemes for data center networks,” _IEEE Commun. Mag._ , vol. 56, no. 6, pp. 179–185, 2018.


Further author information: (Send correspondence to Maxime Quesnel)
E-mail<EMAIL_ADDRESS>
# A Simulator-based Autoencoder
for Focal Plane Wavefront Sensing

Maxime Quesnel  Montefiore Institute of Electrical Engineering and Computer
Science, University of Liège  Space sciences, Technologies and Astrophysics
Research (STAR) Institute, University of Liège  Gilles Orban de Xivry  Space
sciences, Technologies and Astrophysics Research (STAR) Institute, University
of Liège  Olivier Absil  Space sciences, Technologies and Astrophysics
Research (STAR) Institute, University of Liège  Gilles Louppe  Montefiore
Institute of Electrical Engineering and Computer Science, University of Liège

###### Abstract

Instrumental aberrations strongly limit high-contrast imaging of exoplanets,
especially when they produce quasi-static speckles in the science images. With
the help of recent advances in deep learning, we have developed in previous
works an approach that applies convolutional neural networks (CNN) to estimate
pupil-plane phase aberrations from point spread functions (PSF). In this work
we take a step further by incorporating into the deep learning architecture
the physical simulation of the optical propagation occurring inside the
instrument. This is achieved with an autoencoder architecture, which uses a
differentiable optical simulator as the decoder. Because this unsupervised
learning approach reconstructs the PSFs, knowing the true phase is not needed
to train the models, making it particularly promising for on-sky applications.
We show that the performance of our method is almost identical to a standard
CNN approach, and that the models are sufficiently stable in terms of training
and robustness. We notably illustrate how we can benefit from the simulator-
based autoencoder architecture by quickly fine-tuning the models on a single
test image, achieving much better performance when the PSFs contain more noise
and aberrations. These early results are very promising and future steps have
been identified to apply the method on real data.

###### keywords:

Deep learning, convolutional neural networks, autoencoders, focal plane
wavefront sensing, phase retrieval

##  1 INTRODUCTION

Imaging exoplanets requires instruments to have the capability of probing high
contrasts between planets and their host star. Telescopes incorporating a
coronagraph and adaptive optics (AO) tremendously help in this matter, but
residual wavefront aberrations, in particular non-common path aberrations
(NCPA) between the scientific and wavefront sensing arms, still greatly limit
the detection of exoplanets[1]. Focal-plane wavefront sensing (FPWFS)
methods[2] are particularly beneficial to measure NCPAs since the signature of
these aberrations can be found in the focal-plane images.

Deep learning with convolutional neural networks (CNN) has recently been
implemented for focal-plane wavefront sensing[3, 4, 5, 6, 7, 8] with good
results. CNNs offer great performance and robustness as well as fast
inference, making them particularly adequate for NCPA correction during
observing nights. One limitation of supervised learning using CNNs however is
that the models require a ground truth to be trained. As building models with
on-sky data is particularly challenging due to the lack of precise knowledge
about the NCPAs themselves, it may be difficult to build reliable labelled
datasets. Unsupervised learning approaches are thus particularly compelling,
since they work without knowing the true phase aberrations. By reconstructing
the input data, deep learning autoencoders are a well proven unsupervised
approach. For focal-plane wavefront sensing the targets are the NCPAs instead
of the point spread functions (PSF), and one way to extract this phase
information with an autoencoder is to add physical knowledge to the model.

Several works in the field of optics have incorporated a physical simulator
into their deep learning architecture. Bostan et al. 2020[9] and Wang et al.
2020[10] for instance proposed each an unsupervised approach with an untrained
neural network that reconstructs images on-the-fly based on an optical model.
Peng et al. 2020 [11] also incorporates an optical wave propagation model in
the architecture to improve the quality of holographic images. Automatic
differentiation of optical models has also been performed in the context of
wavefront aberrations, notably for PSF reconstruction (Liaudat et al.
2022[12]) as well as for phase retrieval and design (Wong et al. 2021[13]).

Motivated by the benefits a differentiable simulator can offer for an
unsupervised approach, we define an autoencoder architecture that contains a
CNN as the encoder and a optical simulator as the decoder. This allows to
constrain the encoder to predict Zernike coefficients corresponding to the
probed NCPAs. The method is described in Sec. 2. This unsupervised approach is
tested on simulated data, and the results are presented in Sec. 3. We notably
evaluate how the method handles photon noise and atmospheric turbulence
residuals compared to using only a CNN. The possibility of fitting the models
on an observed PSF that lies outside the training data distribution is also
investigated. Developments for the approach to work with real instruments are
finally discussed in Sec. 4.

##  2 SIMULATOR-BASED AUTOENCODER

###  2.1 Data generation

For the present study we consider a simple optical propagation between a pupil
plane and a focal plane, so that the PSF is expressed as:

$PSF(x,y)=|\mathscr{F}[A(x,y)\,e^{\,i\,\Theta(x,y)}]|^{2},$ (1)

where $A(x,y)$ is the pupil function and $\Theta(x,y)$ the pupil-plane phase
aberrations.

We construct the phase maps using Zernike polynomials following the Noll
convention[14], starting from the tip mode, and based on an annular aperture
with a central obstruction of 30% of the total diameter. The sets of Zernike
coefficients are randomly generated to approximate a $1/f^{2}$ power spectral
density profile, and scaled to a given median root-mean-square (rms) wavefront
error (WFE). In our simulations we consider two aberration regimes: NCPAs
around 70 nm rms covering 20 Zernike modes, and 350 nm rms over 100 modes. The
PSFs are computed using the PROPER[15] optical propagation package. Our data
is generated in the $K$ band at $\lambda=2200$ nm, while we consider an
aperture diameter of 10 m, a pixel scale of $\simeq 11.43$ mas/pix, and a
field-of-view of $\simeq 1.47"$. The phase maps as well as the PSFs contain
$64\times 64$ pixels. In addition to the in-focus PSFs, we also generate out-
of-focus PSFs with a defocus of $\lambda/4=550$ nm rms, which are used to lift
the sign ambiguity that exists for the modes of even radial order, as
typically done for focal-plane wavefront sensing[16, 17]. We add photon noise
to the PSFs so that SNR $=\sqrt{N_{ph}}$, with $N_{ph}$ the number of photons.
We also apply a square-root stretching operation to the PSFs, before
normalizing the flux in the range [0,1]. This operation is done inside the
encoder for the autoencoder architecture.

For the experiments done in Sec. 3.2.4, we also generate data containing
atmospheric turbulence residuals. This is achieved by simulating an extreme AO
system using the COMPASS library[18], giving a wavefront error of about 50 nm
rms. To simulate a 1-s exposure in presence of a given amount of static NCPA,
we use a sequence of 10 consecutive AO phase screens, and we sum up the
corresponding PSFs. These AO simulations are described in more details in
Quesnel et al. 2022[8] (submitted).

###  2.2 Autoencoder architecture

Figure 1: The proposed simulator-based autoencoder architecture (SimAE), using
a CNN as the encoder and a differentiable optical simulator as the decoder.
The latent space contains the Zernike coefficients which are used during
training to reconstruct noiseless PSFs, while during inference these
coefficients can be applied for NCPA correction.

In this work we propose an autoencoder architecture, which is composed of an
encoder and a decoder. The encoder is a CNN that maps Zernike coefficients
from PSFs $x$, i.e., it approximates a non-linear function $f$ such that
$z\approx f(x)$. We choose to use a state-of-the-art deep CNN called
EfficientNet [19]. This type of architecture stands out from other ones by
using a new scaling technique, with all dimensions of the CNN (depth, width,
and resolution) being scaled similarly. We choose to use the EfficientNet-B4
flavour, which offers the best trade-off between performance and runtime.

For our unsupervised approach, we also add a decoder which corresponds to a
differential optical simulator, i.e. the decoder reconstructs the input PSFs
based on the predicted Zernike coefficients, allowing to work without labels.
This decoder is the same as the data generator defined in Eq. 1, with the
exception that it produces noiseless PSFs. All the simulator’s parameters are
here fixed, while it may be beneficial in the future to consider learning
parameters of the simulator for in-lab and on-sky applications. Regarding the
latent space that exists between the encoder and the decoder, the number of
Zernike modes has to be predefined. During inference on a trained model, the
encoder containing the CNN allows to obtain the Zernike coefficients that
correspond to the NCPAs of interest. Fig. 1 summarizes our simulator-based
autoencoder architecture, which we call “SimAE” in the following sections.

The loss function, which is used to optimize the neural network’s weights in
order to fit the model to the data, needs to be adapted to our SimAE
architecture. Since the decoder reconstructs noiseless PSFs, the typical mean-
square-error loss that makes the assumption of Gaussian noise is not
appropriate on the PSF residuals. Since we are working on a simplified case
where only photon noise is present in the images, i.e. each pixel follows a
Poisson distribution, we can define our loss so that the log-likelihood of the
Poisson probability distributions is maximized, giving the expression:

$\mathcal{L}_{SimAE}(x;\phi)=-\mathbb{E}_{x\sim
p(x)}\left[\log\left(\frac{\lambda(x;\phi)^{x}}{x!}\exp(-\lambda(x;\phi))\right)\right],$
(2)

where $\lambda(x;\phi)$ are the rates of the Poisson distributions, which
depend on both the input PSFs $x$ sampled from the dataset distribution
$p(x)$, and the parameters $\phi$ of the encoder (CNN). These rates correspond
to the noiseless PSFs that are constructed by the decoder.

The baseline on which we compare our SimAE models is a standard CNN
architecture. The loss function is in this case defined as the root-mean-
square error of the phase residuals. It can be expressed using the Zernike
coefficients as:

$\mathcal{L}_{CNN}(z,\widehat{z}(x;\phi))=\sqrt{\frac{1}{N}{\sum_{i}^{N}(z_{i}-\widehat{z_{i}}(x;\phi))^{2}}},$
(3)

where $z$ and $\widehat{z}$ correspond to the true and predicted Zernike
coefficients respectively, while $N$ is the number of Zernike modes (20 or
100).

##  3 EXPERIMENTS

###  3.1 Protocol

Figure 2: Examples of training and validation losses over the epochs for each
approach: CNN (left) and SimAE (right). In these examples the models are
trained on datasets containing 350 nm rms of aberrations distributed over 100
modes, and a SNR equal to 300.

For a given training, a dataset composed of $10^{5}$ PSF pairs is randomly
split into training (90%) and validation (10%) sets, with a batch size of 64.
For the supervised learning approach with a CNN, each sample also contains the
true NCPA phase maps as labels. The Adam optimizer [20] with a starting
learning rate of $10^{-3}$ is used to update the weights of the CNN. The
learning rate is decreased by a factor two as soon as the validation loss
reaches a plateau over 15 epochs to improve the performance. Efficientnet-b4
weights are initialized with models pretrained on ImageNet. An example of the
resulting training and validation losses for each method can be found in Fig.
2. The models are then systematically evaluated on 1000 new samples.

###  3.2 Results

####  3.2.1 Performance compared to a standard deep CNN

We first compare the performance regarding each Zernike mode. The metric used
is the root-mean-square error (RMSE) per mode, computed over the entire test
set:

$\sigma_{z}=\sqrt{\frac{1}{N_{\rm test}}\,\sum_{i}^{N_{\rm
test}}(\hat{c}_{i}-c_{i})^{2}}\>,$ (4)

where $N_{\rm test}$ is the number of test samples, while $\hat{c}$ and $c$
are the estimated and true Zernike coefficients, respectively.

As shown in Fig. 3, the Zernike coefficients are very well reconstructed, even
at high aberrations contents (350 nm rms over 100 modes). There is no
particular distinction between even and odd modes, and the SimAE and CNN
manifest very similar performance.

Figure 3: Performance per Zernike mode, for a low aberration regime (70 nm
over 20 modes, light shades) and a high aberration regime (350 nm over 100
modes, darker shades), at a SNR of 1000. The CNN and SimAE architectures are
compared (blue and red, respectively).

Since the SimAE is defined to take into account photon noise through its loss
function, it is worth assessing the performance for different SNR levels. We
thus compare the residual errors in terms of rms WFE on the total phase
residuals. The metric is therefore defined for each test sample as:

$\sigma_{\phi}=\sqrt{\frac{1}{N_{\rm pix}}\sum_{i}^{N_{\rm
pix}}(\hat{\phi}_{i}-\phi_{i})^{2}},$ (5)

where $N_{\rm pix}$ is the number of pixels within the pupil area, while
$\hat{\phi}$ and $\phi$ are the estimated and true pupil phases, respectively.

The errors on the phase retrieval for a range of SNR levels are shown in Fig.
4, where each point corresponds to a model trained and evaluated at the same
given SNR. The two approaches show almost identical performance, and a plateau
is reached around a SNR of 3000 in the high aberration regime. This has been
observed and discussed in previous works[7, 8] as well. Only at very low SNRs
the SimAE appear to differ from the CNN approach. This difference could be due
to the loss function of the SimAE that is defined on the PSFs and not on the
NCPAs, so that there is not the implicit regularization occurring with
supersized CNNs, which has been described by Orban de Xivry et al. 2021[7].
The residual errors thus tend to follow the theoretical limit for lower SNRs
with the SimAE.

Figure 4: Performance at different SNR levels for the SimAE (blue) and the
standard CNN (red) models, at low and high aberration regimes (darker and
brighter shades respectively). The points correspond to the median value and
the error bars to the 2-98th percentiles.

####  3.2.2 Robustness of the trained models

Deep learning models may show limitations when evaluated on data lying outside
the training data distribution. It is therefore important to assess the
stability of the models in such conditions, which would probably occur when
applying the algorithm on sky.

Here we compare again the SimAE method with the CNN, this time varying the
testing SNR (Fig. 5, left), on a model trained at a median WFE of 70 nm rms
distributed over 20 modes. Giving data containing less noise to the models
always provides stable predictions, and we even obtain a marginal improvement
with the SimAE when trained at a SNR of 30. If we introduce more noise to the
test data, the error quickly increases, at about the same rate with both
approaches. We also experiment by changing the testing input WFE in the data
(Fig. 5, right). Here the CNN and SimAE showcase identical behaviour, with
constant performance for lower testing WFE than the training one, and strong
degradation for higher testing WFE, although with the SimAE the performance
does not degrade as quickly as with the CNN. Similar trends regarding the CNN
has been found in previous works[7, 8]. Overall the SimAE method offers only
slight improvement in terms of robustness compared to the CNN. To benefit the
most from such a physics-based approach, it would be useful to include more
knowledge about the data in the models, which could be provided by a
variational autoencoder[21] approach for instance.

Figure 5: Left: Performance with changing the SNR of the test data for two
models trained on SNRs of 30 and 300 (indicated by the vertical dotted lines).
The median input WFE is 70 nm with 20 modes. Right: Performance for different
testing input WFE, on models trained at 70 nm and 350 nms (20 and 100 modes
respectively).

####  3.2.3 Transfer learning

Figure 6: Performance of models fine-tuned on each test sample individually,
for varying SNR (left) and input WFE levels (right). The model used for both
plots was initially trained on data with 70 nm rms, 20 modes, and a SNR of
100, as indicated by the vertical dotted lines.

Instead of directly applying a trained model on observed data, we propose to
fine-tune the model using a single “on-sky” image, initializing with the
trained weights of the original model. This is done to illustrate how we could
quickly adapt the model from simulated to on-sky data, which may slightly
differ from the training data. To test this on our SimAE method, we use a
model trained on data containing 70 nm rms and a SNR of 100. The levels of
NCPA and SNR are changed on the test data, and the results are shown in Fig.
6. The models perform much better when transfer learning is done on each test
sample, in particular when their SNR and WFE are higher than for the original
training data. This second training on a single test sample is performed over
200 epochs, which takes about 20 seconds to complete on one GPU (Nvidia
GeForce RTX 2080 Ti). For data relatively close to the ones used for pre-
training, convergence is naturally achieved after much less epochs and the
final prediction can be made in only a few seconds.

####  3.2.4 Atmospheric turbulence residuals

In the previous experiments, the decoder, i.e the differentiable optical
simulator, was able to exactly reconstruct what our data generator produced to
build our training and test sets, except for photon noise, but for which the
loss function was adapted. In order to test the method on data the decoder
cannot totally reconstruct, we introduce atmospheric turbulence residuals to
the input data, without giving this information to the decoder. The way these
AO residuals are added is described in Sec. 2.1. Even with the presence of
this additional source of noise, the training is very stable, and the
evaluation results shown in Fig. 7 indicate that the SimAE suffers a
reasonable degradation in performance compared to a case without turbulence.
Using a CNN with labels that does not account for these AO residuals gives the
same results.

Figure 7: Residual errors as a function of the input WFE, for a model trained
on a dataset with WFEs distributed around 70 nm (20 modes). A SimAE model with
data containing only photon noise (cyan), is compared to CNN and SimAE models
(purple and blue, respectively) trained on data containing also AO residuals,
which are not provided to the decoder during training and evaluation.

##  4 CONCLUSIONS AND FUTURE WORK

We have developed an unsupervised deep learning method for focal-plane
wavefront sensing by incorporating an optical simulator into a CNN-based
autoencoder architecture (SimAE). We observe very similar performance compared
to using a simple CNN, both approaches reaching the expected theoretical limit
over a large range of conditions. The SimAE however has the major advantage of
not requiring labels during the training of the models. One promising
opportunity with an approach like this SimAE is that we can quickly fit the
model to the currently observed PSF using a pre-trained model. We have shown
in this paper that such fine-tuning of the models can greatly improve the
results for data showing different noise or NCPA levels. Dealing with data the
decoder cannot reproduce is naturally a potential limiting factor for the
method. We have trained models with atmospheric turbulence residuals present
on top of NCPAs, and the loss of performance that we obtain stays well
constrained. In practice these turbulence residuals may even be given to the
decoder, by using the wavefront sensor telemetry to learn what the AO
residuals effectively are, with the help of a small neural network for
example.

Several components of our architecture can be improved in the future. We may
first benefit from variational inference, i.e. predicting posterior
distribution for Zernike coefficients as well as to incorporate prior
knowledge about them. This could yield better robustness for the models thanks
to the additional information given during training and inference. Regarding
the decoder, a more complete simulator, e.g. including a coronagraph for high-
contrast imaging and other relevant instrumental components, will be essential
to envision applying the technique on-sky with a real instrument.

###### Acknowledgements.

For our study we have used the optical propagation package PROPER[15] to
generate our datasets. We have trained our deep learning models using the
PyTorch[22] library while Sacred[23] was used to manage our experiments. The
research was supported by the the Wallonia-Brussels Federation (grant for
Concerted Research Actions), and by the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation program (grant
agreement No 819155).

## References

  * [1] Guyon, O., “Extreme Adaptive Optics,” Annual Review of Astronomy and Astrophysics 56, 315–355 (Sept. 2018).
  * [2] Jovanovic, N., Absil, O., Baudoz, P., Beaulieu, M., Bottom, M., Cady, E., Carlomagno, B., Carlotti, A., Doelman, D., Fogarty, K., Galicher, R., Guyon, O., Haffert, S., Huby, E., Jewell, J., Keller, C., Kenworthy, M. A., Knight, J., Kühn, J., Miller, K., Mazoyer, J., N’Diaye, M., Por, E., Pueyo, L., Riggs, A. J. E., Ruane, G., Sirbu, D., Snik, F., Wallace, J. K., Wilby, M., and Ygouf, M., “Review of high-contrast imaging systems for current and future ground-based and space-based telescopes: Part II. Common path wavefront sensing/control and coherent differential imaging,” in [Adaptive Optics Systems VI ], Close, L. M., Schreiber, L., and Schmidt, D., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 10703, 107031U (July 2018).
  * [3] Paine, S. W. and Fienup, J. R., “Machine learning for improved image-based wavefront sensing,” Optics Letters 43, 1235 (Mar. 2018).
  * [4] Andersen, T., Owner-Petersen, M., and Enmark, A., “Neural networks for image-based wavefront sensing for astronomy,” Optics Letters 44, 4618 (Sept. 2019).
  * [5] Andersen, T., Owner-Petersen, M., and Enmark, A., “Image-based wavefront sensing for astronomy using neural networks,” Journal of Astronomical Telescopes, Instruments, and Systems 6, 034002 (July 2020).
  * [6] Quesnel, M., Orban de Xivry, G., Louppe, G., and Absil, O., “Deep learning-based focal plane wavefront sensing for classical and coronagraphic imaging,” in [Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 11448, 114481G (Dec. 2020).
  * [7] Orban de Xivry, G., Quesnel, M., Vanberg, P. O., Absil, O., and Louppe, G., “Focal plane wavefront sensing using machine learning: performance of convolutional neural networks compared to fundamental limits,” Monthly Notices of the Royal Astronomical Society 505, 5702–5713 (Aug. 2021).
  * [8] Quesnel, M., Orban de Xivry, G., Absil, O., and Louppe, G., “A deep learning approach for focal-plane wavefront sensing using vortex phase diversity,” (submitted, 2022).
  * [9] Bostan, E., Heckel, R., Chen, M., Kellman, M., and Waller, L., “Deep phase decoder: self-calibrating phase microscopy with an untrained deep neural network,” Optica 7, 559–562 (Jun 2020).
  * [10] Wang, F., Bian, Y., Haichao, W., Lyu, M., Pedrini, G., Osten, W., Barbastathis, G., and Situ, G., “Phase imaging with an untrained neural network,” Light: Science & Applications 9, 77 (05 2020).
  * [11] Peng, Y., Choi, S., Padmanaban, N., and Wetzstein, G., “Neural holography with camera-in-the-loop training,” ACM Trans. Graph. 39 (nov 2020).
  * [12] Liaudat, T., Starck, J.-L., Kilbinger, M., and Frugier, P.-A., “Rethinking data-driven point spread function modeling with a differentiable optical model,” (2022).
  * [13] Wong, A., Pope, B., Desdoigts, L., Tuthill, P., Norris, B., and Betters, C., “Phase retrieval and design with automatic differentiation: tutorial,” Journal of the Optical Society of America B 38, 2465 (aug 2021).
  * [14] Noll, R. J., “Zernike polynomials and atmospheric turbulence.,” Journal of the Optical Society of America (1917-1983) 66, 207–211 (Mar. 1976).
  * [15] Krist, J. E., “PROPER: an optical propagation library for IDL,” in [Optical Modeling and Performance Predictions III ], Kahan, M. A., ed., 6675, 250 – 258, International Society for Optics and Photonics, SPIE (2007).
  * [16] Gonsalves, R. A., “Phase Retrieval And Diversity In Adaptive Optics,” Optical Engineering 21(5), 829 – 832 (1982).
  * [17] Vievard, S., Bos, S. P., Cassaing, F., Ceau, A., Guyon, O., Jovanovic, N., Keller, C., Lozi, J., Martinache, F., Mary, D., Montmerle-Bonnefois, A., Mugnier, L., N'diaye, M., Norris, B., Sahoo, A., Sauvage, J.-F., Snik, F., Wilby, M. J., and Wong, A., “Overview of focal plane wavefront sensors to correct for the Low Wind Effect on SUBARU/SCExAO,” in [6th International Conference on Adaptive Optics for Extremely Large Telescopes, AO4ELT 2019 ], (June 2019).
  * [18] Ferreira, F., Gratadour, D., Sevin, A., and Doucet, N., “Compass: An efficient gpu-based simulation software for adaptive optics systems,” 2018 International Conference on High Performance Computing & Simulation (HPCS) , 180–187 (2018).
  * [19] Tan, M. and Le, Q., “EfficientNet: Rethinking model scaling for convolutional neural networks,” in [Proceedings of the 36th International Conference on Machine Learning ], Chaudhuri, K. and Salakhutdinov, R., eds., Proceedings of Machine Learning Research 97, 6105–6114, PMLR (09–15 Jun 2019).
  * [20] Kingma, D. P. and Ba, J., “Adam: A method for stochastic optimization,” (2017).
  * [21] Kingma, D. P. and Welling, M., “Auto-encoding variational bayes,” (2013).
  * [22] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S., “Pytorch: An imperative style, high-performance deep learning library,” in [Advances in Neural Information Processing Systems ], Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., eds., 32, 8026–8037, Curran Associates, Inc. (2019).
  * [23] Klaus Greff, Aaron Klein, Martin Chovanec, Frank Hutter, and Jürgen Schmidhuber, “The Sacred Infrastructure for Computational Research,” in [Proceedings of the 16th Python in Science Conference ], Katy Huff, David Lippa, Dillon Niederhut, and Pacer, M., eds., 49 – 56 (2017).