# Result and Congestion Aware Optimal Routing and Partial Offloading in
Collaborative Edge Computing

Jinkun Zhang, Yuezhou Liu, Edmund Yeh  Department of Electrical and Computer
Engineering, Northeastern University, Boston, US
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Collaborative edge computing (CEC) is an emerging paradigm where heterogeneous
edge devices (stakeholders) collaborate to fulfill computation tasks, such as
model training or video processing, by sharing communication and computation
resources. Nevertheless, the optimal data/result routing and computation
offloading strategy in CEC with arbitrary topology still remains an open
problem. In this paper, we formulate a partial-offloading and multi-hop
routing model for arbitrarily divisible tasks. Each node individually decides
the computation of the received data and the forwarding of data/result
traffic. In contrast to most existing works, our model applies for tasks with
non-negligible result size, and enables separable data sources and result
destination. We propose a network-wide cost minimization problem with
congestion-aware cost to jointly optimize routing and computation offloading.
This problem covers various performance metrics and constraints, such as
average queueing delay with limited processor capacity. Although the problem
is non-convex, we provide non-trivial necessary and sufficient conditions for
the global-optimal solution, and devise a fully distributed algorithm that
converges to the optimum in polynomial time, allows asynchronous individual
updating, and is adaptive to changes in network topology or task pattern.
Numerical evaluation shows that our proposed method significantly outperforms
other baseline algorithms in multiple network instances, especially in
congested scenarios.

##  I Introduction

In recent years, we are experiencing an explosive increment in the number of
mobile and IoT devices. Our daily life is increasingly involved with new
mobile applications for work, entertainment, social networking, health care,
etc. Many of them are computation-intensive and time-critical, such as AR/VR
and autonomous driving. Meanwhile, mobile devices running these applications
generate a huge amount of data traffic which is predicted to reach 288EB per
month in 2027 [1]. It becomes impractical to direct all computation requests
and their data to central cloud due to limited backhaul bandwidth and high
latency. Edge computing is then proposed as a promising solution to provide
computation resources and cloud-like services in close proximity to the mobile
devices.

In edge computing, requesters offload their computation to the edge servers. A
new concept, extending this idea, is the so-called collaborative edge
computing (CEC). Besides point-to-point offloading, CEC let multiple
stakeholders (mobile devices, IoT devices, edge servers, or cloud) collaborate
with each other by sharing data, communication resources, and computation
resources to finish computation tasks [2]. On the one hand, CEC improves the
utilization efficiency of resources so that computation-intensive and time-
critical services can be better completed at the edge. Mobile devices equipped
with computation capabilities can collaborate with each other through D2D
communication [3]. Edge servers can also collaborate with each other for load
balancing or further with central cloud to offload demands that they cannot
accommodate [4]. On the other hand, CEC is needed when there is no direct
connection between devices and edge servers. Consider unmanned aerial vehicle
(UAV) swarms or autonomous cars in rural areas, computation-intensive tasks of
UAVs or cars far away from the wireless access point should be collaborated
computed or offloaded through multi-hop routing to the edge server with the
help of other devices [5, 3].

We wish to study a general framework of CEC which enables varies types of
collaboration among stakeholders. In particular, we consider a multi-hop
network with arbitrary topology, where the nodes collaboratively finish
multiple computation tasks. Nodes have heterogeneous computation capabilities
and some are also data sources (sensors and mobile users) that generate data
for computation tasks. Each task has a requester node for the computation
result. We allow partial offloading, where the computation of a task can be
partitioned and conducted by multiple nodes, by considering data-driven and
divisible tasks, i.e., a task can be divided arbitrarily into sub-tasks that
compute different parts of the input data. Finishing a task requires the
routing of data from possibly multiple data sources to multiple nodes for
computation and the routing of results to the requester. We aim for a joint
routing (how to route the data/result) and computation offloading (where to
compute) strategy for all tasks that minimizes the total communication and
computation costs. This problem is non-trivial, but we are able to find the
optimal joint routing and computation strategy in the proposed framework: In
particular, we provide the necessary and sufficient conditions for the global
optimum and propose an distributed and adaptive algorithm that converges to
the global optimum with small information exchange overhead.

Compared with existing edge computing and CEC studies, our work is distinct in
several aspects. First, while most works study only offloading-type
computation where requesters offload their data for computation [5, 6, 7], we
allow data sources and requesters be different nodes to cover computation that
requires fetching data from other nodes (like sensors) besides the requester
itself. Sahni et al. [3, 8] also consider arbitrary data sources, but are
restricted to fully connected networks or predefined routing protocols.
Second, while we study routing for both data and result flows, most existing
works ignore the cost for transmitting results [9, 2, 10] and simply let
results transmitted along the reverse path of data [11, 12]. However, the size
of computation result is not negligible in many applications. For example, the
communication cost of intermediate results is one major concern in federated
or distributed machine learning [13], and result size can even be larger than
data size in applications like decompression and image enhancement. Third, we
model congestion in our framework by considering non-linear communication and
comuputation costs instead of linear costs as in [5, 14]. Finally, the
existence of optimal solution, and distributed and adaptive algorithm makes
our framework possible to be implemented in practical distributed edge
networks with provable performance guarantee and robustness to network
condition changes.

Our detailed contributions are as follows:

  * •

To the best of our knowledge, we formulate the fist framework that jointly
studies partial offloading and routing for both data and result in arbitrary
CEC networks with congestible links.

  * •

We provide the global optimal routing and offloading strategy for this non-
convex problem, by studying the necessary and sufficient optimality condition.

  * •

We devise a distributed, asynchronous, and adaptive algorithm that converges
to the global optimum.

  * •

By extensive experiments, we show the advantages of the proposed algorithm
over several baselines in different network topologies.

The remainder of this paper is organized as follows. We present our model of a
CEC network and formulate the optimization problem in Section II. The
necessary and sufficient optimality conditions are discussed in Section III,
while the proposed algorithm is presented in Section IV. We summarize our
experiment results in Section V and finally conclude in Section VI.

##  II Network Model and Problem Formulation

We begin by presenting our formal model of a collaborative edge computing
network where multiple stakeholders collaborate to finish computation tasks.
Such network is motivated by several real-word applications such as IoT
networks, connected vehicles and UAV swarms. An example system that involves
IoT network on the edge is shown in Fig. 1.

###  II-A Network model

We consider a quasi-static network, represented by a directed and strongly
connected graph $G=(V,E)$ where $V$ is the set of nodes (devices with either
task, data or computation resource to share) and $E$ is the set of links. For
nodes $i,j\in V$, there exists link $(i,j)\in E$ if a feasible connection from
$i$ to $j$ is available for data or result forwarding. Denote $\mathcal{I}(i)$
the set of all incoming nodes $j\in V$ of $i$ such that $(j,i)\in E$, and
$\mathcal{O}(i)$ the set of outgoing nodes $j\in V$ of $i$ such that $(i,j)\in
E$.

Figure 1: Sample system topology involving IoT network on the edge

We consider general computation tasks that map input data to result of non-
negligible size, including image/video compression, message encoding/decoding,
model training111Model training usually involves multiple training epochs, in
each the gradient of the loss function for every data point in the dataset is
computed., etc. A task involves the computation of data generated by multiple
data sources and the computation results being forwarded to a given
destination (the requester). Specifically, we assume the computation performed
in the network are of $M$ different types, and denote $S$ the set of all task.
A task is represented by a pair $(d,m)\in S$, where $d\in V$ is the result
destination and $m\in[M]$ is a specified computation type222We denote $[M]$
the set of integers $\\{1,2,\dots,M\\}$ throughout the paper.. Let
$r_{i}(d,m)\geq 0$ be the exogenous input rates of the data corresponding to
task $(d,m)$ at data source node $i\in V$. With this model we enable multi-
source-single-destination tasks, namely multiple nodes $i$ with $r_{i}(d,m)>0$
are allowed for a task $(d,m)$, e.g., a monitor process involving multiple
sensors. The input rate $r_{d}(d,m)$ at destination could also be positive,
incorporating applications with locally provided data and computation
offloading.

$G=(V,E)$ |  Network graph $G$, set of nodes $V$ and links $E$
---|---
$M$ | Total number of computation types
$t=(d,m)$ |  Task $t$, destination $d$ and computation type $m$
$S$ | Set of all tasks
$r_{i}(d,m)$ |  Input data rate of task $(d,m)$ at node $i$
$\mathcal{O}(i)$, $\mathcal{I}(i)$ |  Out-neighbors and in-neighbors of node $i$
$f^{-}_{ij}(d,m)$ |  Data flow of task $(d,m)$ on link $(i,j)$
$f^{+}_{ij}(d,m)$ |  Result flow of task $(d,m)$ on link $(i,j)$
$g_{i}(d,m)$ |  Flow assigned to computation of $(d,m)$ at $i$
$a_{m}$ |  Result size per unit input data of computation $m$
$t^{-}_{i}(d,m)$ |  Total data flow of $(d,m)$ at $i$
$t^{+}_{i}(d,m)$ |  Total result flow of $(d,m)$ at $i$
$\phi^{-}_{ij}(d,m)$ |  Fraction of $t^{-}_{i}(d,m)$ forwarded to node $j$, $j\neq 0$
$\phi^{-}_{i0}(d,m)$ |  Fraction of $t^{-}_{i}(d,m)$ assigned to computation at $i$
$\phi^{+}_{ij}(d,m)$ |  Fraction of $t^{+}_{i}(d,m)$ forwarded to node $j$
$F_{ij}$ |  Total flow on link $(i,j)$
$g_{i}^{m}$ |  Computation input amount of task type $m$ at $i$
$\bf{G}_{i}$ |  Vector of computation input for all task types at $i$
$D_{ij}(F_{ij})$ |  Communication cost (e.g. queueing delay) on $(i,j)$
$C_{i}(\bf{G}_{i})$ |  Computation cost (e.g. CPU time) at node $i$
$D$ | Sum of all communication and computation costs
TABLE I:  Major notations used in network model

###  II-B Routing and computation strategy

The injecting data flows of each task are routed to nodes with computation
resources to be computed. After computation, result flows are generated and
routed to corresponding destination nodes. Since data and computation result
are simultaneously forwarded in the network, we distinguish them by
superscript $-$ and $+$, respectively.

Data and result flows. We consider a hop-by-hop routing scheme: Let
$f_{ij}^{-}(d,m)\geq 0$ denote the data flow of task $(d,m)$ on link $(i,j)$
and $f_{ij}^{+}(d,m)\geq 0$ denote the result flow on the same link. Let
$g_{i}(d,m)\geq 0$ denote the computation flow at node $i$ corresponding to
task $(d,m)$, i.e., the data flow forwarded to the processor of node $i$ for
computation. We consider that the result flow generated by computation for
type $m$ is a non-decreasing function $\gamma_{m}(\cdot)$ of the corresponding
computation flow, given as In particular, we assume that $\gamma_{m}(\cdot)$
is a non-negative weighted sum of a linear function and a sign function. Let
$\mathbbm{1}_{A}=1$ if statement $A$ is true and $0$ otherwise, then

$\gamma_{m}(g_{i}(d,m))=a_{m}g_{i}(d,m)+b\mathbbm{1}_{g_{i}(d,m)>0},$

where $a_{m}\geq 0$ is the ratio of result size versus data, $b\geq 0$ is a
fixed overhead due to task partition or the nature of computation. In this
paper, we focus on case $b=0$. Such assumption fits most of the computation
required in modern applications with usually $a_{m}\leq 1$. Whereas we also
allow $a_{m}>1$, representing special types of computation with result size
larger than the input data size, e.g., video rendering, image super-resolution
or file decompression.333We defer the analysis of case $b\neq 0$ to our future
work. As an expansion, subject to minor modifications, the main mathematical
conclusions of this paper apply for cases where $\gamma_{m}(\cdot)$ is
increasing and convex.

Let $t_{i}^{-}(d,m)$ denote the total data traffic of task $(d,m)$ forwarded
and injected to node $i$, and $t_{i}^{+}(d,m)$ denote the total result traffic
forwarded to and generated at node $i$, given by the following

$\displaystyle t_{i}^{-}(d,m)$
$\displaystyle=\sum\nolimits_{j\in\mathcal{I}(i)}f_{ji}^{-}(d,m)+r_{i}(d,m),$
(1) $\displaystyle t_{i}^{+}(d,m)$
$\displaystyle=\sum\nolimits_{j\in\mathcal{I}(i)}f_{ji}^{+}(d,m)+\gamma_{m}(g_{i}(d,m)).$
(2)

Routing and computation strategy. To describe the computation and forwarding
scheme in a distributed fashion, we assume that each node has two virtual
router, one for data flow and one for result flow, as shown in Fig. 2. Let
$\phi_{ij}^{-}(d,m)$, $\phi_{ij}^{+}(d,m)\in[0,1]$ be the fraction of data or
result flow of task $(d,m)$ at node $i$ forwarded to node $j$. For a coherent
notation, we also use $\phi_{i0}^{-}(d,m)\in[0,1]$ to denote the fraction of
computation flow at node $i$, and use
$\boldsymbol{\phi}^{-}=[\phi_{ij}^{-}(d,m)]_{i\in
V,j\in\\{0\\}\cup\mathcal{O}(i),(d,m)\in S}$,
$\boldsymbol{\phi}^{+}=[\phi_{ij}^{+}(d,m)]_{(i,j)\in E,(d,m)\in S}$ and
$\boldsymbol{\phi}=(\boldsymbol{\phi}^{-},\boldsymbol{\phi}^{+})$ to represent
the system-wide routing and computation strategy. To ensure all tasks are
fulfilled, all injected data must be computed and the result must be delivered
to the destination, given by the following flow conservation: for all
$(d,m)\in S$ and $i\in V$, the data traffic is either computed or forwarded,

$\displaystyle f_{ij}^{-}(d,m)=t_{i}^{-}(d,m)\phi_{ij}^{-}(d,m),$ (3)
$\displaystyle g_{i}(d,m)=t_{i}^{-}(d,m)\phi_{i0}^{-}(d,m),$ (4)
$\displaystyle\sum\nolimits_{j\in\left\\{0\right\\}\bigcup\mathcal{O}(i)}\phi_{ij}^{-}(d,m)=1,$
(5)

and for the result traffic, the destination is a sink:

$\displaystyle f_{ij}^{+}(d,m)=t_{i}^{+}(d,m)\phi_{ij}^{+}(d,m),$ (6)
$\displaystyle\sum_{j\in\mathcal{O}(i)}\phi_{ij}^{+}(d,m)=\begin{cases}1,\quad\text{if
}i\neq d,\\\ 0,\quad\text{if }i=d.\end{cases}$ (7)

Figure 2: Example of data and result flows for nodes $j\to i\to k$

###  II-C Communication and computation cost

For mathematical formulation, we ignore the request messages sent by the
requester to data sources, as the size of such messages are typically
negligible compared to data or result, and could be delivered using a separate
channel.

Instead of sharp bandwidth or computation capacity constraint in [15] or
linear costs in [16], we assign convex costs for D2D communication and local
computation depending on corresponding flow rate, a more general assumption
reflecting network congestion status. In particular, define the total flow on
link $(i,j)\in E$ as

$\textstyle F_{ij}=\sum_{(d,m)\in
S}\left(f_{ij}^{+}(d,m)+f_{ij}^{-}(d,m)\right),$

we assume the communication cost on $(i,j)$ is $D_{ij}(F_{ij})$, where
$D_{ij}(\cdot)$ is increasing, continuously differentiable and convex.
Similarly, denote $\boldsymbol{G}_{i}$ the data rate vector of the amount of
computation for each type performed at node $i$, i.e.,

$\textstyle\boldsymbol{G}_{i}=\left[g_{i}^{1},g_{i}^{2},\cdots,g_{i}^{M}\right],$

where $g_{i}^{m}$ is the total computation amount of type $m$,

$\textstyle g_{i}^{m}=\sum_{d:(d,m)\in S}g_{i}(d,m).$

We denote $C_{i}(\boldsymbol{G}_{i})$ the cost for node $i$ to fulfill the
computation load $\boldsymbol{G}_{i}$, where function $C_{i}(\cdot)$ is an
increasing, continuously differentiable and convex multivariable function,
i.e., increasing on every coordinate and jointly convex in
$g_{i}^{1},\cdots,g_{i}^{M}$.

Such assumption of cost functions incorporates a variety of existing
performance metrics. For example,
$D_{ij}(F_{ij})={F_{ij}}/\left({\mu_{ij}-F_{ij}}\right)$ gives the average
number of packets waiting for or under transmission at link $(i,j)$, provided
that $\mu_{ij}$ is the service rate in an M/M/1 queue model[17] and
$F_{ij}<\mu_{ij}$. One could also approximate the sharp capacity constraint
$F_{ij}\leq C_{ij}$ (e.g., in [15]) by a smooth convex function that goes to
infinity when approaching to capacity limit $C_{ij}$.

Note that we measure the computation cost $C_{i}(\boldsymbol{G}_{i})$ as a
function of data rate, which is adopted in network function virtualization
(NFV) studies[18], but most previous computation offloading researches [5][19]
assign data amount and computation workload separately for a task. Our
formulation can be reduced to the latter if choose certain cost function. For
instance, consider the scenario where users make requests with input data size
$d_{m}$ (bits) and computation workload $c_{m}$ (CPU cycles). By setting
$C_{i}(\boldsymbol{G_{i}})=\sum_{m\in[M]}c_{m}g_{i}^{m}/d_{m}$, we measure CPU
cycles as cost. Or by setting $v_{i}^{m}$ to be the computation speed of type
$m$ at $i$ and
$C_{i}(\boldsymbol{G_{i}})=\sum_{m\in[M]}c_{m}g_{i}^{m}/(d_{m}v_{i}^{m})$, we
measure CPU runtime as cost.

Note that for a network with heterogeneous computation resources, our
formulation is even more flexible than that in [2][19], where there a task
always assigns the same computation workload wherever it is computed. In fact,
our model captures the fact that in practical network edge, the workload for a
certain task may be very different depending on where to perform it, e.g.,
some parallelizable computation is easier at nodes equipping GPU, but slower
at others.

###  II-D Joint routing and computation offloading problem

In this paper, we aim at minimizing the overall cost of edges and devices for
both communication and computation,

$\displaystyle\min_{\boldsymbol{\phi}}\quad$ $\displaystyle D=\sum_{(i,j)\in
E}D_{ij}(F_{ij})+\sum_{i\in V}C_{i}(\boldsymbol{G}_{i})$ (8) such that (1) to
(7) hold

Note that problem (8) is not convex in $\boldsymbol{\phi}$. We will
demonstrate by example that standard gradient-based methods solving for KKT
points highly likely generate sub-optimal solutions.

##  III Sufficient Optimality Condition

In this section, we provide necessary and sufficient conditions for the global
optimum of problem (8). Our analysis follows [20], while we make non-trivial
extensions for considering both data and result flows, as well as in-network
computation. Remind that we focus on the case with $b=0$, where the size of
computation result is in proportion to the data size, i.e.,
$\gamma_{m}(g_{i}(d,m))=a_{m}g_{i}(d,m)$.

We start by giving closed-form expressions to the derivatives of the total
cost $D$. For an increment of exogenous data flow $r_{i}(d,m)$, the increase
of $D$ (i.e. the marginal cost) is caused by two aspects, (1) the cost of
forwarding extra data flow to $i$’s outgoing neighbors $j\in\mathcal{O}(j)$,
and (2) the cost of assigning extra computation load at $i$’s computation
unit. Note that the first aspect could be further decomposed into two terms,
the extra cost on the out-link $(i,j)$ and the extra cost at the next-hop $j$,
and similarly could the second aspect. Thus formally, the marginal cost at
node $i$ is given as

$\displaystyle\frac{\partial D}{\partial r_{i}(d,m)}$
$\displaystyle=\sum_{j\in\mathcal{O}(i)}\phi_{ij}^{-}(d,m)\left[D_{ij}^{\prime}(F_{ij})+\frac{\partial
D}{\partial r_{j}(d,m)}\right]$ (9)
$\displaystyle+\phi_{i0}^{-}(d,m)\left[\frac{\partial
C_{i}(\boldsymbol{G_{i}})}{\partial g_{i}^{m}}+a_{m}\frac{\partial D}{\partial
t_{i}^{+}(d,m)}\right],$

where ${\partial D}/{\partial t_{i}^{+}(d,m)}$ denotes the marginal cost
corresponding to an increment of result traffic at $i$. Similarly to (9), this
marginal of the result traffic is a weighted sum of extra costs at out-links
and at next-hope nodes, given as

$\displaystyle\frac{\partial D}{\partial
t_{i}^{+}(d,m)}=\sum_{j\in\mathcal{O}(i)}\phi_{ij}^{+}(d,m)\left[D_{ij}^{\prime}(F_{ij})+\frac{\partial
D}{\partial t_{j}^{+}(d,m)}\right].$ (10)

Note that (9) and (10) can be calculated recursively, whereas we defer the
detailed mechanism to Section IV. Meanwhile, with an increment of
$\phi^{-}_{ij}$ or $\phi^{+}_{ij}$, the extra cost could also be decomposed in
the similar way as in (9) or (10), respectively.

$\displaystyle\frac{\partial D}{\partial\phi_{ij}^{-}(d,m)}$
$\displaystyle=\begin{cases}t_{i}^{-}(d,m)\left[D_{ij}^{\prime}(F_{ij})+\frac{\partial
D}{\partial r_{j}(d,m)}\right],\text{if }j\neq 0\\\
t_{i}^{-}(d,m)\left[\frac{\partial C_{i}(\boldsymbol{G_{i}})}{\partial
g_{i}^{m}}+a_{m}\frac{\partial D}{\partial t_{i}^{+}(d,m)}\right],\text{ if
}j=0\end{cases}$ (11) $\displaystyle\frac{\partial
D}{\partial\phi_{ij}^{+}(d,m)}$
$\displaystyle=t_{i}^{+}(d,m)\left[D_{ij}^{\prime}(F_{ij})+\frac{\partial
D}{\partial t_{j}^{+}(d,m)}\right].$ (12)

Then we could solve (8) by minimizing the Lagrangian

$\displaystyle L(\boldsymbol{\phi},\lambda^{-},\lambda^{+})=D$ (13)
$\displaystyle-\sum_{i\in V}\sum_{(d,m)\in
S}\lambda_{i,(d,m)}^{-}(\sum_{j\in\left\\{0\right\\}\cup\mathcal{O}(i)}\phi_{ij}^{-}(d,m)-1)$
$\displaystyle-\sum_{i\in V}\sum_{(d,m)\in
S}\lambda_{i,(d,m)}^{+}(\sum_{j\in\mathcal{O}(i)}\phi_{ij}^{+}(d,m)-\mathbbm{1}_{i\neq
d})$

subject to the constraints $\phi_{ij}^{-}(d,m)\geq 0$ and
$\phi_{ij}^{+}(d,m)\geq 0$. By setting the derivative of $L$ to $0$, a KKT
necessary condition of the global minimizer to (8) is given by Lemma 1.

###### Lemma 1

Let $b=0$ and $\boldsymbol{\phi}^{-},\boldsymbol{\phi}^{+}$ be the global
solution that minimizes (8), then for all $i$ and $(d,m)$, and for all
$j\in\left\\{0\right\\}\cup\mathcal{O}(i)$ w.r.t. data flow or
$j\in\mathcal{O}(i)$ w.r.t. result flow,

$\displaystyle\frac{\partial
D}{\partial\phi_{ij}^{-}(d,m)}\begin{cases}=\min\limits_{k\in\left\\{0\right\\}\cup\mathcal{O}(i)}\frac{\partial
D}{\partial\phi_{ik}^{-}(d,m)},\quad\text{if }\phi_{ij}^{-}(d,m)>0,\\\
\geq\min\limits_{k\in\left\\{0\right\\}\cup\mathcal{O}(i)}\frac{\partial
D}{\partial\phi_{ik}^{-}(d,m)},\quad\text{if
}\phi_{ij}^{-}(d,m)=0,\end{cases}$ $\displaystyle\frac{\partial
D}{\partial\phi_{ij}^{+}(d,m)}\begin{cases}=\min\limits_{k\in\mathcal{O}(i)}\frac{\partial
D}{\partial\phi_{ik}^{+}(d,m)},\quad\text{if }\phi_{ij}^{+}(d,m)>0,\\\
\geq\min\limits_{k\in\mathcal{O}(i)}\frac{\partial
D}{\partial\phi_{ik}^{+}(d,m)},\quad\text{if
}\phi_{ij}^{+}(d,m)=0.\end{cases}$

Note that the condition in Lemma 1 is not a sufficient condition for
optimality. A toy example for such non-sufficiency is provided in Fig.3: The
only task is $(4,m)$ with input data only occur at node $1$,the
routing/offloading strategy and marginal costs are shown on figure. It is easy
to verify that condition in Lemma 1 is satisfied. However, by increasing
$\phi^{-}_{24}$ and decreasing $\phi^{-}_{23}$, the input marginal ${\partial
D}/{\partial r_{2}}$ will decrease and thus ${\partial
D}/{\partial\phi^{-}_{12}}$ will decrease. In this case, the objective $D$
could be improved by increasing $\phi^{-}_{12}$ and decreasing
$\phi^{-}_{14}$.

Figure 3: A non-optimal situation satisfying the condition in Lemma 1

The underlying intuition for such non-sufficiency is that the condition in
Lemma 1 automatically holds if $t_{i}^{-}(d,m)=0$ and $t_{i}^{+}(d,m)=0$, no
matter what the routing/computation strategy is. Nevertheless, given
$t_{i}^{-}(d,m)$ and $t_{i}^{+}(d,m)$ exist identically in (11) and (12)
respectively for all $j$, we remove them and devise an augmented condition
specified in Theorem 1, which instead is shown to be sufficient for global
optimality.

###### Theorem 1

Let $b=0$ and $\boldsymbol{\phi}^{-},\boldsymbol{\phi}^{+}$ be feasible to
(8), if the following holds: for all $i$ and $(d,m)$, and for all
$j\in\left\\{0\right\\}\cup\mathcal{O}(i)$ w.r.t data flow or
$j\in\mathcal{O}(i)$ w.r.t. result flow,

$\displaystyle\delta_{ij}^{-}(d,m)\begin{cases}=\min\limits_{k\in\left\\{0\right\\}\cup\mathcal{O}(i)}\delta_{ik}^{-}(d,m),\quad\text{if
}\phi_{ij}^{-}(d,m)>0\\\
\geq\min\limits_{k\in\left\\{0\right\\}\cup\mathcal{O}(i)}\delta_{ik}^{-}(d,m),\quad\text{if
}\phi_{ij}^{-}(d,m)=0\end{cases}$
$\displaystyle\delta_{ij}^{+}(d,m)\begin{cases}=\min\limits_{k\in\mathcal{O}(i)}\delta_{ik}^{+}(d,m),\quad\text{if
}\phi_{ij}^{+}(d,m)>0\\\
\geq\min\limits_{k\in\mathcal{O}(i)}\delta_{ik}^{+}(d,m),\quad\text{if
}\phi_{ij}^{+}(d,m)=0\end{cases}$

where the augmented marginals $\delta_{ij}^{-}(d,m)$ and
$\delta_{ij}^{+}(d,m)$ are

$\displaystyle\delta_{ij}^{-}(d,m)$
$\displaystyle=\begin{cases}D_{ij}^{\prime}(F_{ij})+\frac{\partial D}{\partial
r_{j}(d,m)},&\quad\text{if }j\in\mathcal{O}(i)\\\ \frac{\partial
C_{i}(\boldsymbol{G_{i}})}{\partial g_{i}^{m}}+a_{m}\frac{\partial D}{\partial
t_{i}^{+}(d,m)},&\quad\text{if }j=0\end{cases}$ (14)
$\displaystyle\delta_{ij}^{+}(d,m)$
$\displaystyle=D_{ij}^{\prime}(F_{ij})+\frac{\partial D}{\partial
t_{j}^{+}(d,m)},$

then $(\boldsymbol{\phi}^{-},\boldsymbol{\phi}^{+})$ is a global optimal
solution to (8).

###### Proof:

See Appendix. ∎

Theorem 1 is the main theoretical result of this paper, and in fact a
practical criterion for algorithm implementation. As a simple illustration of
the difference between Theorem 1 and Lemma 1, we further assume the network in
Fig. 3 has linear communication costs. It turns out that for any routing
scheme satisfying Theorem 1, we must have $\phi^{-}_{12}(4,m)=1$ and
$\phi^{-}_{24}(4,m)=1$, which precisely indicates the shortest path $1\to 2\to
4$ for data flow.

##  IV Distributed and Adaptive Algorithm

In this section, we introduce a distributed and adaptive algorithm that
converges to the global optimal solution of (8) specified by Theorem 1, based
on scaled gradient projection. We allow nodes update their routing-computation
strategies in an autonomous and asynchronous manner, and adapt to the changes
of input rate and network topology. Our method follows Xi and Yeh[21], and
further distinguishes data and result flows by improving control messages
exchanging protocol.

###  IV-A Algorithm overview and loop-free property

We first introduce the _loop-free_ property of a global strategy
$\boldsymbol{\phi}$. For a task $(d,m)$, there is a _data path_ from node $i$
to node $j$ ($j\neq i$) if there is a sequence of node $n_{1},\cdots,n_{L}$
such that $(n_{l},n_{l+1})\in E$ and $\phi^{-}_{n_{l}n_{l+1}}(d,m)>0$ for
$l=1,\cdots,L-1$, with $n_{1}=i$ and $n_{L}=j$, where $L$ is the _hop number_
for this data path. We say $\boldsymbol{\phi}$ has a _data loop_ if there
exists task $(d,m)$ and node $i$, $j$ such that $i$ has a data path to $j$,
and vice versa. Similarly, we can define _result path_ and _result loop_ for
result flows. Then, we say strategy $\boldsymbol{\phi}$ is _loop-free_ if it
has neither data loop nor result loop. Loop-free is a fundamental requirement
for $\boldsymbol{\phi}$ to be feasible and is guaranteed throughout the
algorithm. Because given a slight increase in input rate, a loop may cause the
on-loop flow to build up to infinity. 444We allow loops concatenated by a data
path and a result path of the same task, which occurs in scenarios where the
destination is the data source.

We denote respectively by $\boldsymbol{\phi}^{-}_{i}(d,m)$ and
$\boldsymbol{\phi}^{+}_{i}(d,m)$ node $i$’s strategy vector
$\left[\phi^{-}_{ij}(d,m)\right]_{j\in\left\\{0\right\\}\cup\mathcal{O}(i)}$
and $\left[\phi^{+}_{ij}(d,m)\right]_{j\in\mathcal{O}(i)}$, and by
$\boldsymbol{\delta}^{-}_{i}(d,m)$ and $\boldsymbol{\delta}^{+}_{i}(d,m)$ the
vectors of augmented marginals (defined in (14)),
$\left[\delta^{-}_{ij}(d,m)\right]_{j\in\left\\{0\right\\}\cup\mathcal{O}(i)}$
and $\left[\delta^{+}_{ij}(d,m)\right]_{j\in\mathcal{O}(i)}$. We assume the
network starts with a feasible and loop-free state ${\boldsymbol{\phi}}^{0}$.
At $t$-th iteration, each node $i$ updates its strategy corresponding to task
$(d,m)$ with the following scaled gradient projection variant

$\displaystyle\boldsymbol{\phi}^{-}_{i}(d,m)^{t+1}=\underset{\boldsymbol{v}\in\mathcal{D}_{i}^{-}(d,m)^{t}}{\arg\min}\boldsymbol{\delta}^{-}_{i}(d,m)^{t}\cdot(\boldsymbol{v}-\boldsymbol{\phi}^{-}_{i}(d,m)^{t})$
(15)
$\displaystyle\quad+(\boldsymbol{v}-\boldsymbol{\phi}^{-}_{i}(d,m)^{t})^{T}M^{-}_{i}(d,m)^{t}(\boldsymbol{v}-\boldsymbol{\phi}^{-}_{i}(d,m)^{t}),$

and $\boldsymbol{\phi}^{+}_{i}(d,m)^{t+1}$ is updated similarly with “$-$”
replaced by “$+$”, where $M^{-}_{i}(d,m)^{t}$, $M^{+}_{i}(d,m)^{t}$ are
symmetric and positive semi-definite scaling matrices designed to achieve good
convergence properties. $\mathcal{D}_{i}^{-}(d,m)^{t}$ is the feasible set of
$\boldsymbol{\phi}^{-}_{i}(d,m)^{t+1}$ given by
$\boldsymbol{\phi}_{i}\geq\boldsymbol{0}$ and

$\displaystyle\sum\nolimits_{j\in\\{0\\}\cup\mathcal{O}(i)}\phi_{ij}=1;\quad\phi_{ij}=0\text{
for }\forall j\in\mathcal{B}_{i}^{-}(d,m)^{t},$

and feasible set $\mathcal{D}_{i}^{+}(d,m)^{t}$ is defined similarly, where
$\mathcal{B}_{i}^{-}(d,m)^{t}$ and $\mathcal{B}_{i}^{+}(d,m)^{t}\subseteq V$
are the _blocked nodes_ of $i$ relevant to data and result of task $(d,m)$ to
guarantee the feasibility and loop-free property. We next describe in detail
the estimation of marginals, and how to obtain matrices $M^{-}_{i}(d,m)^{t}$,
$M^{+}_{i}(d,m)^{t}$ and sets $\mathcal{B}_{i}^{-}(d,m)^{t}$,
$\mathcal{B}_{i}^{+}(d,m)^{t}$. We emphasis that our proposed algorithm is not
pure gradient-based, as the gradients are replaced by the augmented marginals
$\boldsymbol{\delta}^{-}_{i}(d,m)$ and $\boldsymbol{\delta}^{+}_{i}(d,m)$,
corresponding to Theorem 1.

###  IV-B Marginal cost estimation by broadcast

Each node $i$ needs to compute its augmented marginal cost vectors
$\boldsymbol{\delta}^{-}_{i}(d,m)$ and $\boldsymbol{\delta}^{+}_{i}(d,m)$
following (14). Node $i$ can directly estimate $D^{\prime}_{ij}(F_{ij})$ and
$C^{\prime}_{i}(\boldsymbol{G_{i}})$ while sending and receiving any message
on link $(i,j)$ or performing computation at local computation unit. To obtain
$\frac{\partial D}{\partial r_{i}(d,m)}$ and $\frac{\partial D}{\partial
t^{+}_{i}(d,m)}$, we introduce a two-stage broadcast protocol:

1) To calculate $\frac{\partial D}{\partial t^{+}_{i}(d,m)}$, node $i$ first
waits until receives messages carrying $\frac{\partial D}{\partial
t^{+}_{j}(d,m)}$ from all downstream nodes $j\in\mathcal{O}(i)$ such that
$\phi^{+}_{ij}(d,m)>0$, and estimated $D^{\prime}_{ij}(F_{ij})$. Then
calculates its own $\frac{\partial D}{\partial t^{+}_{i}(d,m)}$ according to
(10) and broadcasts this to all upstream nodes $k\in\mathcal{I}(i)$ such that
$\phi^{+}_{ki}(d,m)>0$.

2) For the $\frac{\partial D}{\partial r_{i}(d,m)}$, a similar procedure is
used to compute and broadcast $\frac{\partial D}{\partial r_{i}(d,m)}$
according to (9). Note that node $i$ must obtain $\frac{\partial D}{\partial
t^{+}_{i}(d,m)}$ before calculating and broadcasting $\frac{\partial
D}{\partial r_{i}(d,m)}$.

With the loop-free property held, such broadcast starting at destination $d$
(with $\frac{\partial D}{\partial t^{+}_{d}(d,m)}=0$) is guaranteed to
traverse throughout the network. Though data flows have no explicit (and
fixed) sink, the loop-free property guarantees the broadcast of stage 2) could
successfully start with the last node of each data path.

###  IV-C Blocked nodes and scaling matrices

To achieve the feasibility and loop-free property, we consider sets
$\mathcal{B}_{i}^{-}(d,m)^{t},\mathcal{B}_{i}^{+}(d,m)^{t}\subseteq V$, to
nodes in which node $i$ is forbidden to forward data or result of task
$(d,m)$, respectively. By Theorem 1 combined with expression (9) and (10), at
a global optimal strategy, the input marginals $\frac{\partial D}{\partial
r_{i}(d,m)}$ or $\frac{\partial D}{\partial t^{+}_{i}(d,m)}$ should be
monotonically decreasing along any data-path or result-path. We thus mandate
that node $i$ should not increase flow rate to a neighbor $j$ that either (1)
has higher input marginal, or (2) could form a data/result-path containing
some link $(p,q)$ and $q$ has higher input marginal than $p$. We denote by
$\mathcal{B}_{i}^{-}(d,m)^{t},\mathcal{B}_{i}^{+}(d,m)^{t}$ the sets of such
neighbor $j$. Note that $\mathcal{B}_{i}^{-}(d,m)^{t}$ and
$\mathcal{B}_{i}^{+}(d,m)^{t}$ are defined separately according to
$\frac{\partial D}{\partial r_{i}(d,m)}$ and $\frac{\partial D}{\partial
t^{+}_{i}(d,m)}$. Then, the loop-free property is maintained throughout the
algorithm if sets of block nodes are practiced in each iteration. The
feasibility is also guaranteed since when current state approaches a
link/processor capacity, the corresponding marginal will grow to infinity,
preventing any flow increase. The readers are referred to [20] for detail.

The scaling matrices $M^{-}_{i}(d,m)^{t}$ and $M^{+}_{i}(d,m)^{t}$ are
introduced to improve the convergence speed while guaranteeing convergence
from arbitrary initial points [21]. Specifically,

$\displaystyle
M_{i}^{+}(d,m)^{t}=\frac{t_{i}^{+}(d,m)^{t}}{2}\times\text{diag}\\{A_{ij}(D^{0})+$
$\displaystyle\left|\mathcal{O}(i)\backslash\mathcal{B}_{i}^{+}(d,m)^{t}\right|h_{j}^{+}(d,m)^{t}A(D^{0})\\}_{j\in\mathcal{O}(i)\backslash\mathcal{B}_{i}^{+}(d,m)^{t}},$

where $D^{0}$ is the overall cost at initial state, $h_{j}^{+}(d,m)^{t}$ is
the maximum hop number among all existing result paths from $j$ to destination
$d$, operator $diag$ forms a diagonal matrix, and

$\displaystyle A_{ij}(D^{0})$
$\displaystyle=\sup_{D<D^{0}}D^{\prime\prime}_{ij}(F_{ij}),\quad A(D^{0})$
$\displaystyle=\max_{(i,j)\in E}A_{ij}(D^{0}).$

The definition of $M_{i}^{-}(d,m)^{t}$ is almost a repetition as above, but in
terms of the data flow.

###  IV-D Asynchronous convergence and complexity

Our algorithm allows nodes to update their variables asynchronously, or with a
non-perfect synchronization due to practical constraints such as the broadcast
delay in a large-scale network. To formulate this asynchrony, we assume that
at $t$-th iteration, only one node $i$ updates either its
$\boldsymbol{\phi}_{i}^{-}(d,m)$ or $\boldsymbol{\phi}_{i}^{+}(d,m)$ for one
task $(d,m)\in S$, and let

$\displaystyle T_{\boldsymbol{\phi}^{-}_{i}(d,m)}$
$\displaystyle=\left\\{t\big{|}\text{ node }i\text{ update its
}\boldsymbol{\phi}_{i}^{-}(d,m)\text{ at iteration }t\right\\},$

and similarly as $T_{\boldsymbol{\phi}^{+}_{i}(d,m)}$, then Theorem 2 holds.

###### Theorem 2

Assume $b=0$ and the network is static, if

$\displaystyle\lim_{t\to\infty}\left|T_{\boldsymbol{\phi}^{-}_{i}(d,m)}\right|=\infty,\quad\lim_{t\to\infty}\left|T_{\boldsymbol{\phi}^{+}_{i}(d,m)}\right|=\infty,$

then the constructed sequence
$\left\\{{\boldsymbol{\phi}}^{t}\right\\}_{t\to\infty}$ converges to a point
${\boldsymbol{\phi}}^{*}$, where ${\boldsymbol{\phi}}^{*}$ is feasible and
loop-free, and the condition in Theorem 1 holds (proof see [21]).

We assume that the variables of all nodes are updated one round every time
slot of duration $T$, and every broadcast message described in Section IV-B is
sent once in a slot. There are $2|E|$ transmissions of control messages
corresponding to a task in one slot, and thus totally $2|S||E|$ transmissions,
with on average $2|S|/T$ per link/second and at most $2\bar{d}|S|$ for each
node, where $\bar{d}$ is the largest out-degree, and a singe broadcast message
has $O(1)$ size. Moreover, let $\bar{h}$ be the maximum path hop, and $t_{c}$
be the maximum time for control message transmission, the broadcast procedure
yields a delay of at most $2\bar{h}t_{c}$. The variable size for individual
node optimization problem is at most $2\bar{d}|S|$, where each problem,
although is a scaled gradient projection, could be efficiently solved by
various commercial solvers since the scaling matrix is PSD and diagonal, and
the constraint set is simplex.

##  V Numerical Evaluation

In this section, we evaluate the scaled gradient projection algorithm, i.e.,
SGP proposed in Section IV by simulation. We implement several baseline
algorithms and compare the performance of those against SGP over different
networks and parameter settings. Note that we set $b=0$ in all experiments.

Network | Parameters
---|---
Topology |  $|V|$ |  $|E|$ |  $|S|$ |  $|\mathcal{R}|$ | Link |  $\bar{d}_{ij}$ | Comp |  $\bar{s}_{i}$
Connected-ER |  $20$ |  $40$ |  $15$ |  $5$ | Queue |  $10$ | Sum-Queue |  $12$
Balanced-tree |  $15$ |  $14$ |  $20$ |  $5$ | Queue |  $20$ | Sum-Queue |  $15$
Fog |  $19$ |  $30$ |  $30$ |  $5$ | Queue |  $20$ | Sum-Queue |  $17$
Abilene |  $11$ |  $14$ |  $10$ |  $3$ | Queue |  $15$ | Sum-Queue |  $10$
LHC |  $16$ |  $31$ |  $30$ |  $5$ | Queue |  $15$ | Sum-Queue |  $15$
GEANT |  $22$ |  $33$ |  $40$ |  $7$ | Queue |  $25$ | Sum-Queue |  $20$
SW |  $100$ |  $320$ |  $120$ |  $10$ | (both) |  $20$ | (both) |  $20$
Other |  $d_{\max}=30$, $d_{\min}=10$, $s_{\max}=30$, $s_{\min}=2$
Parameters |  $M=5$, $r_{\min}=0.5$, $r_{\max}=1.5$
TABLE II: Simulated Network Scenarios

We summarize the simulation scenarios in Table II. We evaluate the algorithms
in the following different network topologies: Connected-ER is a connectivity-
guaranteed Erdős–Rényi graph, generated by uniformly-randomly creating links
with probability $p=0.1$ on a linear network concatenating all nodes.
Balanced-tree is a complete binary tree. Fog is a sample topology for fog-
computing, where nodes on the same layer are linearly linked in a balance tree
[22]. Abilene is the topology of the predecessor of _Internet2 Network_ [23].
GEANT is a pan-European data network for the research and education community
[23]. SW (small-world) is a ring-like graph with additional short-range and
long-range edges [24].

Figure 4: Normalized total cost for network scenarios in Table II

(a) Topology _Connected-ER_ , link width equal to link capacity $d_{ij}$ and
node size equal to computation capacity $s_{i}$

(b) Convergence trajectory of GP and SGP subject to server failure at S1

(c) Total cost versus scaled input rate

(d) $L_{\text{data}}$, $L_{\text{result}}$ and their ratio versus $a_{m}$

Figure 5: Topology and results in scenario _Connected-ER_

Table II also summarizes the number of nodes $|V|$ and edges $|E|$, as well as
the number of tasks $|S|$ in each network. We set $a_{m}$ to be exponential
with mean value $0.5$ and truncated into interval $[0.1,5]$, considering that
most computations have $a_{m}$ smaller that $1$, but special types like video
rendering have relatively larger $a_{m}$. Each task is randomly assigned with
one computation type and one destination node, along with $|\mathcal{R}|$
random active data source (i.e. $r_{i}(d,m)>0$). The input rate $r_{i}(d,m)$
of each active data source is chosen u.a.r. in $[r_{\min},r_{\max}]$. Link is
the type of link cost $D_{ij}(\cdot)$, where _Linear_ denotes a linear link
cost with unit cost $d_{ij}$, i.e. $D_{ij}(F_{ij})=d_{ij}F_{ij}$, and _Queue_
denotes a queueing delay with link capacity $d_{ij}$, i.e.
$D_{ij}(F_{ij})=\frac{F_{ij}}{d_{ij}-F_{ij}}$. Comp is the type of computation
cost $C_{i}(\bf{G}_{i})$, where _Sum-Linear_ denotes a weighted sum of linear
cost for each type, i.e.
$C_{i}(\boldsymbol{G_{i}})=s_{i}\sum_{m}c_{im}g_{i}^{m}$, and _Sum-Queue_
denotes a queueing delay-like computation cost with capacity $s_{i}$, i.e.
$C_{i}(\boldsymbol{G}_{i})=\frac{\sum_{m}c_{im}g_{i}^{m}}{s_{i}-\sum_{m}c_{im}g_{i}^{m}}$,
where the weights $c_{im}$ is u.a.r. drawn from $[1,5]$. The parameters
$d_{ij}$ are u.a.r. drawn from $[0,2\bar{d}_{ij}]$ and truncated into
$[d_{\min},d_{\max}]$. Parameter $s_{i}$ are exponential random variables with
mean $\bar{s}_{i}$ truncated into $[s_{\min},s_{\max}]$ for _Sum-Queue_ , or
uniform with mean $\bar{s}_{i}$ for _Sum-Linear_.

We implement the following baseline algorithms. Since this paper is the first
to study joint routing and computation partial offloading in arbitrary network
topologies with congestion-dependent cost and non-negligible result size, we
make adaptation to these baselines to fit our model.

GP(Gradient Projection): similar to SGP but with the scaling matrices
$M^{-}_{i}(d,m)$ and $M^{+}_{i}(d,m)$ being identity matrix multiplied by
universal stepsize $0.1$. GP and SGP lead to the same global strategy but with
different convergence speed.

SPOO(Shortest Path Optimal Offloading): fixes the routing variables
$\phi^{-}_{ij}$,$\phi^{+}_{ij}$ with $j\neq 0$ to the shortest path (measured
with marginal cost at $F_{ij}=0$, accounting for the propagation delay without
queueing effect), and studies the optimal offloading along these paths.
Similar strategy is considered in [11] with linear-topology and partial
offloading. Note that when destination and data source are the same node, SPOO
is restricted to local computation.

LCOR(Local Computation Optimal Routing): computes at the data sources (or with
minimum offloading if pure local computation is not feasible), optimally route
the result to destinations using scaled gradient projection in [25].

LPR(Linear Program Rounded): the joint path-routing and offloading method by
[6], which does not consider partial offloading, congestible links and result
flow. To adapt LPR’s linear link costs to our schemes, we use the marginal
cost at zero flow. To ensure sufficient communication resources for the result
flow, we assign a saturate-factor of $0.7$ for queueing delay costs, i.e., the
data flow could not exceed $0.7$ times real capacity. Shortest path routing is
used for result flow.

Fig.4 compares the total cost $D$ of different algorithms at the steady state
over networks in Table II (we omit GP as it has the same steady state
performance with SGP), where the bar heights of each scenario are normalized
according to the worst algorithm. We test both linear cost and queueing delay
with other parameters fixed in topology SW, labeled as SW-linear and SW-queue.
Our proposed algorithm SGP significantly outperforms other baselines in all
simulated scenarios, with more than $30\%$ improvement on average over LPR,
which also jointly optimizes routing and task offloading but does not consider
partial offloading and congestible links. The difference of case SW-linear and
SW-queue suggests that our proposed algorithm promises a considerable
improvement to SOTA especially when the networks are congestible. Note that
LCOR and SPOO reflects the optimal objective for routing and offloading
subproblems, respectively. The gain of jointly optimizing over both strategies
could be inferred by comparing SGP against LCOR and SPOO. For example, LCOR
performs extremely bad in topology _Balanced-tree_ , because no routing could
be optimized in a tree topology.

We also perform refined experiments in Connected-ER, with the network topology
and capacity shown in Fig.5(a). There are $4$ major servers as labeled, and we
assume server S1 fails (communication and computation capability disabled,
stop sending data or making requests) at the $100$-th iteration. We compare
the convergence speed of GP and SGP in Fig.5(b) subject to such server
failure. SGP takes much less iterations to converge and adapt to topology
change, showing the advantages of the sophisticatedly designed scaling
matrices.

Fig.5(c) shows the change of total cost subject to universally scaled input
rates $r_{i}(d,m)$, with other parameters fixed. The performance advantage of
SGP has a rapid growth as the network getting more congested, especially
against LPR.

To further illustrate why SGP outperforms baselines significantly with
congestion-dependent cost, we define $L_{\text{data}}$ and $L_{\text{result}}$
as the average travel distance (hop number) of data packages from input to
computation, and that of result packages from generation to being delivered,
respectively.

In Fig. 5(d), we compare $L_{\text{data}}$, $L_{\text{result}}$ for SGP over
different universal $a_{m}$ with other parameters fixed. The trajectories
suggest that the average computation offloading distance grows with $a_{m}$.
Namely, for tasks generating more result with unit input data, SGP tends to
compute them nearer to the destination. Considering that when $a_{m}\gg 1$,
the network is highly congested mainly due to the result flow, thus the
optimal strategy is to offload computation closer to destinations in order to
reduce the result transmission distance and mitigate such congestion, in which
case $L_{\text{data}}$ is large and $L_{\text{result}}$ is small. Note that
the trajectories also imply that the speed of growth of $L_{\text{data}}$ and
descent of $L_{\text{result}}$ are low when $a_{m}$ is small. This is because
when little result flow is generated, the network is lightly congested and the
computation cost dominates the total cost, our algorithm tends to offload
large portion of computation to the servers. However when $a_{m}$ is
sufficiently large, the transmission cost overwhelms computation cost, then
SGP speeds up shifting computation sites closer to destination. The above
behavior demonstrates the underlying optimality of our proposed method, namely
reaching a “balance” among the cost for data forwarding, result forwarding and
computation, and therefore optimizes the overall cost. As a comparison, LPR
does not consider these aspects, the solution hardly changes with $a_{m}$,
implying the cost of forwarding results could grow extremely high with large
$a_{m}$.

##  VI Conclusion

We propose a novel joint routing and computation offloading model
incorporating the result flow, partial offloading and multi-hop routing for
both data and result. This is also the first flow model analysis of
computation offloading adopting congestion-dependent link cost and arbitrary
network topology.We propose a total cost minimization problem to decide
optimal routing-computation strategy. We optimally solve this non-convex
problem by providing necessary and sufficient optimality conditions, and
devise a fully distributed and scalable algorithm that reaches the global
optimal.

## References

  * [1] Ericsson. Ericsson mobility report (2021, Nov.). [Online]. Available: https://www.ericsson.com/en/reports-and-papers/mobility-report
  * [2] Y. Sahni, J. Cao, L. Yang, and Y. Ji, “Multi-hop multi-task partial computation offloading in collaborative edge computing,” _IEEE Transactions on Parallel and Distributed Systems_ , vol. 32, no. 5, pp. 1133–1145, 2020.
  * [3] Y. Sahni, J. Cao, S. Zhang, and L. Yang, “Edge mesh: A new paradigm to enable distributed intelligence in internet of things,” _IEEE access_ , vol. 5, pp. 16 441–16 458, 2017.
  * [4] K. Zhu, W. Zhi, X. Chen, and L. Zhang, “Socially motivated data caching in ultra-dense small cell networks,” _IEEE Network_ , vol. 31, no. 4, pp. 42–48, 2017.
  * [5] Z. Hong, W. Chen, H. Huang, S. Guo, and Z. Zheng, “Multi-hop cooperative computation offloading for industrial iot–edge–cloud computing environments,” _IEEE Transactions on Parallel and Distributed Systems_ , vol. 30, no. 12, pp. 2759–2774, 2019.
  * [6] B. Liu, Y. Cao, Y. Zhang, and T. Jiang, “A distributed framework for task offloading in edge computing networks of arbitrary topology,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 4, pp. 2855–2867, 2020\.
  * [7] H. Al-Shatri, S. Müller, and A. Klein, “Distributed algorithm for energy efficient multi-hop computation offloading,” in _2016 IEEE International Conference on Communications (ICC)_. IEEE, 2016, pp. 1–6.
  * [8] Y. Sahni, J. Cao, and L. Yang, “Data-aware task allocation for achieving low latency in collaborative edge computing,” _IEEE Internet of Things Journal_ , vol. 6, no. 2, pp. 3512–3524, 2018.
  * [9] Q. Luo, W. Shi, and P. Fan, “Qoe-driven computation offloading: Performance analysis and adaptive method,” in _2021 13th International Conference on Wireless Communications and Signal Processing (WCSP)_. IEEE, 2021, pp. 1–5.
  * [10] W. Shi, J. Zhang, R. Zhang, and K. Hu, “An area-based offloading policy for computing offloading in mec-assisted wireless mesh network,” in _2019 IEEE/CIC International Conference on Communications in China (ICCC)_. IEEE, 2019, pp. 507–511.
  * [11] X. He, R. Jin, and H. Dai, “Multi-hop task offloading with on-the-fly computation for multi-uav remote edge computing,” _IEEE Transactions on Communications_ , 2021.
  * [12] C. Funai, C. Tapparello, and W. Heinzelman, “Computational offloading for energy constrained devices in multi-hop cooperative networks,” _IEEE Transactions on Mobile Computing_ , vol. 19, no. 1, pp. 60–73, 2019.
  * [13] B. McMahan, E. Moore, D. Ramage, S. Hampson, and B. A. y Arcas, “Communication-efficient learning of deep networks from decentralized data,” in _Artificial intelligence and statistics_. PMLR, 2017, pp. 1273–1282.
  * [14] Z. Hong, H. Huang, S. Guo, W. Chen, and Z. Zheng, “Qos-aware cooperative computation offloading for robot swarms in cloud robotics,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 4, pp. 4027–4041, 2019.
  * [15] B. Liu, K. Poularakis, L. Tassiulas, and T. Jiang, “Joint caching and routing in congestible networks of arbitrary topology,” _IEEE Internet of Things Journal_ , vol. 6, no. 6, pp. 10 105–10 118, 2019.
  * [16] J. Ren, G. Yu, Y. He, and G. Y. Li, “Collaborative cloud and edge computing for latency minimization,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 5, pp. 5031–5044, 2019.
  * [17] D. Bertsekas and R. Gallager, _Data networks_. Athena Scientific, 2021.
  * [18] J. Zhang, A. Sinha, J. Llorca, A. Tulino, and E. Modiano, “Optimal control of distributed computing networks with mixed-cast traffic flows,” in _IEEE INFOCOM 2018-IEEE Conference on Computer Communications_. IEEE, 2018, pp. 1880–1888.
  * [19] Z. Chen, Q. Ma, L. Gao, and X. Chen, “Edgeconomics: Price competition and selfish computation offloading in multi-server edge computing networks,” in _2021 19th International Symposium on Modeling and Optimization in Mobile, Ad hoc, and Wireless Networks (WiOpt)_. IEEE, 2021, pp. 1–8.
  * [20] R. Gallager, “A minimum delay routing algorithm using distributed computation,” _IEEE transactions on communications_ , vol. 25, no. 1, pp. 73–85, 1977.
  * [21] Y. Xi and E. M. Yeh, “Node-based optimal power control, routing, and congestion control in wireless networks,” _IEEE Transactions on Information Theory_ , vol. 54, no. 9, pp. 4081–4106, 2008.
  * [22] K. Kamran, E. Yeh, and Q. Ma, “Deco: Joint computation, caching and forwarding in data-centric computing networks,” in _Proceedings of the Twentieth ACM International Symposium on Mobile Ad Hoc Networking and Computing_ , 2019, pp. 111–120.
  * [23] D. Rossi and G. Rossini, “Caching performance of content centric networks under multi-path routing (and more),” _Relatório técnico, Telecom ParisTech_ , vol. 2011, pp. 1–6, 2011.
  * [24] J. Kleinberg, “The small-world phenomenon: An algorithmic perspective,” in _Proceedings of the thirty-second annual ACM symposium on Theory of computing_ , 2000, pp. 163–170.
  * [25] D. Bertsekas, E. Gafni, and R. Gallager, “Second derivative algorithms for minimum delay distributed routing in networks,” _IEEE Transactions on Communications_ , vol. 32, no. 8, pp. 911–919, 1984.

## Appendix

### Proof of Theorem 1

For simplicity, we consider the non-destination nodes in this proof, namely we
assume $\sum_{j}\phi_{ij}^{+}=1$ for all $i$, while the derivation is
applicable to destination nodes. We have

$\displaystyle\frac{\partial D}{\partial r_{i}(d,m)}$
$\displaystyle=\sum_{j\in\left\\{0\right\\}\cup\mathcal{O}(i)}\phi_{ij}^{-}(d,m)\delta_{ij}^{-}(d,m)$
$\displaystyle=\sum_{j:\phi_{ij}>0}\phi_{ij}^{-}(d,m)\lambda_{idm}^{-}$
$\displaystyle=\lambda_{idm}^{-},$

and thus

$\displaystyle\delta_{ij}^{-}(d,m)\geq\frac{\partial D}{\partial
r_{i}(d,m)},\forall i\in
V,j\in\left\\{0\right\\}\cup\mathcal{O}(i),\forall(d,m)\in S.$ (16)

Similarly we have

$\displaystyle\delta_{ij}^{+}(d,m)\geq\frac{\partial D}{\partial
t_{i}^{+}(d,m)},\quad\forall i\in V,j\in\mathcal{O}(i),\forall(d,m)\in S.$
(17)

To prove $\phi=(\boldsymbol{\phi}^{-},\boldsymbol{\phi}^{+})$ minimizes
$D=\sum_{(i,j)\in E}D_{ij}(F_{ij})+\sum_{i\in V}C_{i}(\boldsymbol{G_{i}})$,
let $\phi^{*}=(\phi^{-*},\phi^{+*})\neq\phi$ be another set of variable, and
with corresponding forwarding and computation flows
$F_{ij}^{*},\forall(i,j)\in E$ and $\boldsymbol{G_{i}}^{*},\forall i\in V$.
Given both $\phi$ and $\phi^{*}$ are valid forwarding and computing scheme, we
know $(F_{ij},\boldsymbol{G_{i}})$ and $(F_{ij}^{*},\boldsymbol{G_{i}}^{*})$
are in the feasible set of the flow model problem (18), which is a convex
polytope.

$\displaystyle\min_{f^{-},f^{+},g}\quad$ $\displaystyle\sum_{(i,j)\in
E}D_{ij}(F_{ij})+\sum_{i\in V}C_{i}(\boldsymbol{G_{i}})$ (18) such that (3) to
(7) hold, $\displaystyle g_{i}(d,m)\geq 0,\quad\forall i\in V,k\leq M,(d,m)\in
S$ $\displaystyle\begin{cases}f_{ij}^{-}(d,m)\geq 0,\\\ f_{ij}^{+}(d,m)\geq
0,\end{cases}\,\forall(i,j)\in E,(d,m)\in S$

Due to the convexity of the feasible set, for any $\mu\in[0,1]$,
$\left((1-\mu)F_{ij}+\mu
F_{ij}^{*},(1-\mu)\boldsymbol{G_{i}}+\mu\boldsymbol{G_{i}}^{*}\right)$ is also
feasible for (18), we then let

$\displaystyle D(\mu)$ $\displaystyle=\sum_{(i,j)\in
E}D_{ij}((1-\mu)F_{ij}+\mu F_{ij}^{*})$ $\displaystyle+\sum_{i\in
V}C_{i}((1-\mu)\boldsymbol{G_{i}}+\mu\boldsymbol{G_{i}}^{*}).$

Since $D$ is convex in $F_{ij}$ and $\boldsymbol{G}_{i}$, we know $D(\mu)$ is
convex in $\mu$. Thus combining with the arbitrary choice of $\phi^{*}$, the
sufficiency in Theorem 1 is proved if $\frac{dD(\mu)}{d\mu}$ is non-negative
at $\mu=0$. That is, we will show the following is non-negative

$\displaystyle\frac{dD(\mu)}{d\mu}\bigg{|}_{\mu=0}$
$\displaystyle=\sum_{(i,j)\in E}D^{\prime}_{ij}(F_{ij})(F_{ij}^{*}-F_{ij})$
(19) $\displaystyle+\sum_{i\in V}\sum_{m\leq M}\frac{\partial
C_{i}(\boldsymbol{G_{i}})}{\partial g_{i}^{m}}(g_{i}^{m*}-g_{i}^{m})$

Starting with the data flow, multiply both side of (16) by
$\phi_{ij}^{-*}(d,m)$ and sum over $j\in\left\\{0\right\\}\cup\mathcal{O}(i)$,
we have

$\displaystyle\frac{\partial C_{i}(\boldsymbol{G_{i}})}{\partial
g_{i}^{m}}\phi_{i0}^{-*}(d,m)+\sum_{j\in\mathcal{O}(i)}D_{ij}^{\prime}(F_{ij})\phi_{ij}^{-*}(d,m)$
(20) $\displaystyle\geq$ $\displaystyle\frac{\partial D}{\partial
r_{i}(d,m)}-a_{m}\frac{\partial D}{\partial
t_{i}^{+}(d,m)}\phi_{i0}^{-*}(d,m)$
$\displaystyle-\sum_{j\in\mathcal{O}(i)}\frac{\partial D}{\partial
r_{j}(d,m)}\phi_{ij}^{-*}(d,m),$

then multiply both side by
$t_{i}^{-*}(d,m)=\sum_{j\in\mathcal{I}(i)}f_{ji}^{-*}(d,m)+r_{i}(d,m)$, we
have

$\displaystyle\frac{\partial C_{i}(\boldsymbol{G_{i}})}{\partial
g_{i}^{m}}g_{i}^{*}(d,m)+\sum_{j\in\mathcal{O}(i)}D_{ij}^{\prime}(F_{ij})f_{ij}^{-*}(d,m)$
$\displaystyle\geq$ $\displaystyle t_{i}^{-*}(d,m)\frac{\partial D}{\partial
r_{i}(d,m)}-a_{m}\frac{\partial D}{\partial
t_{i}^{+}(d,m)}t_{i}^{-*}(d,m)\phi_{i0}^{-*}(d,m)$ $\displaystyle-$
$\displaystyle\sum_{j\in\mathcal{O}(i)}\frac{\partial D}{\partial
r_{j}(d,m)}t_{i}^{-*}(d,m)\phi_{ij}^{-*}(d,m),$

further sum over $(d,m)\in S$ and $j\in V$, we get

$\displaystyle\sum_{i\in V}\sum_{m\leq M}\frac{\partial
C_{i}(\boldsymbol{G_{i}})}{\partial g_{i}^{m}}g_{i}^{m*}+\sum_{(i,j)\in
E}D_{ij}^{\prime}(F_{ij})F_{ij}^{-*}$ (21) $\displaystyle\geq\sum_{i\in
V}\sum_{(d,m)\in S}t_{i}^{-*}(d,m)\frac{\partial D}{\partial r_{i}(d,m)}$
$\displaystyle-\sum_{i\in V}\sum_{(d,m)\in S}a_{m}\frac{\partial D}{\partial
t_{i}^{+}(d,m)}t_{i}^{-*}(d,m)\phi_{i0}^{-*}(d,m)$
$\displaystyle-\sum_{(d,m)\in S}\sum_{i\in
V}\sum_{j\in\mathcal{O}(i)}\frac{\partial D}{\partial
r_{j}(d,m)}t_{i}^{-*}(d,m)\phi_{ij}^{-*}(d,m),$

where $F_{ij}^{-*}=\sum_{(d,m)\in S}f_{ij}^{-*}(d,m)$.

Meanwhile, by the flow conservation (3) to (7), we know that for all $j\in
V,(d,m)\in S$,

$\displaystyle\sum_{i\in\mathcal{I}(j)}t_{i}^{-*}(d,m)\phi_{ij}^{-*}(d,m)=t_{j}^{-*}(d,m)-r_{j}(d,m).$

Substitute above into the very last term in (21) and cancel, we get

$\displaystyle\sum_{i\in V}\sum_{m\leq M}\frac{\partial
C_{i}(\boldsymbol{G_{i}})}{\partial g_{i}^{m}}g_{i}^{m*}+\sum_{(i,j)\in
E}D_{ij}^{\prime}(F_{ij})F_{ij}^{-*}$ (22) $\displaystyle\geq\sum_{i\in
V}\sum_{(d,m)\in S}r_{i}(d,m)\frac{\partial D}{\partial r_{i}(d,m)}$
$\displaystyle-\sum_{i\in V}\sum_{(d,m)\in S}a_{m}g_{i}^{m*}\frac{\partial
D}{\partial t_{i}^{+}(d,m)}.$

Next, about the flow of computation result, multiply both side of (17) by
$\phi_{ij}^{+*}(d,m)$ and sum over $j\in\mathcal{O}(i)$, we have

$\displaystyle\sum_{j\in\mathcal{O}(i)}D_{ij}^{\prime}(F_{ij})\phi_{ij}^{+*}(d,m)$
(23) $\displaystyle\geq$ $\displaystyle\frac{\partial D}{\partial
t_{i}^{+}(d,m)}-\sum_{j\in\mathcal{O}(i)}\frac{\partial D}{\partial
t_{j}^{+}(d,m)}\phi_{ij}^{+*}(d,m).$

Multiply both side by
$t_{i}^{+*}(d,m)=\sum_{j\in\mathcal{I}(i)}f_{ji}^{+*}(d,m)+a_{m}g_{i}^{*}(d,m)$,
sum over $(d,m)\in S$ and $j\in V$, we get

$\displaystyle\sum_{(i,j)\in E}D_{ij}^{\prime}(F_{ij})F_{ij}^{+*}$ (24)
$\displaystyle\geq\sum_{i\in V}\sum_{(d,m)\in S}t_{i}^{+*}(d,m)\frac{\partial
D}{\partial t_{i}^{+}(d,m)}$ $\displaystyle-\sum_{(d,m)\in S}\sum_{i\in
V}\sum_{j\in\mathcal{O}(i)}t_{i}^{+*}(d,m)\frac{\partial D}{\partial
t_{j}^{+}(d,m)}\phi_{ij}^{+*}(d,m).$

By (3) to (7), we have for all $j\in V,(d,m)\in S$,

$\displaystyle\sum_{i\in\mathcal{I}(j)}t_{i}^{+*}(d,m)\phi_{ij}^{+*}(d,m)=t_{j}^{+*}(d,m)-a_{m}g_{j}^{*}(d,m).$

Substituting above into the very last term in (24) and canceling, we get

$\displaystyle\sum_{(i,j)\in
E}D_{ij}^{\prime}(F_{ij})F_{ij}^{+*}\geq\sum_{i\in V}\sum_{(d,m)\in
S}a_{m}g_{i}^{m*}\frac{\partial D}{\partial t_{i}^{+}(d,m)}.$ (25)

Summing up both side of (22) and (25), we have

$\displaystyle\sum_{i\in V}\sum_{m\leq M}\frac{\partial
C_{i}(\boldsymbol{G_{i}})}{\partial g_{i}^{m}}g_{i}^{m*}+\sum_{(i,j)\in
E}D_{ij}^{\prime}(F_{ij})F_{ij}^{*}$ (26) $\displaystyle\geq\sum_{i\in
V}\sum_{(d,m)\in S}\frac{\partial D}{\partial r_{i}(d,m)}r_{i}(d,m).$

Note that the equality would always hold in (20) and (23) if we substitute
$\phi^{*}$ with $\phi$ in the above reasoning, as a consequence of (9) and
(10). Thus we have the following analogue of (26),

$\displaystyle\sum_{i\in V}\sum_{m\leq M}\frac{\partial
C_{i}(\boldsymbol{G_{i}})}{\partial g_{i}^{m}}g_{i}^{m}+\sum_{(i,j)\in
E}D_{ij}^{\prime}(F_{ij})F_{ij}$ (27) $\displaystyle=\sum_{i\in
V}\sum_{(d,m)\in S}\frac{\partial D}{\partial r_{i}(d,m)}r_{i}(d,m).$

Abstracting (27) from (26), we show (19) and complete the proof.


# Probing the Consistency of Cosmological Contours for Supernova Cosmology

P. Armstrong  https://orcid.org/0000-0003-1997-3649 Mt Stromlo Observatory,
The Research School of Astronomy and Astrophysics, Australian National
University, ACT 2601, Australia  [ H. Qu
https://orcid.org/0000-0003-1899-9791 Department of Physics and Astronomy,
University of Pennsylvania, Philadelphia, PA 19104, USA  D. Brout
https://orcid.org/0000-0001-5201-8374 Department of Astronomy, Boston
University, 725 Commonwealth Ave., Boston, MA 02215, USA  T. M. Davis
https://orcid.org/0000-0002-4213-8783 School of Mathematics and Physics, The
University of Queensland, Brisbane, QLD 4072, Australia  R. Kessler
https://orcid.org/0000-0003-3221-0419 Kavli Institute for Cosmological
Physics, University of Chicago, Chicago, IL 60637, USA  A. G. Kim
https://orcid.org/0000-0001-6315-8743 Physics Division, Lawrence Berkeley
National Laboratory, Berkeley, CA 94720, USA  C. Lidman
https://orcid.org/0000-0003-1731-0497 Mt Stromlo Observatory, The Research
School of Astronomy and Astrophysics, Australian National University, ACT
2601, Australia  M. Sako  https://orcid.org/0000-0003-2764-7093 Department of
Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA
B. E. Tucker  https://orcid.org/0000-0002-4283-5159 Mt Stromlo Observatory,
The Research School of Astronomy and Astrophysics, Australian National
University, ACT 2601, Australia

###### Abstract

As the scale of cosmological surveys increases, so does the complexity in the
analyses. This complexity can often make it difficult to derive the underlying
principles, necessitating statistically rigorous testing to ensure the results
of an analysis are consistent and reasonable. This is particularly important
in multi-probe cosmological analyses like those used in the Dark Energy Survey
and the upcoming Legacy Survey of Space and Time, where accurate uncertainties
are vital. In this paper, we present a statistically rigorous method to test
the consistency of contours produced in these analyses, and apply this method
to the Pippin cosmological pipeline used for Type Ia supernova cosmology with
the Dark Energy Survey. We make use of the Neyman construction, a frequentist
methodology that leverages extensive simulations to calculate confidence
intervals, to perform this consistency check. A true Neyman construction is
too computationally expensive for supernova cosmology, so we develop a method
for approximating a Neyman construction with far fewer simulations. We find
that for a simulated data-set, the 68% contour reported by the Pippin pipeline
and the 68% confidence region produced by our approximate Neyman construction
differ by less than a percent near the input cosmology, however show more
significant differences far from the input cosmology, with a maximal
difference of 0.05 in $\Omega_{M}$, and 0.07 in $w$. This divergence is most
impactful for analyses of cosmological tensions, but its impact is mitigated
when combining supernovae with other cross-cutting cosmological probes, such
as the Cosmic Microwave Background.

###### keywords:

cosmological parameters, type Ia supernovae, astrostatistics

P<EMAIL_ADDRESS>Department of Astronomy and
Astrophysics, University of Chicago, Chicago, IL 60637, USA Centre for
Gravitational Astrophysics, College of Science, The Australian National
University, ACT 2601, Australia National Centre for the Public Awareness of
Science, Australian National University, Canberra, ACT 2601, Australia The ARC
Centre of Excellence for All-Sky Astrophysics in 3 Dimensions (ASTRO 3D)

##  1 Introduction

For much of the history of supernova cosmology, parameter estimation was
performed via Bayesian methods which maximise the likelihood by minimising a
$\chi^{2}$ function between the observed distance modulus of type Ia
supernovae (SNe Ia) and the distance modulus predicted by cosmological theory.

This method of parameter estimation is shown in Equation 4 of
RiessFilippenko1998 and Equation 4 of PerlmutterAldering1999, which detail the
discovery of the accelerated expansion of the Universe using type Ia
supernovae (SNe Ia).

Since these early efforts, parameter estimation has expanded in complexity to
account for additional systematic uncertainties ConleyGuy2011, and to leverage
large simulated data-sets to correct for contamination of core collapse
supernovae KunzHlozek2012, and observational biases KesslerScolnic2017. This
increased complexity is facilitated by cosmological pipelines to perform
accurate parameter estimation. Due to the complex structure of modern
cosmological pipelines, it is no longer possible to analytically define a
likelihood function which describes the analysis being performed by these
pipelines. As such, it is difficult to rigorously test the final cosmological
contours produced by these cosmological pipelines, and to validate that the
reported uncertainties are accurate.

There have been a number of attempts at developing alternative Bayesian
frameworks which do not suffer from a non-analytic likelihood function. One
such example is Approximate Bayesian Computation (JenningsWolf2016;
JenningsMadigan2017, ABC;), which uses realistic simulations to perform
likelihood free parameter inference, at the cost of dramatically increased
computation time due to the large number of simulations required. Another
alternative Bayesian framework is Bayesian hierarchical models (BHM), which
was implemented for supernova cosmology in (HintonDavis2019, Steve;),
(RubinAldering2015, UNITY;), and (MandelThorp2021, BayeSN;). BHMs utilise
multiple layers of connected parameters, allowing for a more complex
analytical likelihood function to be defined.

Though these alternative frameworks have significant advantages over the
$\chi^{2}$ minimisation methods, there has not been wide-spread adoption of
these techniques, and many modern cosmological analyses, such as the Dark
Energy Survey (DarkEnergySurveyCollaborationAbbott2016, DES;), PANTHEON+
BroutScolnic2022, and simulations of the upcoming Legacy Survey of Space and
Time (LSST; LSSTScienceCollaborationAbell2009; SanchezKessler2022;
MitraKessler2022) still use the simpler $\chi^{2}$ method.

As modern analyses still use cosmological pipelines which rely on the
$\chi^{2}$ methodology, it is important to rigorously test these cosmological
pipelines. While each individual component of the SN Ia analysis pipeline is
well-tested LaskerKessler2019; KesslerNarayan2019; PopovicBrout2021;
ToyWiseman2023; TaylorJones2023; KelseySullivan2023; VincenziSullivan2023, a
complete end-to-end consistency check is still necessary to understand the
effects and assumptions that propagate between each individual step of the
pipeline, and account for any systematic issues that may arise.

In this paper, we present a new methodology to validate the reported
cosmological contours of current pipelines. For this effort, we utilize Pippin
HintonBrout2020, which automates a number of key components of the SuperNova
ANAlysis framework (KesslerBernstein2009, SNANA;) used in DES, LSST’s Dark
Energy Survey Collaboration, and PANTHEON+. Pippin and SNANA provide
substantial functionality, including simulations, light-curve fitting,
photometric classification training and evaluation, SNe Ia standardization and
bias corrections, and cosmological fitting.

Previous efforts to construct a confidence region include those published in
BroutScolnic2019, who used 200 simulated samples to demonstrate that the
distribution of best fitting cosmologies produced by their pipeline is
consistent with the average of many cosmological contours, which they took as
an estimate of the confidence region (CR). This estimate of the CR, while
informative, is not statistically rigorous.

To rigorously estimate the confidence region, we make use of the Neyman
construction Neyman1937, a Frequentist methodology that leverages simulations
to produce a confidence region. The Neyman construction does not assume the CR
is Gaussian or elliptical, and is robust to small sample sizes. We compare the
Frequentist confidence region produced from the Neyman construction with the
Bayesian contour produced by our cosmological pipeline in order to test for
consistency. Producing a CR and comparing it to cosmological contours is a
powerful, rigorous, and independent method of evaluating the output of a
cosmological pipeline.

In Section 2 we describe the formalism employed by the Pippin cosmological
pipeline. Section 3 describes our Neyman construction methodology. Finally,
the results of applying our methodology to the Pippin cosmological pipeline
are presented in Section 4.

##  2 Simulation & Analysis with Pippin

Here we briefly describe the simulation of realistic data-sets of
spectroscopically confirmed SN Ia, and the analysis procedure used to produce
a cosmological contour. Throughout this analysis we make use of SNANA for
simulation and analysis, integrated into the Pippin pipeline.

###  2.1 Simulating a Supernova Data-set

We use the SALT2 GuyAstier2007 framework within SNANA for simulating SNe Ia.
This framework models type Ia supernovae with five parameters: redshift, day
of peak rest-frame $B$-band brightness, stretch, color, and apparent peak
brightness in rest-frame $B$-band. SNANA produces simulated observed fluxes by
randomly selecting these model parameters from associated probability
distributions. SNANA also applies host-galaxy extinction, k-correction, and
galactic extinction. For this analysis we use the SALT2 model produced by
TaylorLidman2021, which was trained on a sample of 420 SNe Ia spanning a
redshift range of $\sim$0.1 to $\sim$0.9 with improved zero-point calibration
offsets and Milky Way extinction compared to previous SALT2 models.

####  2.1.1 Bias Correction Simulations

In addition to data-like simulations, we also simulate much larger data-sets
to correct observational biases. As part of this analysis, we investigate the
impact of the cosmology on these bias corrections. Our principal analysis uses
only a single bias correction simulation with the input cosmology set to our
nominal input cosmology ($\Omega_{M}=0.3$, $w=-1.0$), but we repeat our
analysis using many bias correction simulations, with input cosmologies equal
to the input cosmology of the data-set they are correcting, to see if this
affects the Neyman construction.

###  2.2 Analysis

The supernovae in each simulated data-set are fit to determine the SALT2
parameters: amplitude ($x_{0}$), stretch ($x_{1}$), and colour ($c$). From
here, the distance modulus of each SN Ia can be computed via the Tripp
equation Tripp1998

$\mu=m_{B}+\alpha{}x_{1}-\beta{}c+\mathcal{M}-\Delta\mu_{bias}$ (1)

Here $\alpha$ and $\beta$ are global stretch and colour nuisance parameters,
$\mathcal{M}$ is a global offset, and $\Delta\mu_{bias}=\mu-\mu_{true}$ is a
distance bias correction, where $\mu_{true}$ is the true distance modulus.

Pippin makes use of the BEAMS with Bias Correction (KesslerScolnic2017, BBC;)
framework to produce a Hubble Diagram (HD) that has been corrected for both
selection effects and contamination. BBC uses the detailed simulations
described in Section 2.1.1 alongside the BEAMS KunzHlozek2012 method to
correct for both distance biases, contamination, and selection effects
KesslerBrout2019. It then uses a cosmology-independent method
(MarrinerBernstein2011, SALT2mu;) to fit for global nuisance parameters and
standardise the SNe Ia magnitudes.

In order to fit $\alpha$, $\beta$, and $\mathcal{M}$, BBC adopts the
likelihood $\mathcal{L}=\prod_{i=1}^{N}\mathcal{L}_{i}$ where

$\displaystyle\mathcal{L}_{i}$
$\displaystyle=P_{Ia,i}D_{Ia,i}+(1-P_{Ia,i})D_{CC,i}$ (2)

Here $P_{Ia,i}$ is the photometric classification probability for the $i$th
supernova to be an SN Ia. This is usually calculated via a photometric
classifier such as SuperNNova MollerdeBoissiere2020, or Scone QuSako2021,
however for our analysis we do not simulate contamination, so $P_{Ia,i}=1$.
$D_{Ia,i}$ encodes the influence of SNe Ia on the likelihood, including
corrections for observational biases in the data-set. Details of $D_{Ia,i}$
are presented in KesslerScolnic2017. The $D_{CC,i}$ component encodes the
effects of contaminants; however, since our simulations are contaminant free,
it is unimportant for this analysis.

The end result of the BBC framework is a redshift-binned HD. BBC can also
provide an unbinned HD which BroutHinton2021 shows can result in smaller
systematic uncertainties, but is more computationally expensive. Since our
analysis only includes statistical uncertainties, we gain no benefit from
using an unbinned HD, therefore we only use the default, binned HD.

This binned HD is passed to a cosmological fitter to produce the final
cosmological contours. In this analysis we make use of WFit, which measures a
$\chi^{2}$ likelihood over a grid within the parameter space. WFit has the
advantage of being much faster than other methods, although is only suitable
for simple test cases such as the one used in this paper, and not necessarily
suitable for final survey cosmological analysis. We allow the parameter space
of $\Omega_{M}$ to vary below 0, something which is not usual for cosmological
analyses, as our analysis requires WFit to explore large sections of the
$\Omega_{M}$, $w$ parameter space, and we do not wish to artificially truncate
the likelihood surface we produce.

###  2.3 Producing an experiment data-set

Our methodology can be used to validate the contour produced by any
cosmological pipeline, and is not dependent on the details of the data-set
investigated by the cosmological pipeline. As such, we test our methodology on
a simple, simulated dataset which mimics the 3 year DES data-set
BroutSako2019, including the cadence, spectroscopic selection, and
observational noise DAndreaSmith2018 of this data-set. We assume a flat, cold
dark matter ($wCDM$) cosmology with $H_{0}=70$km/s/Mpc, $\Omega_{M}=0.3$, and
$w=-1.0$. The DES 3 year data-set includes a previously released low-$z$
sample from several sources. We simplify the low-$z$ simulation by generating
a DES-like sample for $0.0\leq{}z\leq{}0.08$ with the same statistics as the
low-$z$ sample. Additionally, we only simulate SNe Ia and do not consider
contamination from core-collapse supernovae, so that we can keep our analysis
as simple as possible. The true DES data-set includes data from a variety of
telescopes, as well as misclassified core-collapse SNe, so if our methodology
were to be used to test the DES analysis, these details will need to be
included in all simulations. The redshift distribution of our DES and low-$z$
simulated sample is presented in Figure 1. An example of a simulated
lightcurve is presented in Figure 2.

We analyse this simulated data-set with Pippin in order to produce the
cosmological contour that we aim to validate (shown in Figure 3), and to
calculate the best fitting cosmology:

$\displaystyle\begin{split}\Omega_{M}^{best}&=0.32^{+0.054}_{-0.075}\\\
w^{best}&=-1.00\pm 0.16\end{split}$ (3)

Figure 1: The redshift distribution of the simulated DES sample and simulated
low-$z$ sample. Figure 2: An example of a simulated lightcurve in our
simulated sample which lies at $z=0.4$. Figure 3: The cosmological contour
produced by Pippin for our simulated data-set. The aim of our methodology is
to test the consistency of this contour. The central panel shows the 2-D 68%
and 95% contours, whilst the top and right panel show the marginalised, 1-D
contour for $\Omega_{M}$ and $w$ respectively. Here
$\Omega_{M}^{best}=0.320^{+0.054}_{-0.075}$ and $w^{best}=-1.00\pm 0.16$

##  3 Methods

Here we describe our methodology for estimating the confidence region of a
cosmological data-set using Neyman construction. To aid the reader, a glossary
of terms used throughout is provided in Table 1.

Parameter |  Description
---|---
$\Omega_{M}^{best}$, $w^{best}$ |  The best-fitting output cosmology for our experiment data-set. Described in Section 2.3.
$\Omega_{M}^{\prime}$, $w^{\prime}$ |  A strategically selected input cosmology from which 150 realisations will be drawn. Described in Section 3.1.
$\vec{\Omega_{M}}$, $\vec{w}$ |  A distribution of 150 best-fitting cosmologies, produced by processing the 150 realisations of $\Omega_{M}^{\prime}$, $w^{\prime}$ with Pippin. Described in Section 3.1.
$w^{*}(\Omega_{M})$ |  A one dimensional function that approximates $\vec{\Omega_{M}}$, $\vec{w}$. Found by fitting a Gaussian process through $\vec{\Omega_{M}}$, $\vec{w}$. Described in Section 3.1.1.
Table 1: A glossary of terms used in our methodology, which are defined
throughout the text

For a given input cosmology, the Neyman construction provides a prescription
for using simulations to calculate the percentile contour, or the boundary of
a confidence region, that input cosmology lies on. By calculating these
percentile contours for a grid of input cosmologies, the confidence region can
be estimated as the set of input cosmologies which lie on percentile contours
less than or equal to the desired confidence level. The extensive compute time
of supernova simulations makes it difficult to densely sample the parameter
space, so we instead create an approximate Neyman construction by
strategically choosing a small number of input cosmologies at representative
locations on the 68% contour.

In Section 3.1 we describe how to calculate the percentile contour for a
single cosmological input. In Section 3.2 we describe how we find cosmological
inputs which lie on the 68% percentile contour, and how we estimate the
confidence region.

###  3.1 Calculating the percentile contour

To calculate the percentile contour at any given cosmology
($\Omega_{M}^{\prime}$, $w^{\prime}$, which as noted above is chosen to lie at
representative locations of the 68% contour), we simulate 150 SNe Ia data-sets
with $\Omega_{M}^{\prime}$, $w^{\prime}$ as the true cosmological input, using
the procedure described in Section 2.3. Each data-set is produced with a
different random seed, allowing for statistical fluctuations between
realisations. We analyse each of these data-sets with Pippin to produce a
distribution of best fitting cosmologies ($\vec{\Omega_{M}}$, $\vec{w}$) in
the space of measured $\Omega_{M}$ and $w$. The Neyman construction predicts
that the percentile contour for $\Omega_{M}^{\prime}$, $w^{\prime}$ is the
percentage of these best fitting cosmologies encompassed by a coverage ellipse
of this distribution. The coverage ellipse is defined to be centered on
$\Omega_{M}^{\prime}$, $w^{\prime}$ and to intersect $\Omega_{M}^{best}$,
$w^{best}$, which were calculated in Section 2.3, Equation 3. This ellipse
represents the probability of a data-set with true cosmology
$\Omega_{M}^{\prime}$, $w^{\prime}$ having a best fitting cosmology equal to
$\Omega_{M}^{best}$, $w^{best}$. Figure 4 shows an example of this calculation
for $\Omega_{M}^{\prime}=0.188$, $w^{\prime}=-0.783$, which lies at one
extreme end of the 68% contour we are testing. Figure 5 presents an example
Hubble Diagram for both $\Omega_{M}^{best}$, $w^{best}$ and
$\Omega_{M}^{\prime}$, $w^{\prime}$.

Figure 4: An example of using simulations to calculate the percentile contour
for $\Omega_{M}^{\prime}$, $w^{\prime}$, where $\Omega_{M}^{best}$, $w^{best}$
represent the best fitting cosmology for our test data-set. We simulate 150
data-sets using $\Omega_{M}^{\prime}$, $w^{\prime}$ as the input, and process
each data-set with Pippin to find the best fitting cosmology. The coverage
ellipse is defined to intersect $\Omega_{M}^{best}$, $w^{best}$. The
percentile contour for $\Omega_{M}^{\prime}$, $w^{\prime}$ is the percentage
of best fitting cosmologies contained within this coverage ellipse. Figure 5:
Top Panel: Hubble Diagram for $\Omega_{M}^{best}=0.3$, $w^{best}=-1.0$ and
$\Omega_{M}^{\prime}=0.188$, $w^{\prime}=-0.783$. This includes both simulated
distance moduli and the analytic distance moduli based on the input cosmology.
Bottom Panel: Difference between the analytic distance moduli of
$\Omega_{M}^{best}$, $w^{best}$ and $\Omega_{M}^{\prime}$, $w^{\prime}$.

####  3.1.1 Fitting the coverage ellipse

The distribution of best fitting cosmologies about $\Omega_{M}^{\prime}$,
$w^{\prime}$ usually follows the "banana" distribution that is typical of
supernova cosmology, and is due to the inherent degeneracy between
$\Omega_{M}$ and $w$. This contour shape makes it difficult to determine an
accurate coverage ellipse around this distribution. To determine a coverage
ellipse, we first fit a Gaussian Process (GP) through $\vec{\Omega_{M}}$,
$\vec{w}$ to produce the one-dimensional function $w^{*}(\Omega_{M})$, which
approximates $w$ as a function of $\Omega_{M}$ in the plane of
$\vec{\Omega_{M}}$, $\vec{w}$. Next we subtract $w^{*}(\vec{\Omega_{M}})$ from
$\vec{w}$ to transform the distribution of best fitting cosmologies into a
more elliptical distribution. We fit a coverage ellipse to this transformed
distribution which is centered on $\Omega_{M}^{\prime}$, $w^{\prime}$ and
intersects $\Omega_{M}^{best}$, $w^{best}$.

The percentile contour for $\Omega_{M}^{\prime}$, $w^{\prime}$ is then the
percentage of the best fitting cosmologies covered by this coverage ellipse.
Figure 6 shows an example of this transformation and ellipse fitting technique
for $\Omega_{M}^{\prime}=0.188$, $w^{\prime}=-0.783$, the same cosmology as
shown in Figure 4.

The uncertainty in the computed percentage is estimated by performing 1000
bootstrap resamples of the distribution of best fitting cosmologies, and
typically results in an uncertainty in the percentile contour of $\pm 4\%$.

Figure 6: Top panel: A GP fit ($w^{*}(\vec{\Omega_{M}})$) to the best fitting
output cosmologies of the 150 realisations with input cosmology:
$\Omega_{M}^{\prime}=0.188$, $w^{\prime}=-0.783$. Bottom panel: The same
distribution of maximum likelihood output cosmologies, transformed by
subtracting $w^{*}(\vec{\Omega_{M}})$ from $\vec{w}$. This transformed
distribution is more elliptical than the original distribution, and is more
appropriate for fitting coverage ellipses. We show one such coverage ellipse
in the bottom panel, scaled to intersect with the experiment cosmology input.
In this example 46% of the simulations are covered by the ellipse, so this
$\Omega_{M}^{\prime}$, $w^{\prime}$ lies on the 46% percentile contour.

###  3.2 Estimating the confidence region

We now have a statistically rigorous method of calculating the percentile
contour over the cosmological parameter space. Using this method, we could
compute a Neyman construction by computing the percentile contour across a
grid that covers the cosmological parameter space, and from this determine a
confidence region. However, each evaluation of the percentile contour requires
150 simulated data-sets, which is computationally expensive. To more
efficiently determine a confidence region, we develop an approximate Neyman
construction method which requires far fewer simulations. Instead of
evaluating the percentile contour over the entire cosmological parameter
space, we use the bisection method to iteratively find the input cosmologies
that lie on the 68% percentile contour. These cosmologies define the edge of
the 68% confidence region.

We first calculate the percentile contour for cosmologies at several
representative locations on the 68% contour. We select two input cosmologies
which are at the furthest extent of the 68% contour, and two input cosmologies
which are at the closest region of the 68% contour to $\Omega_{M}^{best}$,
$w^{best}$. This enables us to probe the consistency of the 68% contour across
its entire span. The cosmological inputs used in defining the approximate
Neyman construction, and the percentile contour for each input are shown in
Figure 7, and detailed in Table 2. Inputs 1a, 1b, 1c, and 1d were chosen to
lie on the 68% contour of the original experiment cosmology posterior. A
second set of inputs (2a, 2b, 2c, and 2d) were chosen to compensate for how
far the previous set of coverage ellipses differed from 68% coverage, as
described below.

Figure 7: The input $\Omega_{M}$ and $w$ values for the experiment
($\Omega_{M}^{best}$, $w^{best}$) and the approximate Neyman construction
($\Omega_{M}^{\prime}$, $w^{\prime}$).

If the input cosmologies lie at a percentile contour of $<68\%$ confidence, we
select a new cosmology further from $\Omega_{M}^{best}$, $w^{best}$, and
conversely select a cosmology closer to $\Omega_{M}^{best}$, $w^{best}$ if the
initial percentile contour is $>68\%$ confidence. This allows us to find an
input cosmology which lies on a percentile contour within one standard
deviation of 68%, as measured by bootstrap resampling. Though the iterative
method typically converges to a percentile contour within one standard
deviation of $68\%$ within 2 or 3 iterations, converging to exactly 68% would
take significantly more iterations. As such, once we are within one standard
deviation of 68%, we linearly interpolate or extrapolate to find an input
cosmology which lies exactly on the 68% percentile contour.

Figure 8 shows this approximate Neyman Construction technique for
$\Omega_{M}^{\prime}=0.188$, $w^{\prime}=-0.783$, the same cosmology shown in
Figures 4 and 6. This cosmology was found to lie on the $46\%\pm 4\%$
percentile contour, and after iteration the cosmology
$\Omega_{M}^{\prime}=0.145$, $w^{\prime}=-0.725$ was found to lie on the
$65\%\pm 4\%$ percentile contour, which is within one standard deviation of
68%. We linearly extrapolated from these two input cosmologies to find that
$\Omega_{M}^{\prime}=0.138$, $w^{\prime}=-0.716$ lay on the 68% percentile
contour. Figure 9 shows the coverage ellipses fit to each cosmological input.

Figure 8: Example of finding the edge of the 68% confidence region. $\Omega_{M}^{\prime}$, $w^{\prime}$ 1a was defined with an input cosmology on the extreme end of the Pippin 68% contour, and was found to lie on the $46\%\pm 4\%$ percentile contour. $\Omega_{M}^{\prime}$, $w^{\prime}$ $2a$ was found iteratively, and lies on the $65\%\pm 4\%$ percentile contour. Linearly extrapolating from these two cosmological inputs gives us $68\%$ $\Omega_{M}^{\prime}$, $w^{\prime}$. Cosmology | Input | Input | Percentile
---|---|---|---
Input | $\Omega_{M}$ | $w$ | Contour
Experiment | 0.3 | -1.0 | -
1a | 0.188 | -0.783 | $46\%\pm 4\%$
1b | 0.38 | -1.25 | $74\%\pm 4\%$
1c | 0.307 | -0.977 | $47\%\pm 4\%$
1d | 0.292 | -1.02 | $62\%\pm 4\%$
2a | 0.145 | -0.725 | $65\%\pm 4\%$
2b | 0.37 | -1.204 | $66\%\pm 4\%$
2c | 0.3075 | -0.9765 | $67\%\pm 4\%$
2d | 0.2917 | -1.0205 | $70\%\pm 4\%$
Table 2: The input $\Omega_{M}$ and $w$ values for the experiment and
approximate Neyman construction input cosmologies, as well as the percentile
contour each cosmological input lies on. Figure 9: Coverage ellipses fit to
the maximum likelihood distribution of each Neyman input cosmology,
transformed to $\left\\{\Omega_{m},w-w^{*}(\Omega_{M})\right\\}$. The coverage
ellipse is defined to be centered on the Neyman input, scaled such that it
contains the experiment cosmology input. The title of each plot shows the
percentage of maximum likelihood output cosmologies covered by the ellipse,
this is our numerical estimate of likelihood. The uncertainty in this estimate
is calculated via 1000 bootstrap resamples and is 4%.

##  4 Results

In the previous section we describe how we computed the location of the 68%
confidence region at four strategically selected cosmological inputs in the
$\Omega_{M}-w$ plane using an approximate Neyman construction. The results are
presented in Table 3 and shown in Figure 10.

We see the largest difference between the experiment cosmology 68% contour and
the 68% confidence region in the first input cosmology, which lies in the top
left quadrant of the cosmological contour. We find a shift of $\sim$0.05 in
$\Omega_{M}$ and a shift of $\sim$0.07 in $w$. The second input cosmology,
which lies in the bottom right quadrant has a shift of $\sim$0.007 in
$\Omega_{M}$ and $\sim$0.03 in $w$. These input cosmologies lie at the
furthest extent of the 68% cosmological contour from $\Omega_{M}^{best}$ and
$w^{best}$. In contrast, the third and fourth input cosmologies, which lie
much closer to the experiment input cosmology, have a shift of $\lesssim$0.001
in both $\Omega_{M}$ and $w$. The increase in discrepancy as we probe
parameter space that is further from the experiment input cosmology is
expected, as the BBC bias correction is only performed at a single point in
cosmological parameter space, the best-fitting cosmology. The bias correction
depends on the input cosmology, thus the BBC bias correction produces a
contour which is accurate close to the best-fitting cosmology, but induces an
offset in the contour at parameter space further from the best-fitting
cosmology.By contrast, our Neyman construction method applies this bias
correction in a cosmologically dependent manner across the entire parameter
space, removing this offset.

This offset is important to consider, especially when combining supernova
contours with other cosmological probes such as the cosmic microwave
background. Fortunately, it is the region of posterior space close to the
input cosmology which overlaps with the contours of other cross-cutting
probes. As such this offset is unlikely to significantly affect multi-probe
cosmological analyses.

Where this offset could be significant is in the investigation of cosmological
tensions, where accurate uncertainties of the tails are vital to successfully
assess the significance of the tension. In these cases, it may be useful to
use a method like our approximate Neyman construction to produce more accurate
and statistically rigorous measure of the uncertainty in a cosmological fit.

Figure 10: Comparison between the 68% confidence region determined from our approximate Neyman construction, and the 68% contour of the experiment cosmology. The confidence region is consistent with the contour close to the input cosmology, but displays an offset at the extreme ends of the contour. This offset is likely due to the bias correction method used by BBC, which is most accurate close to the input cosmology. Cosmology | $\Omega_{M}^{\prime}$ | $w^{\prime}$ |  $\Omega_{M}$ Absolute |  $w$ Absolute
---|---|---|---|---
Input |  |  | Difference | Difference
1 | 0.138 | -0.716 | 0.05 | 0.07
2 | 0.373 | -1.216 | 0.007 | 0.03
3 | 0.308 | -0.976 | 0.001 | 0.001
4 | 0.2918 | -1.02 | 0.0002 | 0.0
Table 3: Comparison between the 68% confidence region determined from our
approximate Neyman construction, and the 68% contour of the experiment
cosmology. The absolute difference is the difference between the cosmologies
at the edge of the 68% contour produced by Pippin, and the cosmologies at the
edge of the 68% confidence region produced by our approximate Neyman
construction.

###  4.1 Effect of Bias Correction Input Cosmology

We repeat our analysis with bias correction simulations that use a cosmology
that is nearer to the data-set they are attempting to correct, rather than a
single bias correction shared amongst all realisations. Figure 11, and Table 4
show the results of this repeat analysis. Our results are very similar to the
case when we shared the bias correction simulation amongst all realisations,
with the first and second cosmological input deviating the most. Overall,
these results reinforce our suggestion that the offset present between the
cosmological contour and confidence region is caused by the BBC bias
correction, however the choice of bias correction does not significantly
affect our consistency test. If computational cost is a concern, using only
one bias correction simulation shared amongst all realisations will
significantly reduce the computational cost of our approximate Neyman
construction method, without significantly reducing the quality of the
consistency test.

Figure 11: As per Figure 10, but varying the bias correction simulation to match the input cosmology. Very similar results are found, indicating that the cosmology used for the bias correction is not significantly impacting the results. Cosmology | $\Omega_{M}^{\prime}$ | $w^{\prime}$ |  $\Omega_{M}$ Absolute |  $w$ Absolute
---|---|---|---|---
Input |  |  | Difference | Difference
1 | 0.125 | -0.698 | 0.063 | 0.085
2 | 0.371 | -1.21 | 0.009 | 0.004
3 | 0.308 | -0.977 | 0.001 | 0.0
4 | 0.292 | -1.02 | 0.0 | 0.0
Table 4: As for Table 3, but varying the bias correction simulation to match
the input cosmology.

##  5 Conclusions

In this paper we present a statistically rigorous method for checking the
consistency of contours produced in a cosmological analysis. To achieve this,
we implement an approximate Neyman construction which requires far less
computation than a true Neyman construction. This approximate Neyman
construction is then used to define the 68% confidence region for a single
cosmological realisation. We use this confidence region to test the
consistency of the 68% contour produced by the BBC framework, as integrated in
Pippin, although this method can be used to test the consistency of any
cosmological parameter estimation method. This represents the first time the
BBC framework has been tested with a statistically rigorous methodology.

Our analysis showed that, for a DES-3YR like dataset, Pippin is producing
reasonable, consistent parameter estimates. There was some discrepancy between
the CR and the cosmological contour when considering the farthest extent of
the 68% contour. This discrepancy was, at maximum, a shift of $\sim$0.05 in
$\Omega_{M}$, and $\sim$0.07 in $w$, and was likely due to the accuracy of
BBC’s bias correction being best when close to the input cosmology, and
degrading in regions of parameter space which are far from the input
cosmology. It is also important to recognise that this does not correspond to
an equivalent error in the reported maximum posterior cosmological parameters.
When considering cosmological inputs close to the experiment cosmology input,
the confidence region and cosmological contour had near perfect agreement. As
such any overall discrepancy is unlikely to significantly effect the results
of a cosmological analysis, especially when multiple cross-cutting probes are
combined. However, this shift is important when considering analyses concerned
with assessing cosmological tensions - where the precise shape and size of the
contour are vitally important to the analysis.

We see very similar results when each realisation had its own bias correction
simulation, rather than sharing one bias correction simulation amongst all
realisations, indicating that a sensible choice of bias correction is not
likely to significantly effect our consistency checks.

Overall, we believe our method for consistency checking cosmological contours
with an approximate Neyman construction represents an important improvement in
the statistical rigour applied to cosmological analyses, and should become a
standard step in all cosmological analyses. Our methodology can also be used
to rigorously test cosmological contours for other cosmological probes, which
have similarly complex pipelines. We believe this method will be particularly
useful for future analyses, such as the DES 5-year supernova analysis, and the
upcoming LSST survey. We plan to repeat this analysis using simulations that
match the DES 5-year supernova analysis to test the consistency of those
results.

##  6 Acknowledgments

This is a pre-copyedited, author-produced PDF of an article accepted for
publication in the Publications of the Astronomical Society of Australia
following peer review. The version of record is available online at: .

We acknowledge parts of this research were carried out on the traditional
lands of the Ngunnawal, Turrbal, and Jagera people. We pay our respects to
their elders past, present, and emerging. P.A. was supported by an Australian
Government Research Training Program (RTP) Scholarship. H.Q. and M.S. were
supported by DOE grant DE-FOA-0002424 and NSF grant AST-2108094. B.E.T. and
his group were supported by the Australian Research Council Centre of
Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through
project number CE170100013. T.M.D. was supported by an Australian Research
Council Laureate Fellowship FL180100168.

###  6.1 Data availability

The data underlying this article are available in the article and in its
online supplementary material. The code used throughout this article is
available at https://github.com/dessn/BiasValidation.

###  6.2 Software

Pippin HintonBrout2020, SNANA KesslerBernstein2009, ChainConsumer
chainconsumer, Numpy numpy, MatPlotLib matplotlib, SciKit-Learn scikit-learn,
Scipy scipy

###  6.3 Facilities

This work was completed in part with resources provided by the University of
Chicago Research Computing Center.