###### Abstract

The oriented chromatic polynomial of a oriented graph outputs the number of
oriented $k$-colourings for any input $k$. We fully classify those oriented
graphs for which the oriented graph has the same chromatic polynomial as the
underlying simple graph, closing an open problem posed by Sopena. We find that
such oriented graphs can be both identified and constructed in polynomial time
as they are exactly the family of quasi-transitive oriented co-interval
graphs. We study the analytic properties of this polynomial and show that
there exist oriented graphs which have chromatic polynomials have roots,
including negative real roots, that cannot be realized as the root of any
chromatic polynomial of a simple graph.

Chromatic Polynomials of Oriented Graphs

Danielle Cox1 Christopher Duffy2

1Department of Mathematics, Mount Saint Vincent University, Halifax, CANADA
2Department of Mathematics, University of Saskatchewan, Saskatoon, CANADA

An _oriented graph_ arises by assigning directions to the edges of a simple
graph. For an oriented graph, $G$, we let $U(G)$ denote the underlying simple
graph. We say that $G$ is an _orientation_ of $U(G)$. Alternately, an oriented
graph is an irreflexive and anti-symmetric digraph.

The generalization of proper colouring to graph homomorphism provides a path
to define proper colouring for oriented graphs in a way that takes into
account the orientation. For an oriented graph $G$, an _oriented
$k$-colouring_ of $G$ is a homomorphism to a tournament (i.e., an orientation
of a complete graph) on $k$ vertices. One can see that this definition of
oriented colouring is equivalent the following one, which dispenses with the
need to invoke homomorphism. For an oriented graph $G=(V_{G},A_{G})$, a
function $c:V_{G}\to\\{1,2,\dots k\\}$ is an _oriented $k$-colouring_ when

  1. 1.

$c(u)\neq c(v)$ for all $uv\in A_{G}$, and

  2. 2.

for $uv,xy\in A_{G}$, if $c(u)=c(y)$, then $c(v)\neq c(x)$.

This second condition implies directly that non-adjacent vertices at the end
of a directed path of length two (a $2$-dipath) are assigned distinct colours
in any oriented colouring.

Since their introduction by Courcelle in his treatment of monadic second order
logic and graph structure [8], oriented colourings have provided a fertile
area for fundamental research in mathematics and theoretical computer science.
Many of the questions that have interested both applied and theoretical
researchers in the study of graph colourings find an analogue in the study of
oriented graphs. In addition to bounds for a variety of graph families [7, 9,
10, 11, 16], researchers have examined the computational complexity of related
decision problems [2, 13], the notion of clique for oriented graphs [3],
oriented arc-colourings [15], oriented list-colourings [23] and even an
oriented colouring game [14]. An excellent overview of the state-of-the-art is
given in [21].

For an oriented graph $G$, we define the _oriented chromatic polynomial_ to be
the unique interpolating polynomial $f_{o}(G,\lambda)$ so that $f_{o}(G,k)$ is
the number of $k$-colourings of $G$.

Figure 1 gives an oriented graph together with its oriented chromatic
polynomial.

Figure 1: An oriented graph and its oriented chromatic polynomial.

The oriented chromatic polynomial was first introduced by Sopena in [20]. Here
he established some basic properties of $f_{o}(G,k)$, provided a recursive
construction, and exhibited some oriented graphs for which the analytic
behaviour of the oriented chromatic polynomial diverged wildly from possible
behaviour of the chromatic polynomial of any graph.

The definition of the oriented chromatic polynomial follows the from
definition of the _chromatic polynomial_ , $f(\Gamma,\lambda)$, for a graph
$\Gamma$. The polynomial $f(\Gamma,\lambda)$ is defined to be the unique
interpolating polynomial so that $f(\Gamma,k)$ is the number of $k$-colourings
of $\Gamma$.

We call an oriented graph _chromatically invariant_ when
$f_{o}(G,\lambda)=f(U(G),\lambda)$. One can see that the oriented graph given
in Figure 1 is not chromatically invariant; the coefficients of
$f_{o}(G,\lambda)$ do not alternate sign, a feature of every chromatic
polynomial. This further implies, in fact, that there is no graph $\Gamma$
such that $f_{o}(G,\lambda)=f(\Gamma,\lambda)$. We say that a graph $\Phi$ and
an oriented graph $H$ are _chromatically equivalent_ if
$f_{o}(H,\lambda)=f(\Phi,\lambda)$.

The zeros of graph polynomials have been an active area of research for many
years. One such graph polynomial whose zeroes have been studied is the
chromatic polynomial. The real roots of such polynomials are dense in the
interval [32/27,$\infty)$ and there are no real roots in the interval
$(0,1)\cup(1,32/27]$[4, 5]. The complex roots of such polynomials are dense in
the complex plane [19].

Our work proceeds as follows. In the following section we review and
recontextualize the recursive construction of $f_{o}(G,k)$ given by Sopena for
the purposes of providing an explicit formula for the coefficient of
$\lambda^{n-2}$ in $f_{o}(G,k)$. In Section 2 we provide a full classification
of chromatically invariant oriented graphs, show they can be recognized in
polynomial time and further explore the relationship between the oriented
chromatic polynomial and the chromatic polynomial. We provide a partial answer
to the more general question of finding chromatically equivalent pairs of
oriented graphs and graphs. In Section 3 we study the roots of oriented
chromatic polynomials and show that there exist oriented graphs whose oriented
chromatic polynomials have roots that cannot be realized as a root of a
chromatic polynomial. In particular we exhibit oriented graphs whose
polynomials have negative real roots – a feature of no chromatic polynomial.

All graphs considered herein are simple. That is, we do not allow loops or
multiple edges. Further, in graphs with more than one type of adjacency, we
allow at most one type of adjacency between a pair of vertices. We refer the
reader to [6] for graph theoretic definitions and notation.

##  1 The Oriented Chromatic Polynomial

Let $G=(V,E,A)$ be a mixed graph. That is, $G$ is a graph in which a subset
(possibly empty) of the edges have been oriented to be arcs. We say that $c$
is an _oriented colouring_ of $G$ when $c$ is an oriented colouring when $G$
is restricted to the arcs and a proper colouring when $G$ is restricted to the
edges. Notice that if every pair of vertices is either adjacent or at the ends
of a directed path of length $2$, then every vertex must receive a distinct
colour in every colouring. We define the _oriented chromatic polynomial_ of a
mixed graph analogously to that of oriented graphs. In this section we observe
that the oriented chromatic polynomial introduced by Sopena can be generalized
as the oriented chromatic polynomial of mixed graphs. This generalization
allows us to find a closed form for the third coefficient of the oriented
chromatic polynomial of mixed graphs. As every oriented graph is a mixed graph
with an empty arc set, this expression leads us to a closed form for the third
coefficient of the chromatic polynomial of an oriented graph. The expression
for the third coefficient provides us with a tool to study chromatically
invariant oriented graphs.

We begin by providing a recursive formula for the oriented chromatic
polynomial of mixed graphs. Let $G=(V_{G},A_{G},E_{G})$ be a mixed graph. If
every pair of vertices is either adjacent or at the ends of a directed path of
length $2$ (a $2$-dipath), then

$f_{o}(G,\lambda)=\prod_{i=0}^{i=n-1}(\lambda-i)$

Otherwise, let $u$ and $v$ be a pair of vertices that are neither adjacent nor
at the ends of a directed path of length two. In this case we have

$f_{o}(G,\lambda)=f_{o}(G+uv,\lambda)+f_{o}(G_{uv},\lambda),$

where

  * •

$G+uv$ is the mixed graph formed from $G$ by adding an edge between $u$ and
$v$; and

  * •

$G_{uv}$ is the mixed graph formed from $G$ by identifying $u$ and $v$ into a
single vertex, deleting all parallel arcs and edges, and deleting any edge
that is parallel with an arc.

Following the usual convention of having the picture of a graph stand in for
its polynomial, an example of this recursion is provided in Figure 2.

Figure 2: Computing the oriented chromatic polynomial by way of reduction.

The correctness of this reduction follows from the proof of the reduction
provided in [20]. The appendix gives Maple code for generating the oriented
chromatic polynomial of a mixed graph.

For a mixed graph $G$, let $\mathcal{D}_{G}$ be the set of pairs of vertices
that are at the ends of an induced $2$-dipath. In Figure 1 we have
$\mathcal{D}_{G}=\left\\{\\{u_{4},u_{6}\\}\right\\}$.

######  Theorem 1.1.

[20] For any mixed graph $G=(V_{G},E_{G},A_{G})$ with $n$ vertices

  1. 1.

$f(G,\lambda)$ is a polynomial of order $n$ in $\lambda$;

  2. 2.

the coefficient of $\lambda^{n}$ is $1$;

  3. 3.

$f(G,\lambda)$ has no constant term;

  4. 4.

the coefficient of $\lambda^{n-1}$ is $-(|A_{G}|+|E_{G}|+|\mathcal{D}_{G}|)$;
and

  5. 5.

if $G$ has an isolated vertex $x$, then $f_{o}(G)=\lambda\cdot
f_{o}(G-x,\lambda)$.

Let $G=(V_{G},A_{G},E_{G})$ be an mixed graph with non-incident arcs $uv,xy\in
A_{G}$. We say that the arcs $uv,xy$ are _obstructing_ , when

  1. 1.

$u$ and $y$ are not the ends of a $2$-dipath;

  2. 2.

$v$ and $x$ are not the ends of a $2$-dipath and

  3. 3.

$uy,vx\notin E_{G}$.

Let $\mathcal{O}_{G}$ denote the set of unordered pairs of obstructing arcs in
$G$, a mixed graph. In Figure 1 we have

$\mathcal{O}_{G}=\left\\{\\{u_{1}u_{2},u_{4}u_{3}\\},\\{u_{1}u_{2},u_{4}u_{5}\\},\\{u_{1}u_{2},u_{5}u_{6}\\},\\{u_{1}u_{2},u_{7}u_{6}\\}\right\\}.$

######  Lemma 1.2.

If $\mathcal{O}_{G}\neq\emptyset$ or $\mathcal{D}_{G}\neq\emptyset$, then
$f_{o}(G,\lambda)\neq f(U(G),\lambda)$.

###### Proof.

Observe that every colouring of $G$ using $k$ colours is a colouring of
$U(G)$. However if $G$ has an induced $2$-dipath or a pair of obstructing
arcs, then the converse does not hold. As such, there exists $k$ such that
$f_{o}(G,k)<f(U(G),k)$. ∎

Let $c_{i}(f,\lambda)$ be the coefficient of $\lambda^{n-i}$ in
$f(G,\lambda)$.

######  Theorem 1.3.

For a mixed graph $G=(V_{G},A_{G},E_{G})$, we have

$c_{2}(G,\lambda)={|A_{G}|+|\mathcal{D}_{D}|+|E_{G}|\choose
2}-|T_{G}|-|\mathcal{D}_{G}|-|\mathcal{O}_{G}|,$

where $T_{G}$ is the set of induced subgraphs isomorphic to $K_{3}$ in $U(G)$.

###### Proof.

Let $G$ be a minimum counter-example with respect to number of vertices. Among
all such counter examples, let $G$ be the one that maximizes
$|A_{G}|+|E_{G}|+|\mathcal{D}_{G}|$. Note that we may further assume
$\mathcal{D}_{G}=\emptyset$ by adding an edge between every pair of vertices
in $\mathcal{D}_{G}$. The resulting mixed graph has the same set of oriented
colourings, and thus the same oriented chromatic polynomial.

The oriented chromatic polynomial of a mixed complete graph on $n$ is equal to
the chromatic polynomial of a complete graph on $n$ vertices. The third
coefficient of such a graph is given by ${{n\choose 3}\choose 2}-{n\choose 3}$
[18]. Therefore the claim holds for mixed complete graphs. As $G$ is a minimum
counter example, $G$ is not a complete mixed graph. As such there exists
$u,v\in V(G)$ such that $u$ and $v$ are not adjacent, nor at the ends of a
$2$-dipath. Therefore
$c_{2}(G,\lambda)=c_{2}(G+uv,\lambda)+c_{1}(G_{uv},\lambda)$. By the choice of
$G$, the claim holds for both $G+uv$ and $G_{uv}$.

Let $C$ be the set of common neighbours of $u$ and $v$ in $U(G)$. Each of
these common neighbours forms a triangle in $G+uv$. And so
$|T_{G}|=|T_{G+uv}|+|C|$. Further observe that in $G_{uv}$, the arcs/edges
from $u$ and $v$ to a common neighbour $c\in C$ becomes a single adjacency in
$G_{uv}$. Therefore $|A_{G}|+|E_{G}|=|A_{G_{uv}}|+|E_{G_{uv}}|-|C|$.

A pair of obstructing arcs in $G$ is not obstructing in $G+uv$ if and only if
$u$ and $v$ are the head and tail, in some order, of the pair of obstructing
arcs. Let $\mathcal{O}_{G}^{uv}$ be the set of such arcs. No new obstructing
arcs can be created by adding an edge. Therefore
$|\mathcal{O}_{G}|-|\mathcal{O}_{G}^{uv}|$ = $|\mathcal{O}_{G+uv}|$. Also note
that a pair of obstructing arcs in $|\mathcal{O}_{G}^{uv}|$ form a $2$-dipath
in $G_{uv}$ with centre vertex $uv$. All other induced $2$-dipaths in $G$ are
retained in $G_{uv}$, as $u$ and $v$ are not the ends of a $2$-dipath.
Therefore $|D_{G}|+|\mathcal{O}_{G}^{uv}|=|D_{G_{uv}}|$.

$c_{2}(G,\lambda)=c_{2}(G+uv,\lambda)+c_{1}(G_{uv},\lambda)\\\ $ (1)
$={|A_{G}|+|D_{G}|+|E_{G}|+1\choose
2}-(|T_{G+uv}|+|D_{G+uv}|+|\mathcal{O}_{G+uv}|)$
$-(|A_{G_{uv}}|+|E_{G_{uv}}|+|D_{G_{uv}}|)$
$={|A_{G}|+|D_{G}|+|E_{G}|+1\choose
2}-(|T_{G}|-|C|+|D_{G}|+|\mathcal{O}_{G}|-|\mathcal{O}_{G}^{uv}|)$
$-(|A_{G}|+|E_{G}|+|C|+|D_{G}|+|\mathcal{O}_{G}^{uv}|)$

$={|A_{G}|+|D_{G}|+|E_{G}|+1\choose
2}-(|A_{G}|+|D_{G}|+|E_{G}|)-|T_{G}|-|D_{G}|-|\mathcal{O}_{G}|$

$={|A_{G}|+|D_{G}|+|E_{G}|\choose 2}-|T_{G}|-|D_{G}|-|\mathcal{O}_{G}.|$

Thus the claim holds for $G$, contradicting the choice of $G$ as a minimum
counter example. ∎

For the case $E=\emptyset$, we arrive at the desired result for oriented
graphs.

######  Corollary 1.4.

For an oriented graph $G$, we have

$c_{2}(G,\lambda)={|A_{G}|+|\mathcal{D}_{G}|\choose
2}-|T_{G}|-|\mathcal{D}_{G}|-|\mathcal{O}_{G}|.$

We further note that in the case $A=\emptyset$, we arrive at the usual result
for the third coefficient of the chromatic polynomial:
$c_{2}(G,\lambda)={|E_{G}|\choose 2}-|T_{G}|$.

##  2 Oriented Chromatic Equivalence

A folklore construction gives an orientation $G$ of $K_{n,n}$ so that the
resulting oriented graph has chromatic number $2n$ (see [9]). This common
example is used to convince the reader that the oriented chromatic number and
the chromatic number of the underlying simple graph can be arbitrarily far
apart. We note, however that the set of colourings of $G$ using $\lambda\geq
2n$ colours is exactly that of colourings of $K_{2n}$ using $\lambda$ colours.
And so though the chromatic number of $U(G)$ differs greatly to $G$, there is
still a relationship between colourings of $G$ and colourings of some simple
graph. In this section we find a set of sufficient conditions so that the
$\lambda$-colourings of an oriented graph $G$ are exactly those of some simple
graph $\Gamma$. We conclude this section by using these sufficient conditions
to compute the chromatic polynomial of orientations of stars.

We are interested in the following decision problems:

_CHROM-INVAR_
 _Instance:_ A graph $\Gamma$.
 _Question:_ Is there an orientation $O(\Gamma)$ such that
$f_{o}(O(\Gamma),\lambda)=f(\Gamma,\lambda)$?

_OCHROM-INVAR_
 _Instance:_ An oriented graph $G$.
 _Question:_ Does $f_{o}(G,\lambda)=f(U(G),\lambda)$?

_OCHROM-EQUIV_
 _Instance:_ An oriented graph $G$.
 _Question:_ Is there a graph $\Gamma$ such that
$f_{o}(G,\lambda)=f(\Gamma,\lambda)$?

Figure 3 gives an example of an oriented graph $G$ and a graph $\Gamma$ so
that $G$ and $\Gamma$ are chromatically equivalent.

Figure 3: An oriented graph and a graph with the same chromatic polynomial.

Let $G$ be an oriented graph. Let $G^{\star}$ be the mixed graph resulting
from $G$ by adding an edge between $u$ and $v$ whenever $u$ and $v$ are at the
ends of an induced $2$-dipath. We observe the following.

######  Lemma 2.1.

$f_{o}(G,\lambda)=f_{o}(G^{\star},\lambda)$.

###### Proof.

$G$ is a subgraph of $G^{\star}$. Therefore $f_{o}(G,\lambda)\leq
f_{o}(G^{\star},\lambda)$. Every oriented colouring of $G^{\star}$ using $k$
colours is also an oriented colouring of $G$, therefore $f_{o}(G,\lambda)\geq
f_{o}(G^{\star},\lambda)$.
∎

######  Theorem 2.2.

For $G$, an oriented graph, $f_{o}(G,\lambda)=f(U(G^{\star}),\lambda)$ if and
only if $\mathcal{O}_{G}=\emptyset$.

###### Proof.

By Lemma 2.1 it suffices to show
$f_{o}(G^{\star},\lambda)=f(U(G^{\star}),\lambda)$. Notice that $G$ has no
obstructing arcs if and only if $G^{\star}$ has no obstructing arcs. By Lemma
1.2 it suffices to show that if $G^{\star}$ has no obstructing arcs, then
every colouring of $U(G^{\star})$ is an oriented colouring of $G$. Let $c$ be
a colouring of $U(G^{\star})$. Since $c$ is a proper colouring, if $c$ is not
an oriented colouring of $G^{\star}$, then the second condition of oriented
colouring as been violated. However, this not possible as $G^{\star}$ has
neither an induced $2$-dipath nor a pair of obstructing arcs. ∎

Recall the result of Corollary 1.4. If $G$ has no pair of obstructing arcs,
then in $G^{\star}$ we have $\mathcal{O}_{G}=\mathcal{D}_{G}=\emptyset$. And
so $c_{2}(G^{\star},\lambda)={|A_{G}|+|E_{G}|\choose 2}-|T_{G}|$. Notice that
this is exactly the third coefficient of the chromatic polynomial of
$U(G^{\star})$.

######  Corollary 2.3.

An oriented graph $G$ is chromatically invariant if and only if $G$ has no
induced $2$-dipath and $U(G)$ is $2K_{2}$-free.

###### Proof.

Assume $G$ has no induced $2$-dipath and that $U(G)$ is $2K_{2}$-free. It
follows directly that $G$ has no obstructing arcs and that $G=G^{\star}$. The
conclusion follows by Theorem 2.2.

Let $G$ be a chromatically invariant oriented graph. By definition, every
proper $k$-colouring of $U(G)$ is an oriented colouring of $G$. Therefore $G$
has no induced $2$-dipath, nor does $G$ contain a pair of obstructing arcs.
Since $G$ contains no pair of obstructing arcs, if $U(G)$ contains an induced
copy of $2K_{2}$, say $uv,xy$, then without loss of generality, there must be
an induced $2$-dipath between $u$ and $y$. This contradicts that $G$ contains
no induced $2$-dipath. Therefore $G$ has no induced $2$-dipath and $U(G)$ is
$2K_{2}$-free. ∎

Introduced by Ghouila-Houri, oriented graphs that contain no induced
$2$-dipath are called _quasi-transitive_.

######  Theorem 2.4.

[12] A graph $\Gamma$ admits a quasi-transitive orientation if and only if
$\Gamma$ is a comparability graph.

Notice that the family of $2K_{2}$-free comparability graphs is exactly the
family of co-interval graphs. And so combining Corollary 2.3 and Theorem 2.4
yields the following classification.

######  Theorem 2.5.

An oriented graph is chromatically invariant if and only if it is a quasi-
transitive orientation of a co-interval graph.

This theorem fully classifies chromatically invariant oriented graphs, as well
as those graphs that admit a chromatically invariant orientation. This closes
an open problem posed by Sopena in [20] and provides a geometric
interpretation of chromatically invariant oriented graphs as co-interval
graphs.

Theorem 2.5 implies that the decision problems _CHROM-INVAR_ and _OCHROM-
INVAR_ can be restated in terms of co-interval graph recognition. As co-
interval graphs can be identified in linear time [22], we arrive at the
following classification of _CHROM-INVAR_ and _OCHROM-INVAR_.

######  Corollary 2.6.

The decision problems _CHROM-INVAR_ and _OCHROM-INVAR_ are _Polynomial_.

We note that such an orientation of a co-interval graph need not be unique (up
to converse). There are many methods in the literature (for example see [1])
that give a quasi-transitive orientation of a comparability graph. A common
element of these methods is the construction of an auxiliary graph, $Aux(G)$,
so that a $2$-colouring of $G$ corresponds to a quasi-transitive orientation.
Such constructions imply that a comparability graph has a unique quasi-
transitive orientation (up to converse) if and only if $Aux(G)$ is connected.
As such constructions can be carried out in polynomial time, we find that
given a co-interval graph $\Gamma$, one may find in polynomial time an
orientation of $\Gamma$, $O(\Gamma)$, so that
$f(\Gamma,\lambda)=f_{o}(O(\Gamma),\lambda)$.

We now consider an application of Theorem 2.2 and find the oriented chromatic
polynomial of the family of orientations of stars. Let $S_{i,o}$ be the
orientation of a star on $i+o+1$ vertices, with centre vertex $x$, so that $x$
has $i$ in-neighbours and $o$ out-neighbours.

######  Corollary 2.7.

$f_{o}(S_{i,o})=\lambda\cdot f(K_{i,o},\lambda-1)$

###### Proof.

Observe that $S_{i,o}^{\star}$ has no obstructing arcs. Further observe that
$U(S_{i,o}^{\star})$ consists of a copy of $K_{i,o}$ together with a universal
vertex. By Theorem 2.2, we have $f_{o}(S_{i,o})=\lambda\cdot
f(K_{i,o},\lambda-1)$.
∎

Conversely, using the results from Section 1, one can find families of
oriented graphs for which there is no chromatically equivalent graph. Let $G$
be an orientation of $tK_{2}$ for some $t>1$. Recalling the notation of the
previous section we have

  * •

$|A_{G}|=t$;

  * •

$|E_{G}|=|\mathcal{D}_{G}|=0$;

  * •

$|\mathcal{O}_{G}|={t\choose 2}$; and

  * •

$|T_{G}|=0$.

If there exists $\Gamma$ such that $f_{o}(G,\lambda)=f(\Gamma,\lambda)$, then
by Theorem 1.1 and Corollary 1.4, it must be that $\Gamma$ has $2t$ vertices,
$t$ edges and ${t\choose 2}$ copies of $K_{3}$. A simple counting argument
implies that no such $\Gamma$ can exist.

######  Conjecture 2.8.

Let $G$ be an oriented graph. There is a graph $\Gamma$ such that $G$ and
$\Gamma$ are chromatically equivalent if and only if $G$ has no obstructing
pairs of arcs.

When restricted to oriented graphs that contain no pair of obstructing arcs,
_OCHROM-EQUIV_ is Polynomial – every instance is a YES instance. However, for
arbitrary inputs is not clear if _OCHROM-EQUIV_ is contained in NP, as
constructing the chromatic polynomial of a graph is NP-hard. We conjecture,
however, that those oriented graphs with no obstructing arcs are the only
oriented graphs whose $\lambda$-colourings have a one-to-one correspondence
with the $\lambda$-colourings of some graph $\Gamma$. If Conjecture 2.8 is
true, then _OCHROM-EQUIV_ is Polynomial for arbitrary inputs.

Conjecture 2.8 is not true when we allow $G$ to be a mixed graph, even when we
require $A\neq\emptyset$. Let $H$ be the mixed graph formed from a pair of
disjoint arcs by adding an edge between the heads and an edge between the
tails. Using the reduction given in Section 1, we find
$f_{o}(H,\lambda)=\lambda^{4}-4\lambda^{3}+5\lambda^{2}-2\lambda$. Notably,
this is exactly the chromatic polynomial of the example, $G$, given in Figure
2. The oriented graph $G$ has $\mathcal{D}_{G}=\emptyset$. Therefore
$G=G^{\star}$. Further $\mathcal{O}_{G}=\emptyset$. And so by Theorem 2.2 we
have $f(U(G),\lambda)=f_{o}(G,\lambda)=f_{o}(H,\lambda)$. In particular, $H$
has a pair of obstructing arcs, but yet there is graph $\Gamma$ such that
$f_{o}(H,\lambda)=f(\Gamma,\lambda)$. From this example, one may generate an
example on $n$ vertices for any $n\geq 4$ by repeatedly adding universal
vertices to $H$ and $U(G)$.

##  3 Roots of Oriented Chromatic Polynomials

The location of the roots of polynomials has been well studied for a variety
of graph polynomials, such the independence, domination, reliability and
chromatic polynomials. In this section we provide results regarding the roots
of the oriented chromatic polynomial. Chromatic polynomials have roots that
are dense in the complex plane [19]. Their coefficients alternate in sign and
hence have no negative real roots [17]. We show the following:

######  Theorem 3.1.

For every integer $k>0$, there exists an oriented graph $G$ so that
$f_{o}(G,\lambda)$ has a root $k^{\prime}$ so that $k^{\prime}<-k$.

###### Proof.

Let $D_{n}$ be the oriented graph on $n$ vertices that consists of a directed
path $v_{1},v_{2},v_{3},v_{4}$ with $n-4$ leaves, $x$ directed from $x$ to
$v_{4}$.

The oriented chromatic polynomial can be computed in the following way. Vertex
$v_{1}$ has $\lambda$ colour choices, $v_{2}$ has $\lambda-1$ colour choices
and $v_{3}$ has $\lambda-2$ colour choices. Now for $v_{4}$ we have two
options. If $v_{4}$ is the same colour as $v_{1}$ then the $n-4$ leaves have
$\lambda-2$ choices of colour, since it can not be the same colour at $v_{4}$
or $v_{2}$. If $v_{4}$ is a different colour than the others in the path it
has $\lambda-3$ colour choices and the leaves have $\lambda-1$ colour choices.
This means

$f_{o}(D_{n},\lambda)=\lambda(\lambda-1)(\lambda-2)((\lambda-2)^{n-4}+(\lambda-3)(\lambda-1)^{n-4})$

The polynomial $f_{o}(D_{n},\lambda)$ has real roots at $\lambda=0,1,2$. We
show we can obtain arbitrarily large negative real roots by showing that a
real root exists between $\lambda=-n$ and $\lambda=-\ln(n)$ for $n$ even.

Observe,

$f_{o}(D_{n},-n)=(-1)^{3}n(n+1)(n+2)((-1)^{n-4}(n+2)^{n-4}+(-1)^{n-3}(n+3)(n+1)^{n-4})$

It can be shown that $(n+2)^{n-4}<(n+3)(n+1)^{n-4}$ if

$\Big{(}1+\frac{1}{n+1}\Big{)}^{n-4}<n+3.$

The quantity $\Big{(}1+\frac{1}{n+1}\Big{)}^{n-4}$ is bounded above by $e$,
and $e<n+3$ for all values of $n$, therefore $f_{o}(D_{n},-n)>0$ for all even
values of $n$.

Now consider $f_{o}(D_{n},-\ln(n))$.

$\displaystyle f_{o}(D_{n},-\ln(n))$ $\displaystyle=$
$\displaystyle(-1)^{3}(\ln(n))(\ln(n)+1)(\ln(n)+2)*$
$\displaystyle\Big{(}(-1)^{n-4}(\ln(n)+2)^{n-4}+(-1)^{n-3}(\ln(n)+3)(\ln(n)+1)^{n-4}\Big{)}$

Clearly $(\ln(n)+2)^{n-4}>(\ln(n)+3)(\ln(n)+1)^{n-4}$ when

$\Big{(}1+\frac{1}{\ln(n)+1}\Big{)}^{n-4}>\ln(n)+3.$

Let $g(n)=\Big{(}1+\frac{1}{\ln(n)+1}\Big{)}^{n-4}-\ln(n)-3$. The derivative
of this function is
$g^{\prime}(n)=\left(1+\left(\ln\left(n\right)+1\right)^{-1}\right)^{n-4}\left(\ln\left(1+\left(\ln\left(n\right)+1\right)^{-1}\right)-{\frac{n-4}{\left(\ln\left(n\right)+1\right)^{2}n\left(1+\left(\ln\left(n\right)+1\right)^{-1}\right)}}\right)-\frac{1}{n}$.

It is the case that $\lim_{n\to\infty}g^{\prime}(n)=\infty$, thus there exists
$N$ so that for all $n>N$, $g(n)$ is an increasing function, as the derivative
of $g$ is positive and hence $g(n)>0$ and $f_{o}(D_{n},-\ln(n))<0$ for large
values of $n$.

It then follows by the intermediate value theorem that $f_{o}(D_{n},\lambda)$
can have an arbitrarily large negative root. ∎

A chromatic polynomial cannot have root in the interval
$(-\infty,0)\cup(0,1)\cup(1,\frac{32}{27})$ [17]. We have shown that oriented
chromatic polynomials can have negative real roots. In addition, there exist
oriented chromatic polynomials in the interval $(0,1)$, as
$f_{o}(D_{5},\lambda)$ has a root at $\lambda=\frac{3}{2}-\frac{\sqrt{5}}{2}$.
Open problems regarding the roots of oriented chromatic polynomials include:
does there exist an oriented graph whose real roots lie in
$(1,\frac{32}{27})$? What is the closure of the complex roots for the oriented
chromatic polynomial?

##  4 Conclusion

The study of oriented graphs often goes hand-in-hand with that of signed
graphs. Though the methods contained herein will extend to the study of
chromatic polynomials of signed graphs, there will be a marked difference in
the classification of chromatically invariant signed graphs. For example,
letting all the edges of $2K_{2}$ be positive leads to a chromatically
invariant signed graph. However, all possible orientations of $2K_{2}$ leads
to an oriented graph that is not chromatically invariant. Of course, with this
approach, every graph can have edge signs trivially assigned so that the
resulting signed graph has the same chromatic polynomial as the underlying
graph. And so one may require that there is at least one edge of each sign.
With this added restriction it is unclear if chromatically invariant signed
graphs can be identified in polynomial time, as we expect the characterization
to require that signs be given so that there is no $2K_{2}$ where the edges
have different signs. Similarly, the generalization of signed graphs and
oriented graphs to $(m,n)$-mixed graphs should yield a definition of a
chromatic polynomial that obeys the reduction outlined in Section 1.
Consequently we expect the results of Theorems 1.1 and 1.3 to generalize in
the same manner. One may also ask, then, for which graphs $\Gamma$ is there an
assignment of arcs, edges, and corresponding colours, so that the resulting
$(m,n)$-mixed graph is chromatically invariant. We have shown that for
oriented graphs that $\Gamma$ must be a co-interval graph.

## Acknowledgments

The authors thank Gary MacGillivray for discussions regarding the recursive
construction of the oriented chromatic polynomial.

## References

  * [1] J. Bang-Jensen and G. Z. Gutin.  Digraphs: theory, algorithms and applications.  Springer Science & Business Media, 2008.
  * [2] J. Bang-Jensen, P. Hell, and G. MacGillivray.  The Complexity of Colouring by Semicomplete Digraphs.  SIAM Journal on Discrete Mathematics, 1(3):281–298, 1988.
  * [3] J. Bensmail, C. Duffy, and S. Sen.  Analogues of cliques for $(m,n)$-colored mixed graphs.  Graphs and Combinatorics, 33(4):735–750, 2017.
  * [4] B.Jackson.  A zero-free interval for the chromatic polynomials of graphs.  Combinatorics, Probability and Computing, 2:325–336, 1993.
  * [5] B.Jackson.  The zero-free interval for the chromatic polynomials of graphs.  Combinatorics, Probability and Computing, 6:497–506, 1997.
  * [6] J. Bondy and U. Murty.  Graph Theory.  Number 244 in Graduate Texts in Mathematics. Springer, 2008.
  * [7] O. Borodin, A. Kostochka, J. Nešetřil, A. Raspaud, and E. Sopena.  On the Maximum Average Degree and the Oriented Chromatic Number of a Graph.  Discrete Mathematics, 206:77–89, 1999.
  * [8] B. Courcelle.  The Monadic Second Order Logic of Graphs VI: On several representations of graphs by relational structures.  Discrete Applied Mathematics, 54:117–149, 1994.
  * [9] C. Duffy, G. MacGillivray, and E. Sopena.  Oriented colourings of graphs with maximum degree three and four.  Discrete Mathematics.  To appear.
  * [10] J. Dybizbański and A. Szepietowski.  The oriented chromatic number of Halin graphs.  Information Processing Letters, 114(1):45–49, 2014.
  * [11] G. Fertin, A. Raspaud, and A. Roychowdhury.  On the Oriented Chromatic Number of Grids.  Inform. Proc. Letters, 85:261–266, 2003.
  * [12] A. Ghouila-Houri.  Caractérisation des graphes non orientés dont on peut orienter les arětes de manière à obtenir le graphe d’une relation d’ordre.  C. R. Acad. Sci. Paris, 254:1370 – 1371, 1962.
  * [13] W. F. Klostermeyer and G. MacGillivray.  Pushing vertices and oriented colorings.  Bull. Inst. Combin. Appl, 40:49–58, 2004.
  * [14] J. Nešetril and E. Sopena.  On the oriented game chromatic number.  Electron. J. Combin, 8(2):R14, 2001.
  * [15] P. Ochem, A. Pinlou, and E. Sopena.  On the Oriented Chromatic Index of Oriented Graphs.  Journal of Graph Theory, 57(4):313–332, 2008.
  * [16] A. Raspaud and E. Sopena.  Good and Semi-Strong Colorings of Oriented Graphs.  Information Processing Letters, 51:171–174, 1994.
  * [17] R. Read and W. Tutte.  Selected Topics in Graph Theory 3, Chapter: Chromatic Polynomial.  Number 244. Academic Press Ltd., 1988.
  * [18] R. C. Read.  An introduction to chromatic polynomials.  Journal of Combinatorial Theory, 4(1):52–71, 1968.
  * [19] A. Sokal.  Chromatic roots are dense in the whole complex plane.  Combinatorics, Probability and Computing, 13(2):221–261, 2004.
  * [20] E. Sopena.  Computing chromatic polynomial of oriented graphs.  Université de Bordeaux 1/CNRS. Laboratoire Bordelais de Recherche en Informatique [LaBRI], 1994.
  * [21] E. Sopena.  Homomorphisms and Colourings of Oriented Graphs: An updated survey.  Discrete Mathematics, 339(7):1993–2005, 2016.
  * [22] J. P. Spinrad.  Efficient graph representations.  American Mathematical Society, 2003.
  * [23] Z. Tuza and M. Voigt.  Oriented list colorings of graphs.  Journal of Graph Theory, 36(4):217–229, 2001.

## Appendix

See pages - of MapleCode


# Accelerating Parallel Write via Deeply Integrating Predictive Lossy
Compression with HDF5

Sian Jin⋆, Dingwen Tao⋆, Houjun Tang⋆, Sheng Di2, Suren Byna3, Zarija Lukic3,
Franck Cappello2
⋆School of Electrical Engineering and Computer Science, Washington State
University, Pullman, WA, USA
2Mathematics and Computer Science Division, Argonne National Laboratory,
Lemont, IL, USA
3Computational Research Division, Lawrence Berkeley National Lab, Berkeley,
CA, USA Corresponding author: Dingwen Tao (dingwen.tao@wsu.edu).

###### Abstract

Lossy compression is one of the most efficient solutions to reduce storage
overhead and improve I/O performance for HPC applications. However, existing
parallel I/O libraries cannot fully utilize lossy compression to accelerate
parallel write due to the lack of deep understanding on compression-write
performance. To this end, we propose to deeply integrate predictive lossy
compression with HDF5 to significantly improve the parallel-write performance.
Specifically, we propose analytical models to predict the time of compression
and parallel write before the actual compression to enable compression-write
overlapping.  We also introduce an extra space in the process to handle
possible data overflows resulting from prediction uncertainty in compression
ratios.  Moreover, we propose an optimization to reorder the compression tasks
to increase the overlapping efficiency. Experiments with up to 4,096 cores
from Summit show that our solution improves the write performance by up to
4.5$\times$ and 2.9$\times$ over the non-compression and lossy compression
solutions, respectively, with only 1.5% storage overhead (compared to original
data) on two real-world HPC applications.

##  I Introduction

Large-scale scientific simulations on HPC systems play an important role in
today’s science and engineering domains. Such simulations can generate
extremely large amounts of data that are highly compute and storage intensive.
For example, one Nyx [1] cosmological simulation with a resolution of
$4096\times 4096\times 4096$ cells can generate up to 2.8 TB of data for a
single snapshot; a total of 2.8 PB of disk storage is needed, assuming the
simulation runs 5 times with 200 snapshots dumped per simulation.  Nowadays,
the ever-increasing computation power can be utilized to run the simulations.
However, managing such large amounts of data remains a major challenge.  It is
impractical to save all the generated raw data to disk due to: (1) the limited
storage capacity even for supercomputers, and (2) the I/O bandwidth required
to save these data can create bottlenecks in the transmission [2, 3, 4].

Lossy compression has been identified as one of the major data reduction
techniques to address this issue. Specifically, a new generation of error-
bounded lossy compression techniques, such as SZ [5, 6, 7], ZFP [8], and MGARD
[9], have been widely used in the scientific community [6, 5, 8, 7, 10, 11,
12, 4, 13, 14]. Compared to lossless compression that typically achieves only
$2\times$ of compression ratio [15] on scientific data, error-bounded lossy
compressors provide much higher compression ratio with controllable loss of
accuracy.

Scientific applications running on Supercomputers typically use parallel I/O
libraries such as Hierarchical Data Format 5 (HDF5) [16] for managing their
data. In specific, HDF5 is considered to provide high parallel I/O
performance, portability of data, and rich APIs for managing data on these
systems. It has been heavily used at supercomputing facilities for storing,
reading, and querying scientific datasets [17, 18]. This is because HDF5 has
specific designs and performance optimizations for popular parallel file
systems such as Lustre [19, 20]. In addition, instead of using a general
database in distributed storage, these datasets have their specific data
management approach based on the parallel file system [19, 20]. Moreover, HDF5
also provides users dynamically loaded filters [21] such as lossless and lossy
compression [22], which can automatically store and access data in compressed
formats. Thus, it allows scientific applications to store and access the data
in compressed formats. Parallel I/O in HDF5 with lossy compression filters can
not only significantly reduce data size, but also improve the overall I/O
performance.

Figure 1: Compression bit-rate distribution on a Nyx dataset with 512
partitions. Every partition uses the same compression configuration.

However, the existing implementation of HDF5 with compression filters cannot
fully utilize the benefits of reductions on data size and parallel-write time.
This is because to write the data from different processes to a shared file,
one must specify the offset of each data partition before writing the data.
Thus, the parallel write cannot start until all compression tasks finish and
the compressed sizes are communicated among all processes. Moreover, the
compressed size of each data partition is distributed across a wide range of
bit-rates (i.e., bits/value) and restricts any simple pre-allocation strategy.
Figure 1 shows the distribution of compressed bit-rates on a Nyx dataset. When
using lossless compression filters, the parallel-write performance could be
even lower than the original non-compression solution due to this high
overhead [23].

To solve this issue, we propose a parallel write solution that integrates
predictive lossy compression with the asynchronous I/O feature in HDF5, which
overlaps I/O latency with compression to significantly improve the parallel
write performance [24, 19]. Inspired by the previous research of Jin et al,
[25] that estimates the compression ratio of prediction-based lossy
compression with little overhead, we propose to predict the time of
compression and parallel write before the actual compression, and leverage the
asynchronous I/O feature to overlap compression with write. We also introduce
an extra space to handle the uncertainty of the prediction. Moreover, we
propose an optimization algorithm to reorder each process’s compression tasks
to increase the overlapping efficiency. The contributions of this paper are
summarized as follows.

  * •

We extend the prediction model for compression ratio to predict the
throughputs of compression and parallel write for prediction-based lossy
compression such as SZ.

  * •

We propose a new compression-write scheme with HDF5 that can efficiently write
the compressed data from different processes to a shared file by overlapping
compression with write based on our prediction models.

  * •

We optimize the execution order of compression tasks in each process to
achieve higher parallel-write performance.

  * •

We evaluate our proposed solution on two real-world HPC applications with up
to 4,096 cores on Summit supercomputer and up to 512 cores on Bebop cluster.
Experiments demonstrate that our solution improves the parallel-write
performance by up to $4.5\times$ and $2.9\times$ compared to the HDF5 write
without compression and with the SZ lossy compression filter (called “H5Z-SZ”)
[22], respectively, with only $1.5\%$ storage overhead.

In Section II, we discuss the research background. In Section III, we present
our design of parallel write with compression. In Section IV, we present our
evaluation results. In Section V, we conclude our work and discuss future
work.

##  II Research Background

In this section, we present the background information on HDF5 and lossy
compression and discuss the challenges.

###  II-A Parallel I/O Libraries for HPC Applications

Figure 2: Scientific data management with compression.

HPC applications generate and analyze massive amounts of data. A critical
requirement of these applications is the capability to access and manage this
data efficiently on HPC systems. Parallel I/O becomes the key technology to
enable efficiently moving data between compute nodes and storage considering
the complex storage hierarchy including node-local persistent memory, burst
buffers, and disk-based storage. For example, HDF5 [16], netCDF [26], and
Adaptable IO System (ADIOS) [27] are the most widely used high-performance I/O
libraries for HPC applications. However, these I/O libraries still suffer from
handling extremely large files (e.g. petabytes and beyond) due to inevitably
limited I/O bandwidths. Therefore, compression techniques are often adopted by
them [28]. For instance, H5Z-SZ [22] is a data filter for integrating SZ
compression into HDF5.

Figure 2 shows the abstraction of different layers in the I/O systems, where
compression is an individual layer in these systems. Specifically, compression
is normally performed between generating and storing the data  (e.g., H5Z-SZ).
It is worth noting that the compression tasks of all processes and the
parallel writing of compressed data to a shared file must be performed
sequentially; in other words, there must be a synchronization between these
two steps. This is because parallel writing of data from different processes
to a shared file requires the size/offset of each data partition, but
different processes may have vastly different data sizes after compression
(even with the same size before compression). Thus, the current compression-
write scheme cannot fully utilize the high compression ratio provided by lossy
compression, especially for large-scale HPC applications with multiple data
fields.

In this paper, considering HDF5 is well received by the HPC community as a
system supporting parallel I/O, we mainly focus our performance evaluation on
HDF5 without loss of generality.  In addition, to improve the performance and
productivity, a recent release of HDF5 [19] implements virtual object layer
(VOL) which can redirect I/O operations into VOL connector and allow
asynchronous I/O [24]. This feature enables an application to overlap I/O with
other operations such as compression.  Therefore, we can leverage this
capability to deeply integrate and overlap predictive lossy compression with
parallel write to improve the parallel write performance. Moreover, we focus
on the parallel write to a large shared file due to three main factors: (1) it
reduces scientists’ workload to manage multiple files for storage, post-hoc
analysis, and visualization; (2) it reduces the performance overhead of
opening/closing multiple files and the storage overhead of metadata for many
small files; and (3) partial processes (e.g., up to 4096 processes in [29]) of
a large-scale simulation with subfiling still writes to a shared file.

###  II-B Error-Bounded Lossy Compression

Lossy compression can compress data with extremely high compression ratio by
losing non-critical information in the reconstructed data. Two types of most
important metrics to evaluate the performance of lossy compression are: (1)
compression ratio, i.e., the ratio between original data size and compressed
data size, or bit-rate, i.e., the number of bits on average for each data
point on average (e.g., 32/64 for single/double-precision floating-point data
before compression); and (2) data distortion metrics such as peak signal-to-
noise ratio (PSNR) to measure the reconstructed data quality compared to the
original data. In recent years, a new generation of high accuracy lossy
compressors for scientific data have been proposed and developed for
scientific floating-point data, such as SZ [6, 5, 7] and ZFP [8]. These lossy
compressors provide parameters that allow users to control the loss of
information due to lossy compression precisely. Unlike traditional lossy
compressors such as JPEG [30] which are designed for images (in integers), SZ
and ZFP are designed to compress floating-point data and can provide a strict
error-controlling scheme based on user’s requirements. Generally, lossy
compressors provide multiple compression modes, such as error-bounding mode.
Error-bounding mode requires users to set an error type, such as point-wise
absolute error bound and point-wise relative error bound, and a bound value
(i.e., $10^{-3}$). The compressor ensures that differences between original
and reconstructed data do not exceed the error bound.

SZ is a prediction-based error-bounded lossy compressor for scientific data.
It has three main steps: (1) predict each data point’s value based on its
neighboring points by using an adaptive, best-fit prediction method; (2)
quantize the difference between the real value and predicted value based on
the user-set error bound; and (3) apply a customized Huffman coding and
lossless compression to achieve a high ratio.

Prior works have studied the impact of lossy compression on reconstructed data
quality and post-hoc analysis, providing guidelines on how to set the
compression configurations for certain applications [13, 31, 32, 33, 34, 35].
For example, a comprehensive framework was established to perform a systematic
analysis on compression configurations with a given dataset and provides the
best-fit solution that satisfies the post-hoc analysis [14]. Moreover, Jin et
al. [25] proposed a theoretical ratio-quality model to efficiently maximize
the compression ratio given the quality constraints of post-hoc analysis. Note
that similar to the previous work [36] on improving communication efficiency
via lossy compression, this work assumes that the compression-configuration is
set up by users based on their requirements of data quality, thus, the above
compression-configuration methods are orthogonal to our solution.

###  II-C Target I/O-Intensive HPC Applications

In this paper, we focus mainly on two I/O-intensive HPC applications—Nyx [18]
and VPIC [37], which have been used for many previous I/O studies [38, 39, 40,
41, 42, 43].

Nyx is an adaptive mesh, hydrodynamics code designed to model astrophysical
reacting flows on HPC systems [1, 18]. This code models dark matter as
discrete particles moving under the influence of gravity. The fluid in gas-
dynamics is modeled using a finite-volume methodology on an adaptive set of
3-D Eulerian grids/meshes. The mesh structure is used to evolve both the fluid
quantities and the particles via a particle-mesh method. For parallelization,
Nyx uses MPI for the long-range force calculation and architecture-specific
programming language for the short-range force algorithms, such as OpenMP and
CUDA. Nyx data uses multiple 3-D arrays to represent field information in grid
structure. According to prior studies [18, 44], it can run up to millions of
cores on leadership supercomputers such as Summit [45].

VPIC (vector particle-in-cell) is a large-scale plasma physics simulation that
can produce an unprecedented amount of data [37]. Collisionless magnetic
reconnection is an important mechanism that releases explosive amounts of
energy as field lines break and reconnect in plasmas. This reconnection also
plays an important role in a variety of astrophysical applications that
involve both hydrogen and electron-positron plasmas, including when the
Earth’s magnetosphere reacts to solar eruptions which can interfere satellite
communication. Simulation of magnetic reconnection with VPIC is inherently a
multi-scale problem that initiated in the small scale around individual
electrons but eventually leads to a large-scale reconfiguration of the
magnetic field. Recent simulations have revealed that electron kinetic physics
is important not only in triggering reconnection, but also in its subsequent
evolution [46]. This means plasma physics scientists find that they need to
model the detailed electron motion, and that modeling poses severe
computational challenges for 3D simulations of reconnection. A full-resolution
magnetosphere simulation is an exascale-class computing problem.

##  III Proposed Design Methodology

In this section, we present our proposed design that deeply integrates the
prediction-based lossy compression into HDF5 to significantly improve the
parallel-write performance.

###  III-A Design Overview

Figure 3: Overview of our proposed solution. Figure 4: Timeline of data
aggregation with 5 processes and 2 data fields.

Figure 3 shows an overview of our proposed framework. Specifically, we first
conduct the ratio and throughput prediction phase for all data partitions on
each process, which includes estimations of compression ratio, compression
throughput, and I/O throughput. Then, an all-gather communication is performed
to distribute the estimated compression ratio of each partition to all
processes. Note that the estimations of compression throughput and I/O
throughput are not distributed in this phase, since each process only requires
these estimations of its own data partition in order to perform the following
optimization. Next, each process computes the offset of its own data partition
of each field for parallel write based on the estimated compressed data size
with an extra space (a mechanism to handle overflow data under unexpected
prediction failure). Note that since each process gathers the same estimation
results from others, the computation of offsets is consistent across all
processes. After that, we perform an optimization on the order of compressing
different data fields in each process. Lastly, we overlap compressions and
writes based on the predicted write offsets and the optimized compression
order (the order of compressing different fields).

Compared to the non-compression solution that writes the entire data to a
shared file, our solution can significantly improve the overall write
performance due to the high compression ratio of lossy compression that
reduces I/O traffic. Compared to the previous lossy compression solution using
HDF5 filter, our solution can overlap compression and write that allows
independent and asynchronous write across processes. Figure 4 shows the
simplified timeline of parallel write from five processes to a shared file on
two data fields with four methods: (1) collective write without compression;
(2) collective write with compression via filter; (3) independent writes
overlapped with compressions; and (4) independent writes overlapped with
reordered compressions.

Specifically, the improvement from method (2) to method (3) is due to the
overlapping. From method (3) to method (4), the scheduling optimizer reordered
some compression tasks in each process to further improve the overlapping
efficiency. For example, the optimizer switches the compression order of the
two data fields in the last process, i.e., the data with smaller compressed
size are compressed later. The improvement of this optimization is more
significant when more data fields are compressed in parallel write as it has
higher chance to improve the overlapping efficiency. We will discuss more
observations and insights in detail in Section IV-D.

In this paper, we focus on accelerating parallel write due to two main
reasons: (1) HPC simulations are mostly write-oriented [47], and (2) parallel
I/O libraries usually have lower performance in write than in read [26, 48,
19, 49, 50].

In the following sections, we first introduce the prediction models for
prediction-based lossy compressors such as SZ. Next, we present the estimation
of I/O throughput based on our models and empirical studies. After that, we
propose our compression-write overlapping design with HDF5. Finally, we
propose our algorithm for compression scheduling.

###  III-B Compressor Throughput Estimation

The estimation on compression ratio in this work is based on recent research
[25], which proposed a ratio-quality model (i.e., compression ratio and
reconstructed data quality model) for prediction-based lossy compressor. It
developed algorithms to predict the compression ratio based on a newly
designed sampling strategy without performing compression, significantly
reducing the time overhead of the compression ratio prediction to less than
$10\%$ (in relative to the compression time). We leverage this idea to enable
the overlapping design with small prediction overhead in our design.

On the other hand, estimating the compression throughput is also essential for
compression order optimization. Generally, higher compression ratio results in
higher compression throughput [13], due to the shorter time for building
Huffman tree and encoding with smaller tree. However, there is no prior work
for an accurate compression-throughput estimation.

Figure 5: Single-core compression throughput with different bit-rates on a Nyx
and a RTM datasets. Evaluated on Bebop.

In this work, we propose to estimate compression throughput based on the
predicted compression ratio for each data partition. Figure 5 shows our
empirical evaluation of compression throughput and compression ratio with
multiple data fields and data types across different error bounds. Note that
we use bit-rate (the average bits used to represent each value) on the x-axis
to better illustrate the size of compressed data. For example, a bit-rate of 4
is equivalent to a compression ratio of $8\times$ when compressing single-
precision floating-point data (32 bits/value). The evaluation is performed on
a single CPU core with SZ lossy compressor, since our target method (4) is to
compress the data in each process/core independently before writing the
compressed data to the shared file. The figure illustrates that (1) the
maximum and minimum compression throughput are similarly bounded across
different data samples (i.e., about 120$\sim$250 MB/s); and (2) the bitrate-
throughput curve for each data sample is highly consistent.

Figure 6: Minimum and maximum compression throughput of a given data partition
based on 30 samples from Baryon density, dark matter density, temperature and
velocity x data fields in a Nyx dataset.

The reasons behind the first observation are: (1) For extremely high error
bounds, despite the small size of the Huffman tree (negligible tree building
time), the prediction and encoding process still passes each point, which
provides an upper bound on the throughput; (2) For extremely low error bounds,
since SZ limits the maximum size of the Huffman tree (the maximum number of
quantization codes), if one point is unpredictable with the maximum number of
quantization codes, its original value will be saved directly, which provides
a lower bound on the throughput. Figure 6 shows the maximum and minimum
compression throughput evaluated on multiple data samples size of $67.1$ MB
from different data fields and types on the same experiment platform.

The reason behind the second observation is that for different data samples,
the quantization codes after prediction are all centrally distributed, and the
overall bit-rates are similar under the same Huffman encoding efficiency,
where the sizes of the Huffman trees are also similar, resulting in stable
tree building time and encoding time.

Based on Figure 5, we use a power function to predict compression throughput

$\displaystyle\small\textstyle T_{comp}$ $\displaystyle=D/S$ (1)
$\displaystyle=(B_{ori}\times n)/(((C_{max}-C_{min})\times
3^{-a})B^{a}+C_{min}),$

where $T_{comp}$ is the estimated compression time; $D$ is the size of the
original data; $S$ is the compression throughput; $B$ is the compressed bit-
rate (i.e., the average bits used to save each value, e.g., a compression
ratio of 16 from 32 bits floating point values results in $B=2$); $B_{ori}$ is
the original bit-rate (i.e., 32 for single-precision floating point values);
$n$ is the number of points in the current data partition; $C_{min}$ and
$C_{max}$ are the minimum and maximum compression throughputs, respectively;
$a$ is a hyper-parameter to describe the shape of the power function (i.e., a
value smaller than 0; the lower the value, the more curved the function). The
compression-throughput estimation is a power function that takes min-max into
account. Note that the number 3 in the equation is based on our experiment
that yields the best result.

This prediction model can be adapted for different machines, where we perform
lossy compression on a sample of dataset offline to evaluate the minimum and
maximum compression throughput and map it based on Equation (1). We show the
accuracy of our proposed model in Section IV-B.

###  III-C Write Time Estimation

The write time estimation is needed when optimizing the compression order,
since the compressed data with larger estimated write time tends to be
processed earlier to maximize the I/O traffic occupancy. However, compared
with the estimation on compression ratio and throughput, the accuracy required
for write-time estimation is considerably lower. This is because one
inaccurate write-time prediction makes a linear shift to all other write-time
predictions within this process,  so the relative time spent for each write
remains unchanged compared to the entire time,  and the optimization would not
be affected (but only increase or decrease the actual write time). As a
result, our goal is not to provide a highly accurate write-time estimation for
each data partition, but to provide a capability to estimate the relative
write time across different data sizes.

Figure 7: Independent write I/O throughput per process with different data
sizes per process. Evaluated on 128 processes.

Figure 7 shows the parallel-write throughput per process with different data
partition sizes. We can observe that the average throughput first increases as
the data size increases and stabilizes after the data size reaches a certain
point. Similar to the compression-throughput estimation, the write-time
estimation is also based on the predicted bit-rate. Note that the compression
throughput is based on the uncompressed data size. Considering that usually
each process handles a similar amount of data during the simulation, the
variance of the compression-time difference is limited to $T_{max}/T_{min}$.
While the write throughput is based on the compressed data size, which can
easily vary over $10\times$ across data partitions according to our
evaluation. This means the write time is mainly dependent on the compression
ratio assuming the write throughput is relatively stable across different
processes without I/O congestion. Thus, we estimate the write time as

$\displaystyle\small\textstyle T_{write}=(B\times n)/C_{thr},$ (2)

where $T_{write}$ is the write time, $B$ is the compressed bit-rate; $n$ is
the number of points in the data partition; and $C_{thr}$ is the stable write
throughput based on an empirical evaluation.

###  III-D Overlapping Compression and Write

Figure 8: Overflow data handling with preserved extra space.

When writing data from different processes to a shared file, the well-known
I/O libraries such as HDF5, MPI-IO, and ADIOS require users to provide an
offset for each data partition. For the original non-compression parallel
write, the offset of each partition can easily be computed based on the data
size. However, with compression, the size of each compressed partition can
vary drastically due to diverse data compressibility and different compression
configurations. Prior to our work, it was impossible to pre-compute the offset
before the compression, thus, it must sequentially process the collective
write after compression. Thanks to the compression-ratio estimation for each
partition [25], we can estimate the offset even before actual compression and
allow overlapping between compression and write via HDF5’s asynchronous I/O.

We note that although the compression-ratio model can provide high accuracy
for prediction-based lossy compression, the prediction does not have a
guarantee bound. This means that we must handle the situation when the actual
compressed data size is larger than the predicted size, since it is
impractical to simply recompute and update the offsets of all following
partitions every time when the prediction failure happens as it would
introduce a large amount of communication overhead. To this end, we propose to
reserve an particular extra space for each prediction to handle the
possibility of compressed data overflow. Moreover, we also propose an overflow
data handling scheme to redirect and store the exceeded portion beyond the
extra space. Figure 8 shows our proposed overflow handling strategy with an
extra data space.

Figure 9: Trade-off between performance overhead and compression size
overhead. Shown the empirical average result based on Nyx and VPIC datasets on
512 processes.

As the compression-ratio estimation error is not bounded, we cannot simply set
a maximum extra space to guarantee no data overflow in practice. Thus, we take
the extra space ratio as one of the tunable parameters in our framework. For
example, setting this ratio $R_{space}$ to 1.5 means that we over preserve 50%
storage space (in relative to the predicted compressed size) when computing
the offset of each partition.

In our evaluation, we find that the compression-ratio model [51] performs
poorly under extremely high compression ratio because of three factors: (1)
when the compression ratio is high, the Huffman encoder can only provide a
maximum compression ratio of $32\times$ (when the original data is single-
precision) and relies on the following lossless compression to further reduce
the encoded data; (2) the compression-ratio model is based on run-length
encoding to analyze the lossless encoding efficiency, which naturally features
lower estimation accuracy compared to the efficiency estimation of Huffman
encoding; and (3) when the compression ratio is low, the encoded data stream
from the Huffman encoder is highly random and is hardly further reduced by
lossless compressors. Thus, when handling the data partition with compression
ratio higher than 32 (bit-rate lower than 1), we set $r_{space}$ to

$\displaystyle\small\textstyle r_{space}=$
$\displaystyle\min(2,1+(R_{space}-1)\times 4),$ $\displaystyle where\hskip
14.22636ptr_{comp}>32.$ (3)

where $r_{comp}$ is the compression ratio.  On one hand, the higher the extra
space ratio is, the larger the storage size of the compressed shared file will
be, due to the larger extra space that is likely to be wasted.  On the other
hand, increasing $r_{space}$ may improve the overall performance due to less
data overflow processing overhead, and vice versa. Overall, above a certain
threshold (e.g., $r_{space}>1.1\times$ in our evaluation), increasing
$r_{space}$ becomes a trade-off between the write performance and the overall
compression efficiency.

More specifically, below this threshold, decreasing $r_{space}$ results in a
significant write performance drop due to a large number of compressed data
overflows. For example, when $r_{space}=1.1\times$, $32.4\%$ data partitions
suffer from data overflow, causing an extra $65.6\%$ time overhead. Based on
our observations from Figure 14, the trade-off between write performance and
compression efficiency is relatively similar across different fields and
datasets. We will demonstrate this with more details in Section IV-C.

As a result, we propose to weight the write performance and compression
efficiency overhead and provide a tunable ratio between them for users. For
any given weights, we provide the extra space ratio based on our proposed
mapping. Figure 9 shows our mapping solution for the extra space ratio (a more
detailed evaluation illustrating the accuracy of this mapping will be shown in
Section IV-C). Note that we only support $r_{space}$ between $[1.1,1.43]$ due
to (1) an extremely high time overhead below 1.1, and (2) a low efficiency of
trading storage for performance after 1.43. We set the default extra space
ratio to 1.25 in this work.

After estimating the offset of each data partition, we must store this offset
information as metadata for the decompression purpose. Compared to the
compressed size of the entire dataset, this metadata size is totally
negligible. For example, when writing a total size of 2.5 TB Nyx dataset from
4,096 processes with SZ, there is only 295 KB metadata to store the offsets
for 9 fields, while the compressed data size is 210 GB.

For handling the compressed data overflow, our goal is to append the exceeded
data to the end of the shared file. Specifically, we continue to optimize and
write the maximum amount of data into the preserved area of the shared file.
Then, after all processes finish writing, we initiate an all-gather operation
across all processes to distribute the size information of overflow data.
After that, some processes calculate the offset of their own excess and write
it to the end of the shared file independently. Note that since only a small
fraction of processes may need to process/save a small amount of overflow
data, the overall time overhead for processing overflow data is small. A more
detailed evaluation will be presented in Section IV-D.

###  III-E Compression Order Optimization

Algorithm 1 Compression Order Optimization

Notation: data fields in current process: $\ell$; compression queue: $Q$;
compression queue after insert and additional data: $Q^{\circ}$; possible
insert locations in a queue: $\beta$; time to compress: $t_{c}$; time to
write:$t_{w}$; predicted compression time: $P_{c}(\ell)$; predicted write
time: $P_{w}(\ell)$

Global: $P_{c}(\ell),P_{w}(\ell)$

1procedure Time($q$)

2 $t_{c},t_{w}\leftarrow 0$

3 for $\ell\leftarrow$ data fields in $q$ do

4 $t_{c}\leftarrow t_{c}+P_{c}(\ell)$

5 $t_{w}\leftarrow P_{w}(\ell)+\max(t_{c},t_{w})$

6 end for

7 return $t_{w}$

8end procedure

9

10procedure SchedulingOptimizator

11 for $\ell\leftarrow$ data fields in current process do

12 for $\beta\leftarrow$ all possible insert location do

13 $Q^{\circ}\leftarrow$ insert $\ell$ to $\beta$

14 if $\textsc{Time}(Q^{\circ})<\textsc{Time}(Q)$ or first $\beta$ then

15 $Q\leftarrow Q^{\circ}$

16 end if

17 end for

18 end for

19 return $Q$

20end procedure

As shown in Figure 4, when writing data to the shared file, overlapping the
compression and I/O can provide a significant performance benefit over the
previous collective-write solution with lossy compression. However, with
multiple data fields in each process, compressing them sequentially in the
original order is not the optimal solution. In fact, since we use the modeling
approach to predict the compressed size, the compression time, and the write
time of each data partition, we can reorder the compression tasks for
different fields to maximize the overlapping without any penalty.

Algorithm 1 shows the pseudocode for optimizing the compression order in each
process. The proposed method is based on the observation that the total
compression time is theoretically fixed regardless of the compression order.
Our optimization focuses on the dependencies and timing of launching write
operations for each compressed data to minimize timeouts compared to
compression. The time complexity of the proposed algorithm is $O(n^{2})$,
whereas the time complexity of our compression is $O(N)$. Considering that $N$
(i.e., the number of values) in one data partition is significantly larger
than $n$ (i.e., the number of data fields), the optimization overhead is
almost negligible compared to the actual compression and write time. For
example, this optimization overhead is only 0.17% of the compression time
under an extreme condition where $N$ is small at 32768 ($32\times 32\times
32$) but $n$ is very large at 100.

Figure 10: An example of extremely unbalanced compression time and write time,
limiting the benefit from our reordering.

Based on our algorithm design, we can expect the optimization to bring benefit
when the balance between compression time and I/O time is relatively stable.
However, we notice that in the extreme scenarios, the compression time and the
write time could be very unbalanced, which can diminish the benefit of our
reordering optimization, as shown in Figure 10. Specifically, there are two
cases: (a) when the write time is significantly longer than the compression
time, or (b) when the compression time is significantly longer than the write
time. In both cases, there is no significant room to improve performance
because of the limited overlap between compression and I/O (more details will
be discussed in Section IV-D).

In addition, we note that when the number of data fields is relatively large,
our optimization can bring greater benefit. This is because the overall
performance is dependent on the worst process with the longest time among all
the processes due to independent asynchronous writes.

##  IV Experimental Evaluation

In this section, we present the evaluation results of our proposed framework
for accelerating parallel write. We first evaluate the accuracy of our
prediction for compression time and write time. Next, we evaluate the extra
space setup under different compression ratios and scales. Finally, we perform
performance evaluation and scaling study of our approach and compare it with
the original solution without compression and the solution using using the
H5Z-SZ filter [22].

###  IV-A Evaluation Setup

TABLE I: Details of Tested Datasets

Name | Description | Scale | Size
---|---|---|---
nyx [18] | Cosmology simulation |  4096$\times$4096$\times$4096 | 2.47 TB
2048$\times$2048$\times$2048 | 206.15 GB
1024$\times$1024$\times$1024 | 25.76 GB
512$\times$512$\times$512 | 3.22 GB
VPIC [52] | Particle simulation | 161,297,451,573 | 4.62 TB

We rigorously implement our approach with HDF5 [19] and SZ3 [32] (a
modularized prediction-based lossy compressor). We conduct our experiments on
two HPC systems: (1) the Summit supercomputer [45] at Oak Ridge National
Laboratory, each node of which is equipped with two IBM POWER9 processors with
42 physical cores and 512 GB DDR4 memory, and (2) the Bebop cluster [53] at
Argonne National Laboratory, each node of which is equipped with two 18-core
Intel Xeon E5-2695v4 CPUs and 128 GB DDR4 memory.

We use different scales of Nyx and VPIC datasets in our evaluation. Table I
shows the details of our tested datasets. According to the previous work [13,
31], using the absolute error bounds of $(0.2,0.4,1e+3,2e+5,2e+5,2e+5)$ for
compressing the six Nyx data fields (i.e., baryon density, dark matter
density, temperature, velocity x, velocity y, velocity z) can satisfy the
post-hoc analysis quality  with an average PSNR (peak signal-to-noise ratio)
at 78.6 dB, resulting in a compression ratio of $\sim$$16\times$.  The
$4096\times 4096\times 4096$ Nyx dataset has three additional fields:
particle_vx, particle_vy and particle_vz.  Similarly, we compress them with a
compression ratio of $16\times$ to satisfy the post-hoc analysis quality. For
the VPIC dataset, we also compress the 8 data fields with a compression ratio
of $13.8\times$, which is suggested by the application developers according to
their post-hoc analysis.  Note that our solution can also be applied to other
scientific datasets with similar performance improvement expected. This is
because: (1) our solution provides the same reconstructed data quality
compared to previous solutions, where many studies [13, 31, 32, 33, 34, 35,
14] show that lossy compression has a high data reduction capability on a
variety of applications; (2) a previous study shows that the accuracy of the
compression-ratio estimation used in our solution is consistently above 90%
across tens of benchmark datasets [25], implying that the proportion of data
overflows due to inaccurately estimated compression ratios is relatively
stable in most scientific datasets; (3) the accuracy of the compression-
throughput estimation is consistently high across different datasets and
fields, as shown in Figure 5, 11 and 12; and (4) the accuracy of
I/O-throughput estimation only relies on the estimated compressed data size as
discussed in (2).

###  IV-B Accuracy of Compression and I/O Throughput Estimation

Figure 11: Accuracy of our compression-time estimation on $512^{3}$ Nyx data
samples (red line is predicted time; black dots are actual time). Figure 12:
Accuracy of our compression-time estimation on $1024^{3}$ Nyx data samples.
Red line is predicted time; black dots are actual time. Offline parameter is
based on the baryon density of $512^{3}$ Nyx dataset.

First, we evaluate the accuracy of our proposed estimation on compression
time. We perform the evaluation on the $512^{3}$ Nyx dataset only using the
baryon density field and conduct the offline compression with the relative
error bound ranging between [1e-1, 1e-8]. Then we calculate $C_{min}$,
$C_{max}$ and $a$ in Equation (1) based on the offline results, so that we can
use Equation (1) to predict the compression time of any given data (i.e.,
$101.7$, $240.6$, and $-1.716$, respectively, in this case). Next, we
distribute the Nyx dataset to 64 processes (each process has a $128\times
128\times 128$ data partition) and perform the write operation. For each field
of data partition in each process, we predict the compression time using
Equation (1). Figure 11 shows the actual compression time versus the predicted
compression time.

Furthermore, we extend our evaluation to the $1024^{3}$ Nyx dataset and 512
processes using the same $C_{min}$, $C_{max}$, and $a$ obtained from the
$512^{3}$ dataset. Figure 12 shows the actual compression time versus the
predicted compression time. We can observe from both figures that our
compression-time estimation has high accuracy even though the offline
experiment only uses baryon density as input. This is consistent with the
result shown in Figure 5, where different data fields and datasets have
similar compression throughputs.

Figure 13: Accuracy of our write time estimation on $1024^{3}$ Nyx data
samples. Red line is predicted time; black dots are actual time.

(a) Nyx

(b) VPIC

Figure 14: Trade-off between performance overhead and storage overhead based
on different extra space ratios on Nyx dataset (6 data fields) and VPIC
dataset (7 data fields) on both Bebop and Summit with 512 processes. The
target compressed bit-rate is 2 bits/value.

Similarly, we evaluate the accuracy of our proposed write-time estimation. We
perform the offline evaluation by writing data from 128 processes into a
shared file multiple times. The size of the data per process is set to 5 MBs,
10 MBs, 20 MBs, 50 MBs, or 100 MBs. We measure the average write throughput.
Based on our observation, further increasing or decreasing the scale would not
significantly affect the average I/O throughput per process; therefore, we
only perform the offline evaluation on one scale with different data sizes to
reduce the offline evaluation overhead. Next, we distribute the $1024^{3}$ Nyx
dataset to 64 processes, perform the write operation, and measure the time of
independent write for each process along with the bit-rate after compression.

Figure 13 shows the comparison between the estimated write time and the actual
write time. Note that the accuracy of low bit-rate is slightly lower than that
of high bit-rate. This is because the compressed data size is very small
(i.e., 0.94 MB in this case) under low bit-rate, which result in a
significantly low write throughput. However, since we optimize the compression
order based only on the absolute write time, such a small amount of write-time
prediction error has hardly any effects on optimization decisions.

###  IV-C Evaluation on Extra Space Ratio

As mentioned in Section III-D, the extra space ratio is one of the most
important parameters that balance the overall write performance and storage
overhead in our framework. We then evaluate this performance-storage trade-off
on both the Nyx and VPIC datasets. For the Nyx dataset, we use the $1024^{3}$
dataset and distribute it to 512 processes for parallel write. For the VPIC
dataset, we downscale the original dataset by sampling 10% data points and
distribute it to 512 processes.

Figure 14 shows the correlation between the write-performance overhead and the
storage overhead with the two datasets on both Bebop and Summit. Note that the
write-performance overhead is compared with the write time without handling
data overflow (excluding the compression time), and the storage overhead is
compared with the ideal compressed data size (without the extra space). We can
observe that for each dataset on a given system, the trade-off curve is highly
similar across different data fields, since the accuracy of the compression-
ratio model is relatively stable on the same dataset. We also note that even
though the write-performance overhead and storage overhead with the same extra
space ratio are different on the two datasets, the trade-off between the two
overheads is very similar. This is because the lower upper bound of the
compression-ratio estimation results in larger number of processes to hold
overflow data and hence higher performance overhead, which also enlarges the
total amount of overflow data and increases the overall storage overhead.

Figure 15: Evaluation on the consistency of the storage and performance
overheads using the same extra space ratio of 1.25 with 512 processes on
Summit. Red shift stands for different time-steps (higher values means earlier
time in the simulation). Figure 16: Performance comparison among our solution
(overlapping and reordering), original non-compression solution, and previous
compression-write solution on $4096^{3}$ Nyx dataset with 512 processes.

(a) Nyx with different compression ratio

(b) VPIC with different compression ratio

(c) Nyx with different scale

(d) VPIC with different scale

Figure 17: Performance improvement of overall parallel-write with our proposed
solution compared to the previous write solution with H5Z-SZ on both Nyx and
VPIC datasets. Dashed red line is the baseline of HDF5 without compression.
(a) and (b) are evaluated with 512 processes on Summit. (c) and (d) are
evaluated with a target bit-rate of 2.

Regarding the results across the two systems, the difference of the
performance overhead is mainly due to the higher I/O bandwidth of Summit over
Bebop, which reduces the overall write time and enlarges the relative
performance overhead. However, the performance-storage tread-off is still
similar across different systems, which means that we can use our offline
study to guide the online decision making for the trade-off between
performance and storage. Figure 14 shows the averages of the performance and
storage overhead in different extra space ratios across the four dataset-
system setups (as mentioned in Section III-D). By default, we use the extra
space ratio of 1.25 to minimize performance overhead while keeping low storage
cost.  In addition, we apply the extra space ratio of 1.25 with multiple time-
steps in a series of Nyx datasets, as shown in Figure 15. It shows that our
solution has consistent storage and performance overhead across time-steps.
Another alternative approach is to sample the dataset from the first time-step
and provide users a more accurate estimation on the trade-off between
performance and storage. We will design a more user-friendly parameter-tuning
mechanism in future.

###  IV-D Overall Performance Improvement

Next, we compare the performance of our proposed solution with the non-
compression solution and the previous compression-filter solution. For the
non-compression solution, the distributed data is written to the shared file
without any compression. We implement it using HDF5 with independent write,
which can significantly increase the performance over collective write [19].
For the previous compression-filter solution, the compression process is
considered as a filter, and the parallel write starts only when every process
finishes its compression on all data partitions and the compressed data sizes
are known to all processes to compute the write offsets. We implement it using
HDF5 with collective write, since compression filters do not support
independent write. Note that HDF5 only provides additional features such as
asynchronous I/O and does not change the collective write.

(a) Nyx with different compression ratio

(b) VPIC with different compression ratio

(c) Nyx with different scale

(d) VPIC with different scale

Figure 18: Performance improvement (overall) and storage overhead of our
solution compared to the previous solution on both Nyx and VPIC datasets. (a)
and (b) are evaluated with 512 processes om Summit. (c) and (d) are evaluated
with a target bit-rate of 2.

Figure 16 shows the comparison of the performance breakdown between our
proposed solution and the two existing solutions. We implement our solution
both with and without the compression reordering, referred to as “overlapping”
and “reordering” in the figure, respectively. Note that the write time bar
shown in the overlapping solution is measured by the time between the end of
the slowest compression and the end of the writing process rather than the
entire write time (due to overlapping). The original/ideal compression ratio
(without the extra space) in this experiment is $17.94\times$, while the
actual compression ratio considering the extra space (i.e., the extra space
ratio of $1.25$) is $14.13\times$.

We can observe that the collective write solution with lossy compression still
outperforms the independent write without compression by $1.87\times$ due to
the high compression ratio and reduced data size. Compared to the previous
compression-write solution, our overlapping solution further improves the
overall write performance by $1.79\times$ due to the asynchronous independent
write. Note that the compression times spent in the two solutions are similar,
indicating that our framework improves the writing efficiency rather than the
compression throughput. Furthermore, with compression order optimization, we
further improve the performance by $1.30\times$ in addition to our overlapping
optimization. The extra write time (gray bar) is significantly reduced due to
the high overlapping efficiency. Overall, comparing our optimized solution to
the non-compression write, we improve the write performance by $4.46\times$
with the compression ratio of $14.1\times$; compared to the previous write
with H5Z-SZ with the compression ratio of $17.9\times$, our optimized solution
improves the write performance by $2.91\times$ with a $26\%$ storage overhead.
Note that this overhead is compared with the compressed data size, if compared
with the original data size, the storage overhead is equivalent to $1.5\%$.

Furthermore, we evaluate our framework with different overall compression
ratios and different scales to illustrate the efficiency of compression
reordering on both the Nyx and VPIC datasets. Figure 17(a) and Figure 17(b)
show that the compression reordering optimization has a poor improvement over
the overlapping-only solution under extremely high and low compression ratio,
the corresponding performance improvement and storage overhead are shown in
Figure 18(a) and Figure 18(b). For extremely high compression ratio, the
compressed data size in each process is tiny, so the write time is
significantly smaller than the compression time, as shown in Figure 17(a) and
Figure 17(b) (the short gray bar versus the orange bar). In this case,
reordering the compression tasks provides little performance benefit, since
the overall compression time cannot be further improved and the potential of
reducing the extra write time is negligible. For extremely low compression
ratio, the compressed data size in each process is large, so the overall
compression time is significantly shorter than the write time, where the
overlapping efficiency is likely to be sufficiently high even without
optimization (low potential to optimize). Both datasets have a similar
compression ratio ($14.1\times$ and $10.9\times$ for Nyx and VPIC,
respectively) that benefits the compression reordering optimization. Note that
although the higher compression ratio almost always indicates the better write
performance, it also means lower data quality with information loss that may
potentially hurt post-hoc analysis. Under the extreme cases when the
compression ratio is very small, the previous write solution with H5Z-SZ shows
even worse performance than the non-compression write.

Finally, we perform the scaling study of our proposed framework. Specifically,
we conduct a weak scaling study that keeps the same data partition size in
each process, i.e., $256^{3}$ and $39,379,260$ for the Nyx and the VPIC
datasets, respectively. Figure 17(c) and Figure 17(d) illustrate that the
compression order optimization over the overlapping solution is relatively
stable on the evaluated scales at $(256,512,1024,2048,4096)$, the
corresponding performance improvement and storage overhead are shown in Figure
18(c) and Figure 18(d). This is because the overall compression time and write
time in each process are mostly stable through different scales. However, a
larger scale slightly benefits our solution because the independent,
asynchronous write typically provides better scalability compared to the
collective write used by the previous compression-write solution [19]. Note
that the optimization time and the overflow time increase with the scale,
because even though the prediction time is stable, larger scale introduces
longer communication time for the all-gather operation.

In conclusion, under the circumstances of satisfying the user-set compression
configuration, our parallel write with lossy compression is scalable and high-
performance, especially when (1) the number of data fields is relatively
large; (2) the overall compression ratio is preferably between $10\times$ and
$20\times$ (with balanced compression time and write time); and (3) the data
amount in each process is not too small (deserving compression). In addition,
users can fine-tune the extra space ratio to make a good performance-storage
trade-off.

##  V Conclusion and Future Work

In this paper, we propose to integrate predictive lossy compression deeply
with HDF5 to significantly improve the parallel-write performance for large-
scale scientific simulations. We introduce a newly designed prediction method
to estimate the compression and write time for each process and use this
information to pre-compute the write offset before the actual compression.
Furthermore, we introduce an extra space design to handle the uncertainty of
the prediction and an optimization algorithm to reorder the compression tasks
in each process. We evaluate our proposed solution on both Bebop cluster and
Summit supercomputer with up to 4,096 cores. The evaluation shows that our
solution can provide a $4.46\times$ performance improvement with a
$14.1\times$ compression ratio compared to the original parallel write without
compression, and provide a $2.91\times$ performance improvement compared to
the previous parallel write with the SZ compression filter with only $20\%$
compression ratio degradation.

In the future, we will extend our solution to other parallel I/O libraries
such as ADIOS [27] and support more lossy compressors such as ZFP [8].
Moreover, we plan to evaluate our solution on more real-world HPC datasets.

## Acknowledgments

This research was supported by the Exascale Computing Project (ECP), Project
Number: 17-SC-20-SC, a collaborative effort of two DOE organizations—the
Office of Science and the National Nuclear Security Administration,
responsible for the planning and preparation of a capable exascale ecosystem,
including software, applications, hardware, advanced system engineering and
early testbed platforms, to support the nation’s exascale computing
imperative. The material was supported by the U.S. Department of Energy,
Office of Science, Advanced Scientific Computing Research (ASCR), under
contracts DE-AC02-06CH11357 and DE-AC02-05CH11231. This work was also
supported by the National Science Foundation under Grants OAC-2003709,
OAC-2042084, OAC-2104023, and OAC-2104024. We gratefully acknowledge the
computing resources provided on Bebop cluster at Argonne and Summit
supercomputer at Oak Ridge.

## References

  * [1] A. S. Almgren, J. B. Bell, M. J. Lijewski, Z. Lukić, and E. Van Andel, “Nyx: A massively parallel amr code for computational cosmology,” _The Astrophysical Journal_ , vol. 765, no. 1, p. 39, 2013.
  * [2] L. Wan, M. Wolf, F. Wang, J. Y. Choi, G. Ostrouchov, and S. Klasky, “Comprehensive measurement and analysis of the user-perceived i/o performance in a production leadership-class storage system,” in _2017 IEEE 37th International Conference on Distributed Computing Systems (ICDCS)_. IEEE, 2017, pp. 1022–1031.
  * [3] ——, “Analysis and modeling of the end-to-end i/o performance on olcf’s titan supercomputer,” in _2017 IEEE 19th International Conference on High Performance Computing and Communications; IEEE 15th International Conference on Smart City; IEEE 3rd International Conference on Data Science and Systems (HPCC/SmartCity/DSS)_. IEEE, 2017, pp. 1–9.
  * [4] F. Cappello, S. Di, S. Li, X. Liang, A. M. Gok, D. Tao, C. H. Yoon, X.-C. Wu, Y. Alexeev, and F. T. Chong, “Use cases of lossy compression for floating-point data in scientific data sets,” _The International Journal of High Performance Computing Applications_ , 2019.
  * [5] D. Tao, S. Di, Z. Chen, and F. Cappello, “Significantly improving lossy compression for scientific data sets based on multidimensional prediction and error-controlled quantization,” in _2017 IEEE International Parallel and Distributed Processing Symposium_. IEEE, 2017, pp. 1129–1139.
  * [6] S. Di and F. Cappello, “Fast error-bounded lossy HPC data compression with SZ,” in _2016 IEEE International Parallel and Distributed Processing Symposium_. IEEE, 2016, pp. 730–739.
  * [7] X. Liang, S. Di, D. Tao, S. Li, S. Li, H. Guo, Z. Chen, and F. Cappello, “Error-controlled lossy compression optimized for high compression ratios of scientific datasets,” 2018.
  * [8] P. Lindstrom, “Fixed-rate compressed floating-point arrays,” _IEEE Transactions on Visualization and Computer Graphics_ , vol. 20, no. 12, pp. 2674–2683, 2014.
  * [9] M. Ainsworth, O. Tugluk, B. Whitney, and S. Klasky, “Multilevel techniques for compression and reduction of scientific data—the univariate case,” _Computing and Visualization in Science_ , vol. 19, no. 5–6, pp. 65–76, 2018\.
  * [10] T. Lu, Q. Liu, X. He, H. Luo, E. Suchyta, J. Choi, N. Podhorszki, S. Klasky, M. Wolf, T. Liu, and Z. Qiao, “Understanding and modeling lossy compression schemes on HPC scientific data,” in _2018 IEEE International Parallel and Distributed Processing Symposium_. IEEE, 2018, pp. 348–357.
  * [11] H. Luo, D. Huang, Q. Liu, Z. Qiao, H. Jiang, J. Bi, H. Yuan, M. Zhou, J. Wang, and Z. Qin, “Identifying latent reduced models to precondition lossy compression,” in _2019 IEEE International Parallel and Distributed Processing Symposium_. IEEE, 2019.
  * [12] D. Tao, S. Di, X. Liang, Z. Chen, and F. Cappello, “Optimizing lossy compression rate-distortion from automatic online selection between sz and zfp,” _IEEE Transactions on Parallel and Distributed Systems_ , 2019.
  * [13] S. Jin, P. Grosset, C. M. Biwer, J. Pulido, J. Tian, D. Tao, and J. Ahrens, “Understanding gpu-based lossy compression for extreme-scale cosmological simulations,” _arXiv preprint arXiv:2004.00224_ , 2020.
  * [14] P. Grosset, C. Biwer, J. Pulido, A. Mohan, A. Biswas, J. Patchett, T. Turton, D. Rogers, D. Livescu, and J. Ahrens, “Foresight: analysis that matters for data reduction,” in _2020 SC20: International Conference for High Performance Computing, Networking, Storage and Analysis (SC)_. IEEE Computer Society, 2020, pp. 1171–1185.
  * [15] S. W. Son, Z. Chen, W. Hendrix, A. Agrawal, W.-k. Liao, and A. Choudhary, “Data compression for the exascale computing era-survey,” _Supercomputing Frontiers and Innovations_ , vol. 1, no. 2, pp. 76–88, 2014\.
  * [16] The HDF Group. (2000-2010) Hierarchical data format version 5. [Online]. Available: http://www.hdfgroup.org/HDF5
  * [17] M. Folk, G. Heber, Q. Koziol, E. Pourmal, and D. Robinson, “An overview of the hdf5 technology suite and its applications,” in _Proceedings of the EDBT/ICDT 2011 Workshop on Array Databases_ , 2011, pp. 36–47.
  * [18] Nyx, https://github.com/AMReX-Astro/Nyx, 2021.
  * [19] S. Byna, M. S. Breitenfeld, B. Dong, Q. Koziol, E. Pourmal, D. Robinson, J. Soumagne, H. Tang, V. Vishwanath, and R. Warren, “Exahdf5: delivering efficient parallel i/o on exascale computing systems,” _Journal of Computer Science and Technology_ , vol. 35, no. 1, pp. 145–160, 2020.
  * [20] S. Pokhrel, M. Rodriguez, A. Samimi, G. Heber, and J. J. Simpson, “Parallel i/o for 3-d global fdtd earth–ionosphere waveguide models at resolutions on the order of~ 1 km and higher using hdf5,” _IEEE Transactions on Antennas and Propagation_ , vol. 66, no. 7, pp. 3548–3555, 2018.
  * [21] The HDF Group. (2000-2010) Hierarchical data format version 5, Filter. [Online]. Available: https://support.hdfgroup.org/HDF5/doc/H5.user/Filters.html
  * [22] H5Z-SZ, https://github.com/disheng222/H5Z-SZ.
  * [23] T. H. team, “Parallel compression improvements in HDF5 1.13.1,” https://www.hdfgroup.org/2022/03/parallel-compression-improvements-in-hdf5-1-13-1/, 2022, [Online; accessed 01-April-2022].
  * [24] H. Tang, Q. Koziol, J. Ravi, and S. Byna, “Transparent Asynchronous Parallel I/O Using Background Threads,” _IEEE Transactions on Parallel and Distributed Systems_ , vol. 33, no. 4, pp. 891–902, 2022.
  * [25] S. Jin, S. Di, S. Byna, D. Tao, and F. Cappello, “Improving prediction-based lossy compression dramatically via ratio-quality modeling,” _arXiv preprint arXiv:2111.09815_ , 2021.
  * [26] J. Li, W.-k. Liao, A. Choudhary, R. Ross, R. Thakur, W. Gropp, R. Latham, A. Siegel, B. Gallagher, and M. Zingale, “Parallel netcdf: A high-performance scientific i/o interface,” in _SC’03: Proceedings of the 2003 ACM/IEEE Conference on Supercomputing_. IEEE, 2003, pp. 39–39.
  * [27] W. F. Godoy, N. Podhorszki, R. Wang, C. Atkins, G. Eisenhauer, J. Gu, P. Davis, J. Choi, K. Germaschewski, K. Huck _et al._ , “Adios 2: The adaptable input output system. a framework for high-performance data management,” _SoftwareX_ , vol. 12, p. 100561, 2020.
  * [28] H. Tang, S. Byna, N. A. Petersson, and D. McCallen, “Tuning parallel data compression and i/o for large-scale earthquake simulation,” in _2021 IEEE International Conference on Big Data (Big Data)_. IEEE, 2021, pp. 2992–2997.
  * [29] S. Byna, M. Chaarawi, Q. Koziol, J. Mainzer, and F. Willmore, “Tuning hdf5 subfiling performance on parallel file systems,” 2021.
  * [30] G. K. Wallace, “The JPEG still picture compression standard,” _IEEE Transactions on Consumer Electronics_ , vol. 38, no. 1, pp. xviii–xxxiv, 1992\.
  * [31] S. Jin, J. Pulido, P. Grosset, J. Tian, D. Tao, and J. Ahrens, “Adaptive configuration of in situ lossy compression for cosmology simulations via fine-grained rate-quality modeling,” in _Proceedings of the 30th International Symposium on High-Performance Parallel and Distributed Computing_ , 2020, pp. 45–56.
  * [32] SZ3: A Modular Error-bounded Lossy Compression Framework for Scientific Datasets, https://github.com/szcompressor/SZ3.
  * [33] X. Liang, S. Di, D. Tao, S. Li, S. Li, H. Guo, Z. Chen, and F. Cappello, “Error-controlled lossy compression optimized for high compression ratios of scientific datasets,” in _2018 IEEE International Conference on Big Data_. IEEE, 2018, pp. 438–447.
  * [34] D. Tao, S. Di, Z. Chen, and F. Cappello, “Significantly improving lossy compression for scientific data sets based on multidimensional prediction and error-controlled quantization,” in _2017 IEEE International Parallel and Distributed Processing Symposium_. IEEE, 2017, pp. 1129–1139.
  * [35] S. Di and F. Cappello, “Fast error-bounded lossy hpc data compression with SZ,” in _2016 IEEE International Parallel and Distributed Processing Symposium_. Chicago, IL, USA: IEEE, 2016, pp. 730–739.
  * [36] Q. Zhou, C. Chu, N. Kumar, P. Kousha, S. Ghazimirsaeed, H. Subramoni, and D. Panda, “Designing high-performance mpi libraries with on-the-fly compression for modern gpu clusters,” in _2021 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 2021, pp. 444–453.
  * [37] K. J. Bowers, B. Albright, L. Yin, B. Bergen, and T. Kwan, “Ultrahigh performance three-dimensional electromagnetic relativistic kinetic plasma simulation,” _Physics of Plasmas_ , vol. 15, no. 5, p. 055703, 2008.
  * [38] S. Byna, J. Chou, O. Rubel, H. Karimabadi, W. S. Daughter, V. Roytershteyn, E. W. Bethel, M. Howison, K.-J. Hsu, K.-W. Lin _et al._ , “Parallel i/o, analysis, and visualization of a trillion particle simulation,” in _SC’12: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis_. IEEE, 2012, pp. 1–12.
  * [39] B. Behzad, H. V. T. Luu, J. Huchette, S. Byna, R. Aydt, Q. Koziol, M. Snir _et al._ , “Taming parallel i/o complexity with auto-tuning,” in _SC’13: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis_. IEEE, 2013, pp. 1–12.
  * [40] B. Behzad, S. Byna, S. M. Wild, M. Prabhat, and M. Snir, “Improving parallel i/o autotuning with performance modeling,” in _Proceedings of the 23rd international symposium on High-performance parallel and distributed computing_ , 2014, pp. 253–256.
  * [41] W. Bhimji, D. Bard, M. Romanus, D. Paul, A. Ovsyannikov, B. Friesen, M. Bryson, J. Correa, G. K. Lockwood, V. Tsulaia _et al._ , “Accelerating science with the nersc burst buffer early user program,” 2016.
  * [42] J. Han, D. Koo, G. K. Lockwood, J. Lee, H. Eom, and S. Hwang, “Accelerating a burst buffer via user-level i/o isolation,” in _2017 IEEE International Conference on Cluster Computing (CLUSTER)_. IEEE, 2017, pp. 245–255.
  * [43] D. Koo, J. Lee, J. Liu, E.-K. Byun, J.-H. Kwak, G. K. Lockwood, S. Hwang, K. Antypas, K. Wu, and H. Eom, “An empirical study of i/o separation for burst buffers in hpc systems,” _Journal of Parallel and Distributed Computing_ , vol. 148, pp. 96–108, 2021.
  * [44] S. Habib, A. Pope, H. Finkel, N. Frontiere, K. Heitmann, D. Daniel, P. Fasel, V. Morozov, G. Zagaris, T. Peterka, V. Venkatram, L. Zarija, S. Saba, and W.-k. Liao, “Hacc: Simulating sky surveys on state-of-the-art supercomputing architectures,” _New Astronomy_ , vol. 42, pp. 49–65, 2016.
  * [45] Summit supercomputer, https://www.olcf.ornl.gov/summit/, 2020.
  * [46] W. Daughton, J. Scudder, and H. Karimabadi, “Fully kinetic simulations of undriven magnetic reconnection with open boundary conditions,” _Physics of Plasmas_ , vol. 13, no. 7, p. 072101, 2006.
  * [47] A. K. Paul, O. Faaland, A. Moody, E. Gonsiorowski, K. Mohror, and A. R. Butt, “Understanding hpc application i/o behavior using system level statistics,” in _2020 IEEE 27th International Conference on High Performance Computing, Data, and Analytics (HiPC)_. IEEE, 2020, pp. 202–211.
  * [48] S. Lee, K.-y. Hou, K. Wang, S. Sehrish, M. Paterno, J. Kowalkowski, Q. Koziol, R. B. Ross, A. Agrawal, A. Choudhary _et al._ , “A case study on parallel hdf5 dataset concatenation for high energy physics data analysis,” _Parallel Computing_ , vol. 110, p. 102877, 2022.
  * [49] S. Lang, P. Carns, R. Latham, R. Ross, K. Harms, and W. Allcock, “I/o performance challenges at leadership scale,” in _Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis_. IEEE, 2009, pp. 1–12.
  * [50] B. Xie, H. Tang, S. Byna, Q. Koziol, and O. Sarp, “Tuning I/O Performance on Summit – HDF5 Write Use Case Study,” https://hps.vi4io.org/_media/events/2020/hpciodc20-hdf5.pdf, 2020, [Online; accessed 01-April-2022].
  * [51] S. Jin, G. Li, S. L. Song, and D. Tao, “A novel memory-efficient deep learning training framework via error-bounded lossy compression,” in _Proceedings of the 26th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming_ , 2021, pp. 485–487.
  * [52] S. Byna, A. Uselton, D. K. Prabhat, and Y. He, “Trillion particles, 120,000 cores, and 350 tbs: Lessons learned from a hero i/o run on hopper,” in _Cray user group meeting_ , 2013.
  * [53] LCF, “Bebop - laboratory computing resource center,” https://www.lcrc.anl.gov/systems/resources/bebop/, (Accessed on 11/17/2021).