#  Sketching and Streaming for Dictionary Compression

Ruben Becker∗ Matteo Canton† Davide Cenzato∗
Sung-Hwan Kim∗ Bojana Kodric∗ and Nicola Prezza∗
∗Ca’ Foscari University of Venice †University of Udine Via Torino 155 Via
delle Scienze 206 30172 Venezia, Italy 33100 Udine, Italy
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We initiate the study of sub-linear sketching and streaming techniques for
estimating the output size of common dictionary compressors such as Lempel-Ziv
’77, the run-length Burrows-Wheeler transform, and grammar compression. To
this end, we focus on a measure that has recently gained much attention in the
information-theoretic community and which approximates up to a polylogarithmic
multiplicative factor the output sizes of those compressors: the normalized
substring complexity function $\delta$. As a matter of fact, $\delta$ itself
is a very accurate measure of compressibility: it is monotone under
concatenation, invariant under reversals and alphabet permutations, sub-
additive, and asymptotically tight (in terms of worst-case entropy) for
representing strings, up to polylogarithmic factors.

We present a data sketch of $O(\varepsilon^{-3}\log
n+\varepsilon^{-1}\log^{2}n)$ words that allows computing a multiplicative
$(1\pm\varepsilon)$-approximation of $\delta$ with high probability, where $n$
is the string length. The sketches of two strings $S_{1},S_{2}$ can be merged
in $O(\varepsilon^{-1}\log^{2}n)$ time to yield the sketch of
$\\{S_{1},S_{2}\\}$, speeding up the computation of _Normalized Compression
Distances_ (NCD). If random access is available on the input, our sketch can
be updated in $O(\varepsilon^{-1}\log^{2}n)$ time for each character right-
extension of the string. This yields a polylogarithmic-space algorithm for
approximating $\delta$, improving exponentially over the working space of the
state-of-the-art algorithms running in nearly-linear time. Motivated by the
fact that random access is not always available on the input data, we then
present a streaming algorithm computing our sketch in $O(\sqrt{n}\cdot\log n)$
working space and $O(\varepsilon^{-1}\log^{2}n)$ worst-case delay per
character. We show that an implementation of our streaming algorithm can
estimate $\delta$ on a dataset of 189GB with a throughput of 203MB per minute
while using only 5MB of RAM, and that our sketch speeds up the computation of
all-pairs NCD distances by one order of magnitude, with applications to
phylogenetic tree reconstruction.

##  1 Introduction

Sketching techniques allow to summarize in sub-linear space information on big
datasets, enabling the approximation of useful statistics such as high-order
moments [1], norms [2], and frequencies [3] (to name a few). Additionally,
most data sketches can be computed on data streams in sub-linear space, making
them attractive in big data scenarios. In this paper, we consider data
sketches summarizing the _information content_ of a string as approximated by
data compression techniques. Previous research on this problem has focused on
_empirical entropy_. Chakrabarti et al. [4] showed that the zero-order
empirical entropy $H_{0}$ of a data stream can be efficiently approximated up
to a multiplicative $(1+\varepsilon)$-factor in poly-logarithmic space, but
any multiplicative approximation of the $k$-th order entropy $H_{k}$ requires
nearly-linear space for $k\geq 1$. In addition to this fact, it is well-known
that $H_{k}$ is a weak measure when the dataset is highly repetitive [5]. As
extensively shown in the literature (see, for example, the survey by Navarro
[6]), dictionary compression measures such as the number $z$ of phrases of the
Lempel-Ziv’77 factorization (used by winzip, 7-zip, gzip, xz), the number $r$
of equal-letter runs in the Burrows-Wheeler transform (used by bzip-2), and
the size $g$ of a smallest context-free grammar generating (only) the text,
are exempt from such a limitation. The information-theoretic quality of these
measures is strengthened by the fact that _Normalized Compression Distances_
based on dictionary compressors yield very precise notions of string
similarity [7]. Sketching and streaming techniques for such measures would
thus speed up tasks such as the computation of all-pairs similarities when the
underlying metric is based on data compression (useful, for example, in the
computation of phylogenetic trees [7]).

Motivated by the above considerations, in this paper we present the first sub-
linear-space sketching and streaming techniques for estimating the output
sizes of dictionary compressors. This result is obtained by describing a data
sketch yielding a $(1\pm\varepsilon)$-approximation of the _normalized
substring complexity_ $\delta=\max_{k\geq 1}\\{d_{k}/k\\}$, where $d_{k}$ is
the number of distinct length-$k$ substrings of the string, a measure
introduced by Raskhodnikova et al. in [8]. As shown by Kociumaka et al. [9]
and Kempa and Kociumaka [10], any of the above dictionary compression measures
is lower-bounded by $\delta$ and upper-bounded by $\delta(\log n)^{c}$, where
$n$ is the string’s length and $c$ is an opportune constant depending on the
compressor. Even better, Bonnie et al. in [11] experimentally showed that
$\delta$, $z$, and $r$ (normalized to the interval $[0,1]$) are almost
indistinguishable on collections of genomic data. As a matter of fact,
$\delta$ is known to be an even more accurate information measure than $z$,
$r$, and $g$: it is monotone under string concatenation, invariant under
reversals and alphabet permutations, sub-additive, and asymptotically tight
(in terms of worst-case entropy) for representing strings, up to
polylogarithmic factors [9]. None of the measures $z,r,g$ possesses
_simultaneously_ all of these properties.

#### Overview of the paper.

After providing all necessary definitions in Section 2, in Section 3 we prove
new properties of the normalized substring complexity $\delta$ and of the
Normalized Compression Distance [7] $\operatorname{NCD}_{\delta}$ based on
$\delta$. In particular, we show that $\delta$ is perfectly sub-additive, that
$\operatorname{NCD}_{\delta}(x,y)$ always lies in $[0,1]$ (according to [12],
this is an indicator that $\delta$ is a compressibility measure of good
quality), and that $\operatorname{NCD}_{\tilde{\delta}}(x,y)$ is an additive
$\Theta(\varepsilon)$-approximation of $\operatorname{NCD}_{\delta}(x,y)$ if
$\tilde{\delta}$ is a multiplicative $(1\pm\varepsilon)$-approximation of
$\delta$. This motivates designing data sketches for $\delta$, a problem that
we solve in Section 4. Our sketch is based on the observation (already noted
in [13] for the particular case $\varepsilon=1$) that $\max_{i\geq
0}\\{d_{\lceil(1+\varepsilon)^{i}\rceil}/\lceil(1+\varepsilon)^{i}\rceil\\}$
is a $(1-\Theta(\varepsilon))$-approximation of $\delta$. We approximate
$d_{k}$, for each sampled length $k=\lceil(1+\varepsilon)^{i}\rceil$, by
keeping a count-distinct sketch [14] for the subset of distinct (Rabin’s
fingerprints [15] of the) length-$k$ substrings. Our sketch uses space
polynomial in $\varepsilon^{-1}\log n$ and supports updates and queries
(returning a $(1\pm\varepsilon)$ approximation of $\delta$), in
$O(\varepsilon^{-1}\log^{2}n)$ time. The sketches of two strings $S_{1}$ and
$S_{2}$ can moreover be merged in $O(\varepsilon^{-1}\log^{2}n)$ time to
obtain the sketch of $\\{S_{1},S_{2}\\}$, from which one can compute an
additive $\varepsilon$-approximation of
$\operatorname{NCD}_{\delta}(S_{1},S_{2})$. In Section 5 we show how to
compute our sketch in sub-linear space on an input stream of length $n$. The
main difficulty in achieving sub-linear space is that, in order to compute the
Rabin’s fingerprints of the stream’s length-$k$ substrings, we need random
access to the $k$-th most recent stream’s character. Since the largest $k$ for
which we need to compute $d_{k}$ is linear in $n$, storing the most recent $k$
characters would require $\Theta(n)$ working space. Our solution relies on the
observation that, if $\hat{k}=\operatorname{argmax}_{k\geq 1}\\{d_{k}/k\\}$ is
small, then we can afford keeping a sliding window of the last $\hat{k}$
stream’s characters. If, on the other hand, $\hat{k}$ is large, then the
stream is highly repetitive so we can compress it in small space while
supporting bookmarked access to its characters. We conclude in Section 6 with
experimental results. Complete proofs can be found in the full version [16].

#### Related work.

Bonnie et al. [11] have already observed that $d_{k}$ can be efficiently
estimated by employing count-distinct sketches, and that this can yield an
heuristic algorithm for estimating $\delta$. Their strategy relies on
estimating $d_{k}/k$ for increasing values of $k$, until a local maximum is
found. While this strategy works well in practice because, as they showed,
$\hat{k}=\operatorname{argmax}_{k\geq 1}\\{d_{k}/k\\}$ tends to be a very
small number, on particular strings (for example, Thue-Morse) $\hat{k}$ is of
the order of $\Theta(n)$ and, as a result, computing all the sketches for
$d_{k}$ requires linear space and quadratic processing time in the worst case.
Moreover, local maxima of $d_{k}/k$ do not always coincide with the global
maximum, so this strategy does not yield any provable approximation of
$\delta$. We are not aware of other works in the literature describing data
sketches for estimating the output sizes of _dictionary_ compressors (the
literature on estimating empirical entropy is, on the other end, much richer:
see [4] and references therein). Our results can be viewed also as a space-
efficient way to approximate measure $\delta$. Christiansen et al. [17] showed
how to compute $\delta$ for a given string $T$ in linear time and space.
Recently, Bernardini et al. [13] provided space-time trade-offs for
computing/approximating $\delta$ in sub-linear working space on top of the
input string. If $O(n\operatorname{polylog}n)$ time is allowed, their
algorithms require $\Theta(n/\operatorname{polylog}n)$ working space, which
they proved to be optimal for computing $\delta$ exactly. Our algorithm, on
the other hand, computes a multiplicative $(1\pm\varepsilon)$-approximation of
$\delta$ using working space polynomial in $\varepsilon^{-1}\log n$.

##  2 Preliminaries

We denote $[n]:=\\{1,\ldots,n\\}$ for any integer $n$ ($[n]=\emptyset$ for
$n\leq 0$). For $a\in\mathbb{R}^{+}$ and a real number $\varepsilon\in[0,1]$,
we write $[(1\pm\varepsilon)a]$ for the interval
$[(1-\varepsilon)a,(1+\varepsilon)a]$. Similarly, we write $[a\pm\varepsilon]$
for the interval $[a-\varepsilon,a+\varepsilon]$.

We assume to be given a string $S$ of length $n>1$ over an alphabet $\Sigma$
of cardinality $\sigma>1$. For $k\geq 1$, we define
$D_{k}(S):=\\{S[i..i+k-1]:i\in[n-k+1]\\}$, i.e., the set of all distinct
substrings of length $k$ of $S$. Notice that $D_{k}(S)=\emptyset$ if $k>n$.
The _$k$ -substring complexity_ $d_{k}(S)$ of $S$ is the cardinality of this
set, i.e., $d_{k}(S):=|D_{k}(S)|$. The _normalized substring complexity_
$\delta$ is defined as follows:

$\delta(S):=\max_{k\geq 1}\\{|D_{k}(S)|/k\\}=\max_{k\geq 1}\\{d_{k}(S)/k\\}.$

We omit the argument from $D_{k}$, $d_{k}$, and $\delta$ in case it is clear
from the context. Here, we also extend this measure to pairs of strings $S$
and $T$. Rather than using $\delta(ST)$, we propose the following natural
definition that does not take into account artificial length-$k$ substrings
crossing the border between $S$ and $T$:

$\delta(S,T):=\max_{k\geq 1}\\{|D_{k}(S)\cup D_{k}(T)|/k\\}.$

This version also gives mathematically cleaner results (e.g., perfect sub-
additivity) and, in any case, differs from $\delta(ST)$ by at most 1. As a
consequence, most of our results (read also below) hold also by replacing
$\delta(S,T)$ with $\delta(ST)$.

The _Normalized Compression Distance_ has been defined by Cilibrasi and
Vitányi [7] as a proxy for the non-computable Normalized Information Distance
[18]. For two strings $S$ and $T$ and an arbitrary compressibility measure $Z$
(for example, the output size of compression software such as gzip and xz), it
is defined as

$\operatorname{NCD}_{Z}(S,T):=\frac{Z(S,T)-\min\\{Z(S),Z(T)\\}}{\max\\{Z(S),Z(T)\\}}.$

Given a uniform prime $q=n^{\Theta(1)}$, the _Rabin’s fingerprint_ [15] of $S$
is defined as $\rho(S)=\sum_{i=1}^{n}S[i]\cdot\sigma^{n-i}\mod q$. Collisions
between substrings of $S$ through $\rho$ happen with low probability, so the
results of our paper hold with high probability. We extensively use the fact
that the fingerprint of the concatenation of two strings $S_{1},S_{2}$ can be
computed in constant time from (i) the fingerprints of $S_{1}$ and $S_{2}$ and
(ii) $\sigma^{|S_{2}|}\mod q$ (see [15]). Given a set $B\subseteq\Sigma^{*}$
of strings, we define $\rho(B)=\\{\rho(s)\ :\ s\in B\\}$.

Given a set $U$, a _count-distinct sketch_ $CD(U)$ is a sub-linear-space data
structure supporting three main operations: $CD(U).add(x)$, which turns the
sketch into $CD(U\cup\\{x\\})$, $CD(U_{1}).merge(CD(U_{2}))$, which turns the
sketch into $CD(U_{1}\cup U_{2})$, and $CD(U).estimate()$, which returns a
$(1\pm\varepsilon)$ approximation of $|U|$. In our work, we use the optimal
count-distinct sketch of Kane et al. [14]. Letting $U\subseteq[u]$, this
sketch uses $O(\varepsilon^{-2}+\log u)$ words of space and computes a
$(1\pm\varepsilon)$ approximation of $|U|$ with high probability of success.
All operations are supported in $O(\log u)$ time111The authors claim
$O(\varepsilon^{-2}+\log u)$ _bits_ of space and $2/3$ success probability,
which can be amplified by taking the median of $\Theta(\log u)$ sketches (thus
yielding the bounds we claim above). In our paper, the universe is composed by
Rabin’s fingerprints and has therefore size $u=n^{\Theta(1)}$..

Assume $S[1]=\$$, where $\$$ is lexicographically smaller than all other
alphabet’s characters and does not appear anywhere else in $S$. The _Burrows-
Wheeler transform_ (BWT) of the reverse $S^{R}$ of $S$ is obtained by sorting
lexicographically all suffixes of $S^{R}$ and then taking, in this order, the
character preceding each suffix. For example, if $S=\$babba$, then the sorted
suffixes and the BWT of $S^{R}$ are shown in Table 1.

suffixes of $S^{R}$ | BWT
---|---
$ | b
a b $ | b
a b b a b $ | _$_
b $ | a
b a b $ | _b_
b b a b $ | a
Table 1: Burrows-Wheeler transform $BWT(S^{R})$ of the string $S^{R}=abbab\$$.

The _LF property_ of the BWT states that the $i$-th occurrence of $c\in\Sigma$
in the BWT corresponds to the position of the $i$-th suffix starting with
$c\in\Sigma$ in Table 1. The _LF function_ is the permutation of $[1,n]$
implementing this observation: for instance, in the above example
$BWT.LF(2)=5$ because character $BWT[2]$ corresponds to the first character
($b$) of the fifth (in lexicographic order) suffix $bab\$$. We denote with $r$
the number of _equal-letter runs_ of the BWT; in the above example, $r=5$
(runs are highlighted in alternating bold/italic). We moreover use the
following result:

######  Lemma 2.1 ([19], Thm. 2).

Letting $S$ be a string and $r$ be the number of equal-letter runs in
$BWT(S^{R})$, there exists a data structure of $O(r)$ words storing
$BWT(S^{R})$ supporting right-extensions of $S$ (i.e. $BWT(S^{R})\rightarrow
BWT((Sa)^{R})$, for any $a\in\Sigma$) in $O(\log|S|)$ time. Within the same
time, the structure supports computing the LF function and retrieving any
character of $BWT(S^{R})$.

##  3 Properties of $\delta$ and $\operatorname{NCD}_{\delta}$

We start by proving some properties of $\delta$. Proofs of some statements are
omitted due to space limitations and can be found in the full version [16].
The first main property that we show is that $\delta$ is both sub-additive and
monotone in the following sense.

######  Lemma 3.1.

For any strings $S$ and $T$,
$\max\\{\delta(S),\delta(T)\\}\leq\delta(S,T)\leq\delta(S)+\delta(T)$.

The proof of the lemma uses the properties of the corresponding maximizers
together with the fact that the union is a superset of both its arguments
(left inequality) and that the union is of smaller cardinality than the sum of
the cardinalities of its arguments (right inequality). We remark that it is a
well-known fact that the monotonicity property holds for the case of
concatenation of the two strings [9]. Using the sub-additivity of $\delta$, we
obtain:

######  Corollary 3.2.

For any strings $S$ and $T$ it holds that
$0\leq\operatorname{NCD}_{\delta}(S,T)\leq 1$.

To see why this holds, assume, w.l.o.g., that
$\max\\{\delta(S),\delta(T)\\}=\delta(S)$. Then,
$\operatorname{NCD}_{\delta}(S,T)=\frac{\delta(S,T)-\delta(T)}{\delta(S)}\geq\frac{\delta(T)-\delta(T)}{\delta(S)}=0$
and
$\operatorname{NCD}_{\delta}(S,T)\leq\frac{\delta(S)+\delta(T)-\delta(T)}{\delta(S)}=1$.
Ming et al. [12] state that common compressors yield a normalized compression
distance between $0$ and $1+\varepsilon$, where the $\varepsilon$ is due to
“imperfections” of the compression algorithm. Above we proved that in the case
of the normalized substring complexity $\delta$, the corresponding
$\varepsilon$ is equal to $0$.

We conclude by showing that a multiplicative approximation of $\delta$ can be
used to obtain an additive approximation of the Normalized Compression
Distance $\operatorname{NCD}_{\delta}$.

######  Lemma 3.3.

Let $\varepsilon\in(0,1)$, $\varepsilon^{\prime}:=\varepsilon/5$ and let $S$
and $T$ be two strings. Assume that $\tilde{\delta}(S)$, $\tilde{\delta}(T)$,
and $\tilde{\delta}(S,T)$ are approximations of $\delta$ in the sense that
$\tilde{\delta}(S)\in[(1\pm\varepsilon^{\prime})\delta(S)]$,
$\tilde{\delta}(T)\in[(1\pm\varepsilon^{\prime})\delta(T)]$, as well as
$\tilde{\delta}(S,T)\in[(1\pm\varepsilon^{\prime})\delta(S,T)]$. Then

$\operatorname{NCD}_{\tilde{\delta}}(S,T)\in[\operatorname{NCD}_{\delta}(S,T)\pm\varepsilon].$

We prove this lemma by using the facts that $\tilde{\delta}$ is a
multiplicative approximation of $\delta$, that $\delta$ is sub-additive (see
Lemma 3.1), and that $\operatorname{NCD}_{\delta}\in[0,1]$ (see Corollary
3.2).

##  4 A data sketch for estimating $\delta$

We introduce our data sketch, then prove that it yields a good approximation
of $\delta$.

######  Definition 4.1 (Sketch for $\delta$).

Let $S$ be a string,
$A:=\\{\lceil\alpha^{i}\rceil:i\in[\lfloor\log_{\alpha}n\rfloor]\\}$ be a set
of sampled lengths for some real number (sample rate) $\alpha>1$, and
$CD_{k}=CD(\rho(D_{k}(S)))$, where $CD$ is the count-distinct sketch described
in Section 2 and $\rho$ is Rabin’s hash function. Our data sketch is defined
as $\kappa(S)=\langle CD_{k}\ :\ k\in A\rangle.$

We define $\kappa(S).estimate()=\max\\{CD_{k}.estimate()/k\ :\ k\in A\\}$.
When extending the stream $S$ with a new character $a$, yielding string $Sa$,
the sketch is updated by calling $CD_{k}.add(\rho(S[|S|-k+2,|S|]a))$ for all
$k\in A$. We denote this operation by $\kappa(S).extend(a)$. Note that, if
constant-time random access is available on $S$ and if $\sigma^{k-1}\mod q$
has been pre-computed for all $k\in A$ (in $O(\varepsilon^{-1}\log^{2}n)$
time), $\rho(S[|S|-k+2,|S|]a)$ can be computed in constant time from
$\rho(S[|S|-k+1,|S|])$; see [15]. Finally,
$\kappa(S_{1}).merge(\kappa(S_{2}))$ returns the sketch
$\kappa(\\{S_{1},S_{2}\\})=\langle CD^{\prime}_{k}\ :\ k\in A\rangle$, where
$CD^{\prime}_{k}=CD^{1}_{k}.merge(CD^{2}_{k})$ and $CD^{i}_{k}$ is the count-
distinct sketch for the (fingerprints of the) length-$k$ substrings of
$S_{i}$, $i\in\\{1,2\\}$. Operation $extend(a)$ is not defined when the sketch
represents a set of strings; this is not an issue, since we will call $merge$
only to estimate $\delta(S_{1},S_{2})$ and $NCD_{\delta}(S_{1},S_{2})$.

With the next theorem we show that $\kappa(S).estimate()$ returns a
multiplicative $(1\pm\varepsilon)$-approximation of $\delta(S)$ (analogous for
$\kappa(\\{S_{1},S_{2}\\}).estimate()$).

######  Lemma 4.2.

Let $S$ be a string of length $n$. Let $\varepsilon>0$,
$\varepsilon^{\prime}=\varepsilon/4$, and $\alpha=1+\varepsilon^{\prime}$.
Assume that $\tilde{d}_{k}\in[(1\pm\varepsilon^{\prime})d_{k}(S)]$ for all
$k\in A:=\\{\lceil\alpha^{i}\rceil:i\in[\lfloor\log_{\alpha}n\rfloor]\\}$,
then $\tilde{\delta}:=\max\\{\tilde{d}_{k}/k:k\in A\\},$ satisfies
$\tilde{\delta}\in[(1\pm\varepsilon)\delta]$.

We show this theorem by quantifying the impact of two types of errors on the
quantity $\delta$. These two types are (1) the error obtained when
approximating the values $d_{k}$ by $\tilde{d}_{k}$ and (2) the error due to
the restriction of the string’s offsets $[n]$ to the set $A$. The error of
type (1) directly implies an error of the same magnitude
$(1\pm\varepsilon^{\prime})$ on $\delta$. We note that this error actually
itself has two sources, namely (1.1) errors due to collisions when computing
Rabin’s fingerprints and (1.2) errors due to the count-distinct sketch when
applied to the fingerprints. Both of these errors are accounted for in the
assumption $\tilde{d}_{k}\in[(1\pm\varepsilon^{\prime})d_{k}(S)]$. The error
of type (2) instead is more subtle to analyse – the main observation here is
that $d_{j+1}\geq d_{j}-1$ for every $j$, as every distinct length-$j$
substring other than possibly $S[n-j+1..n]$ gives at least one distinct
length-$(j+1)$ substring. Now assume that $i\in[n]\setminus A$ and that $a\in
A$ is the minimum element of $A$ larger than $i$. Then applying the previous
observation iteratively yields $d_{a}\geq d_{i}-\beta$, where $\beta=a-i$.
Hence, we can quantify how much $\delta$ gets “perturbed” by restricting to
the subset $A$ of the string’s offsets $[n]$.

From Lemmas 3.3 and 4.2, the sketch of Definition 4.1 yields a multiplicative
$(1\pm\varepsilon)$-approximation of $\delta$ and an additive
$\varepsilon$-approximation of $\operatorname{NCD}_{\delta}$ if $CD$ is the
count-distinct sketch of [14] with error rate $\varepsilon/20$, and the set
$A$ is built with sample rate $\alpha=1+\varepsilon/20$. From [14] and since
$|A|\in\Theta(\varepsilon^{-1}\log n)$, our data sketch uses
$\Theta(\varepsilon^{-3}\log n+\varepsilon^{-1}\log^{2}n)$ words of space and
supports all operations in time $O(\varepsilon^{-1}\log^{2}n)$.

Using repeatedly operation $extend$ on our data sketch we immediately obtain:

######  Theorem 4.3.

For any string $S$ of length $n$ supporting random access in time at most
$O(\log n),$ and any approximation rate $\varepsilon>0$, we can compute a
multiplicative $(1\pm\varepsilon)$-approximation of $\delta(S)$ in
$O(\varepsilon^{-1}n\log^{2}n)$ time using $\Theta(\varepsilon^{-3}\log
n+\varepsilon^{-1}\log^{2}n)$ words of working space on top of the string. The
result is correct with high probability.

##  5 Streaming algorithm

We now show how to compute the sketch of Definition 4.1 in $O(\sqrt{n}\log n)$
words of working space (on top of the sketch) with one pass over the streamed
input string.

Let $S$ denote the current stream, and $S^{R}$ be the reversed stream. We
assume that an upper-bound $n$ to the maximum stream length is known before
the algorithm starts. Let $r$ be the number of equal-letter runs in the
Burrows-Wheeler transform of $S^{R}$. By [10] and by the fact that $\delta$ is
invariant under string reversals, it holds $r\leq 8\delta\log^{2}n$. Our
streaming algorithm works as follows. We keep a sliding window
$S[|S|-K+1,|S|]$ of the last $K$ stream characters, for some parameter $K$ to
be determined later, and at the same time we keep a dynamic run-length BWT
(RLBWT) of $S^{R}$ that we update by appending the stream’s characters using
Lemma 2.1. Before the algorithm starts, in $O(\varepsilon^{-1}\log^{2}n)$ time
we compute $\sigma^{k-1}\mod q$ for all the $|A|\in O(\varepsilon^{-1}\log n)$
sampled substring lengths $k\in A$ in our sketch, using fast exponentiation.

Let $k\in A$ be one of the sampled string lengths in our sketch, and let $a$
be a new character arriving on the stream (so that the new stream is $Sa$). In
order to update our sketch, we need to compute the fingerprint of the last $k$
stream’s characters: $\rho(S[|S|-k+2,|S|]a)$. At any stage of the algorithm,
we keep the Rabin’s fingerprint $\rho(S[2,|S|])$ of the whole stream,
excluding character $S[1]=\$$. If $|Sa|=k+1$, then $\rho(S[|S|-k+2,|S|]a)$ is
equal to the Rabin’s fingerprint of the whole stream. Otherwise, if $|Sa|>k+1$
then $\rho(S[|S|-k+1,|S|])$ has already been computed in the previous steps
and we can use the formula
$\rho(S[|S|-k+2,|S|]a)=(\rho(S[|S|-k+1,|S|])-S[|S|-k+1]\cdot\sigma^{k-1})\cdot\sigma+a\mod
q$. As a result, updating the fingerprint reduces to extracting character
$S[|S|-k+1]$. We use the window $S[|S|-K+1,|S|]$ to extract $S[|S|-k+1]$ for
any $k\leq K$, and the RLBWT to extract $S[|S|-k+1]$ for any $k>K$ using a
bookmarking technique that we sketch in Figure 1 and we describe in full
detail in the full version [16]. This allows us to update the Rabin’s
fingerprints for all sampled substring lengths $k$ and thus to implement
operation $extend(a)$.

We now describe the policy we employ to keep space usage under control. Let
$r^{\prime}$ be the number of equal-letter runs in the BWT obtained by
ignoring (removing) character $\$$. It is easy to see that (i) $r-2\leq
r^{\prime}\leq r$ and (ii) $r^{\prime}$ is non-decreasing upon appending
characters at the end of the stream. As soon as $r^{\prime}\geq
8n(\log^{2}n)/K$, we discard the RLBWT and keep only the sliding window for
the rest of the stream. As a consequence, from this point on we are only able
to extract (fingerprints of) length-$k$ substrings with $k\leq K$. However, we
show that this is enough: if we discard the RLBWT, then it means that
$\delta\geq\frac{r}{8\log^{2}n}\geq\frac{r^{\prime}}{8\log^{2}n}\geq n/K$. Let
$\hat{k}=\operatorname{argmax}_{k\geq 1}\\{d_{k}/k\\}$. Then,
$\hat{k}=d_{\hat{k}}/\delta\leq n/\delta\leq K$ so to compute $\delta$ on the
rest of the stream we can focus only on the length-$k$ substrings with $k\leq
K$.

The sliding window $S[|S|-K+1,|S|]$ uses $K$ words of space. We discard the
RLBWT when $r^{\prime}\geq 8n(\log^{2}n)/K$, so (since $r\leq r^{\prime}+2$)
this structure always uses at most $O(r)\subseteq O(n(\log^{2}n)/K)$ words. As
a consequence, in total we use $O(K+n(\log^{2}n)/K)$ words of space, which is
optimized asymptotically when $K=\sqrt{n}\log n$; then, our algorithm uses at
most $O(\sqrt{n}\log n)$ words of space.

suffixes of $S^{R}$ | BWT
---|---
$ | $

suffixes of $S^{R}$ | BWT
---|---
$ | b
b$ | $

suffixes of $S^{R}$ | BWT |
---|---|---
$ | b | $j=1$
ab$ | $ |
b$ | a |

suffixes of $S^{R}$ | BWT |
---|---|---
$ | b |
aab$ | $ |
ab$ | a |
b$ | a | $j=4$

Figure 1: Example showing how the BWT (of the reversed stream) is updated upon
character right-extensions of the stream, and how the bookmark $j$
corresponding to window length $k=2$ is initialized and updated. Top left:
empty stream ($S=\$$). Top right: a new character $b$ arrives on the stream
($S=\$b$): in the BWT, $ is replaced by $b$ and a new $ is inserted in the
position corresponding to the lexicographic rank $i=2$ of the new suffix
$b\$$. Position $i$ is computed in $O(\log|S|)$ time using the algorithm
described in [19, Thm. 2]. Bottom left: a new character $a$ arrives on the
stream ($S=\$ba$): in the BWT, $ is replaced by $a$ and a new $ is inserted in
the position corresponding to the lexicographic rank $i=2$ of the new suffix
$ab\$$. Since the stream length is equal to $k+1=3$, we initialize the
bookmark $j\leftarrow BWT.LF(i)=BWT.LF(2)=1$. Note that $BWT[j]=b$ indeed
contains character $S[|S|-k+1]=b$. Bottom right: a new character $a$ arrives
on the stream ($S=\$baa$): in the BWT, $ is replaced by $a$ and a new $ is
inserted in the position corresponding to the lexicographic rank $i=2$ of the
new suffix $aab\$$. Since $1=j<i=2$ ($ is inserted after position $j$), $j=1$
is not modified (otherwise, it would have been incremented by 1). Finally, we
update $j$ by advancing it by one position in the text: $j\leftarrow
BWT.LF(j)=BWT.LF(1)=4$. Note that $BWT[j]=a$ indeed contains character
$S[|S|-k+1]=a$. Importantly, the data structure of [19, Thm. 2] uses always a
space proportional to the number $r$ of equal-letter runs of the BWT.

We keep one bookmark (a position in the BWT) for every sampled length $k\in
A$, so our bookmarking technique does not affect the asymptotic working space
if $\varepsilon\geq n^{-1/2}$ (i.e. $|A|\leq\sqrt{n}\log n$). Updating each
bookmark and extracting $S[|S|-k+1]$ from the RLBWT take $O(\log n)$ time by
Lemma 2.1. This running time is absorbed by operation $merge()$ on the count-
distinct sketches, see Section 4. We obtain:

######  Theorem 5.1.

Given an upper-bound $n$ to the stream’s length, we can compute the sketch of
Definition 4.1 in $O(\sqrt{n}\log n)$ words of working space and
$O(\varepsilon^{-1}\log^{2}n)$ worst-case delay per stream character, for any
approximation factor $\varepsilon\geq n^{-1/2}$.

##  6 Implementation and experiments

We implemented a parallel version of our streaming algorithm in C++.222
https://github.com/regindex/substring-complexity We ran experiments on a
server with Intel(R) Xeon(R) W-2245 CPU @ 3.90GHz with 16 threads and 128GB of
RAM running Ubuntu 18.04 LTS 64-bit. Our complete experimental results are
reported in [16]. We used the repetitive real Pizza&Chilli dataset
(P&C)333https://pizzachili.dcc.uchile.cl/repcorpus/real/, large Canterbury
corpus444http://corpus.canterbury.ac.nz/resources/large.tar.gz, and datasets
from AF Project555https://afproject.org. We computed the relative error of our
approximation $\tilde{\delta}$ with respect to $\delta$ for different sampling
densities (i.e. parameter $\alpha$ of Definition 4.1). With the sparsest (less
precise) sampling scheme (option -p 1), $\tilde{\delta}$ always differed from
$\delta$ by up to 5% and the average throughput was of 174 MB per minute using
up to 16 threads (option -t 0). For efficiency reasons, the RLBWT is disabled
by default: in practice this does not affect precision, since
$\hat{k}=\operatorname{argmax}_{k}d_{k}/k$ was always extremely small
($\hat{k}\leq 100$ in all datasets), meaning that the RLBWT is never required.
We also computed ${\tilde{\delta}}$ on a big dataset of 189GB long reads of
Rana
Muscosa666https://trace.ncbi.nlm.nih.gov/Traces/?view=run_browser&acc=SRR11606868.
Our software finished the computation in 15:31 hours with a throughput of
203MB per minute using only about 5MB of internal memory.

Experiments on repetitiveness measures. We studied the effectiveness of
${\tilde{\delta}}$ as a repetitive measure. We compared it to exact $\delta$,
to the number of runs of the BWT $r$, to the number of phrases of the LZ77
parse $z$, and to the output of two popular compressors, xz and 7z. For each
dataset in the repetitive P&C corpus, we computed these five measures for
prefixes of increasing length. We observe that ${\tilde{\delta}}$ not only
follows closely the values of $\delta$, but it also mirrors the trend of the
other four measures. This suggests experimentally that ${\tilde{\delta}}$
computed by our streaming algorithm is a good indicator of repetitiveness and
compressibility.

Figure 2:  From left to right: $\tilde{\delta}$ error distribution on P&C
repetitive corpus, lineplot showing five repetitiveness measures (normalized
to $[0,1]$) computed on increasing prefixes of para.

Experiments on phylogeny reconstruction. We verified that NCD based on the
compression software $\mathtt{xz}$, on $\delta$, and on $\tilde{\delta}$ yield
similar phylogenetic trees with the Gene-trees dataset from AF Project; the
average normalized Robinson-Foulds distances ranged from 0.1 to 0.3,
indicating that the reconstructed trees were very similar. We also measured
the running time to compute all-pair NCDs on 29 sequences of average length
$\sim$81k. This process took only 3 minutes for $\tilde{\delta}$ and 24
minutes for the exact $\delta$, while for $\mathtt{xz}$ it required 42
minutes.

##  7 References

## References

  * [1] Noga Alon, Yossi Matias, and Mario Szegedy,  “The space complexity of approximating the frequency moments,”  in Proc. of 28th ACM STOC, 1996, pp. 20–29.
  * [2] William B. Johnson and Joram Lindenstrauss,  “Extensions of Lipschitz mappings into Hilbert space,”  Contemporary mathematics, vol. 26, pp. 189–206, 1984.
  * [3] Jayadev Misra and David Gries,  “Finding Repeated Elements,”  Sci. Comput. Program., vol. 2, no. 2, pp. 143–152, 1982.
  * [4] Amit Chakrabarti, Graham Cormode, and Andrew Mcgregor,  “A Near-Optimal Algorithm for Estimating the Entropy of a Stream,”  ACM Trans. Algorithms, vol. 6, no. 3, pp. 51:1–21, 2010.
  * [5] Sebastian Kreft and Gonzalo Navarro,  “On compressing and indexing repetitive sequences,”  Theoretical Computer Science, vol. 483, pp. 115–133, 2013.
  * [6] Gonzalo Navarro,  “Indexing Highly Repetitive String Collections, Part I: Repetitiveness Measures,”  ACM Comput. Surv., vol. 54, no. 2, pp. 29:1–29:31, 2022.
  * [7] Rudi Cilibrasi and Paul MB Vitányi,  “Clustering by compression,”  IEEE Transactions on Information theory, vol. 51, no. 4, pp. 1523–1545, 2005.
  * [8] Sofya Raskhodnikova, Dana Ron, Ronitt Rubinfeld, and Adam D. Smith,  “Sublinear Algorithms for Approximating String Compressibility,”  Algorithmica, vol. 65, no. 3, pp. 685–709, 2013.
  * [9] Tomasz Kociumaka, Gonzalo Navarro, and Nicola Prezza,  “Toward a Definitive Compressibility Measure for Repetitive Sequences,”  IEEE Trans. Inf. Theory, vol. 69, no. 4, pp. 2074–2092, 2023\.
  * [10] Dominik Kempa and Tomasz Kociumaka,  “Resolution of the Burrows-Wheeler Transform Conjecture,”  in Proc of 61st IEEE FOCS, 2020, pp. 1002–1013.
  * [11] Jessica K Bonnie, Omar Ahmed, and Ben Langmead,  “DandD: efficient measurement of sequence growth and similarity,”  bioRxiv, pp. 2023–02, 2023.
  * [12] Ming Li, Xin Chen, Xin Li, Bin Ma, and P.M.B. Vitanyi,  “The similarity metric,”  IEEE Transactions on Information Theory, vol. 50, no. 12, pp. 3250–3264, 2004.
  * [13] Giulia Bernardini, Gabriele Fici, Pawel Gawrychowski, and Solon P. Pissis,  “Substring Complexity in Sublinear Space,”  CoRR, vol. abs/2007.08357, 2020.
  * [14] Daniel M. Kane, Jelani Nelson, and David P. Woodruff,  “An Optimal Algorithm for the Distinct Elements Problem,”  in Proc. of 29th ACM PODS, 2010, p. 41–52.
  * [15] Michael O Rabin,  “Fingerprinting by random polynomials,”  Technical report, 1981.
  * [16] Ruben Becker, Matteo Canton, Davide Cenzato, Sung-Hwan Kim, Bojana Kodric, and Nicola Prezza,  “Sketching and Streaming for Dictionary Compression,” 2023,  https://arxiv.org/abs/2310.17980.
  * [17] Anders Roy Christiansen, Mikko Berggren Ettienne, Tomasz Kociumaka, Gonzalo Navarro, and Nicola Prezza,  “Optimal-Time Dictionary-Compressed Indexes,”  ACM Trans. Algorithms, vol. 17, no. 1, pp. 8:1–8:39, 2021.
  * [18] Sebastiaan A. Terwijn, Leen Torenvliet, and Paul M. B. Vitányi,  “Nonapproximability of the normalized information distance,”  J. Comput. Syst. Sci., vol. 77, no. 4, pp. 738–742, 2011.
  * [19] Alberto Policriti and Nicola Prezza,  “LZ77 computation based on the run-length encoded BWT,”  Algorithmica, vol. 80, no. 7, pp. 1986–2011, 2018.

##  Appendix A Missing Proofs

We start with the following two easy observations that we use in our proofs
later on.

######  Observation A.1.

If $S\neq a^{n}$ for $a\in\Sigma$, then $\delta(S)\geq d_{1}/1\geq 2/1=2$.

Note that if at least two distinct letters appear, it follows that $d_{1}>1$
and thus $\delta\geq d_{1}/1\geq 2$. Since it is easy to recognize the case
$S=a^{n}$ for some $a\in\Sigma$ in constant space and constant delay per
character, from now on we assume w.l.o.g. that $\delta\geq 2$.

We continue with the following simple observation that we will use in the
proof of Lemma 4.2.

######  Observation A.2.

It holds that $\hat{k}=\operatorname{argmax}_{k\geq 1}\\{d_{k}(S)/k\\}\leq
n/2$.

###### Proof.

Assume that $\hat{k}>n/2$. Then it is immediate that $d_{\hat{k}}\leq
n-\hat{k}+1<n/2$ (as the right most character in the substring can be at index
at most $n$). We now obtain that $d_{\hat{k}}/\hat{k}<1$, contradicting
Observation A.1. ∎

See 3.1

###### Proof.

Let $k_{S,T}$, $k_{S}$, and $k_{T}$ be such that
$\delta(S,T)=|D_{k_{S,T}}(S)\cup D_{k_{S,T}}(T)|/k_{S,T}$,
$\delta(S)=|D_{k_{S}}(S)|/k_{S}$, and $\delta(T)=|D_{k_{T}}(T)|/k_{T}$. Let,
w.l.o.g., $\delta(S)=\max\\{\delta(S),\delta(T)\\}$. Then,

$\displaystyle\delta(S,T)$ $\displaystyle=\frac{|D_{k_{S,T}}(S)\cup
D_{k_{S,T}}(T)|}{k_{S,T}}\geq\frac{|D_{k_{S}}(S)\cup
D_{k_{S}}(T)|}{k_{S}}\geq\frac{|D_{k_{S}}(S)|}{k_{S}}=\delta(S),$

where the second inequality uses the fact that $k_{S}$ is the maximizer for
$S$. For the second claim,

$\displaystyle\delta(S,T)$ $\displaystyle=\frac{|D_{k_{S,T}}(S)\cup
D_{k_{S,T}}(T)|}{k_{S,T}}\leq\frac{|D_{k_{S,T}}(S)|}{k_{S,T}}+\frac{|D_{k_{S,T}}(T)|}{k_{S,T}}$
$\displaystyle\leq\frac{|D_{k_{S}}(S)|}{k_{S}}+\frac{|D_{k_{T}}(T)|}{k_{T}}=\delta(S)+\delta(T),$

where the second inequality uses the fact that $k_{S}$ and $k_{T}$ are the
respective maximizers for $S$ and $T$. ∎

See 3.3

###### Proof.

We start with the lower bound. Using the definition of
$\operatorname{NCD}_{\tilde{\delta}}(S,T)$, we obtain

$\displaystyle\operatorname{NCD}_{\tilde{\delta}}(S,T)$
$\displaystyle\geq\frac{(1-\varepsilon^{\prime})\cdot\delta(S,T)-(1+\varepsilon^{\prime})\min\\{\delta(S),\delta(T)\\}}{(1+\varepsilon^{\prime})\max\\{\delta(S),\delta(T)\\}}$
$\displaystyle=\frac{1}{1+\varepsilon^{\prime}}\cdot\operatorname{NCD}_{\delta}(S,T)-\frac{\varepsilon^{\prime}}{1+\varepsilon^{\prime}}\cdot\frac{\delta(S,T)+\min\\{\delta(S),\delta(T)\\}}{\max\\{\delta(S),\delta(T)\\}}$
$\displaystyle=\operatorname{NCD}_{\delta}(S,T)-\frac{\varepsilon^{\prime}}{1+\varepsilon^{\prime}}\cdot\Big{(}\frac{\delta(S,T)+\min\\{\delta(S),\delta(T)\\}}{\max\\{\delta(S),\delta(T)\\}}+\operatorname{NCD}_{\delta}(S,T)\Big{)}$
$\displaystyle\geq\operatorname{NCD}_{\delta}(S,T)-\frac{4\varepsilon^{\prime}}{1+\varepsilon^{\prime}}$
$\displaystyle\geq\operatorname{NCD}_{\delta}(S,T)-\varepsilon$

using the sub-additivity of $\delta$ from Lemma 3.1, the fact that
$\operatorname{NCD}_{\delta}(S,T)\leq 1$ according to Corollary 3.2, and the
definition of $\varepsilon^{\prime}$. Similarly, now for the upper bound, we
obtain

$\displaystyle\operatorname{NCD}_{\tilde{\delta}}(S,T)$
$\displaystyle\leq\frac{(1+\varepsilon^{\prime})\cdot\delta(S,T)-(1-\varepsilon^{\prime})\min\\{\delta(S),\delta(T)\\}}{(1-\varepsilon^{\prime})\max\\{\delta(S),\delta(T)\\}}$
$\displaystyle=\frac{1}{1-\varepsilon^{\prime}}\cdot\operatorname{NCD}_{\delta}(S,T)+\frac{\varepsilon^{\prime}}{1-\varepsilon^{\prime}}\cdot\frac{\delta(S,T)+\min\\{\delta(S),\delta(T)\\}}{\max\\{\delta(S),\delta(T)\\}}$
$\displaystyle=\operatorname{NCD}_{\delta}(S,T)+\frac{\varepsilon^{\prime}}{1-\varepsilon^{\prime}}\cdot\Big{(}\frac{\delta(S,T)+\min\\{\delta(S),\delta(T)\\}}{\max\\{\delta(S),\delta(T)\\}}+\operatorname{NCD}_{\delta}(S,T)\Big{)}$
$\displaystyle\leq\operatorname{NCD}_{\delta}(S,T)+\frac{4\varepsilon^{\prime}}{1-\varepsilon^{\prime}}$
$\displaystyle\leq\operatorname{NCD}_{\delta}(S,T)+\varepsilon$

again using the sub-additivity of $\delta$ from Lemma 3.1, the fact that
$\operatorname{NCD}_{\delta}(S,T)\leq 1$ according to Corollary 3.2, the
definition of $\varepsilon^{\prime}$ and the assumption that $\varepsilon<1$.
∎

See 4.2

###### Proof.

We first observe that $\lceil\alpha^{\lfloor\log_{\alpha}n\rfloor}\rceil\leq
n$. To see this, assume otherwise, i.e., that
$\alpha^{\lfloor\log_{\alpha}n\rfloor}=n+x$ for some $x>0$. Then
$x=\alpha^{\lfloor\log_{\alpha}n\rfloor}-n\leq\alpha^{\log_{\alpha}n}-n=0$,
contradicting the assumption that $x>0$. It follows that $A\subseteq[n]$. Now,
for the upper bound notice that
$\tilde{\delta}\leq\max\\{(1+\varepsilon^{\prime})d_{k}/k:k\in
A\\}\leq(1+\varepsilon^{\prime})\cdot\delta\leq(1+\varepsilon)\cdot\delta$.

For the lower bound, let $\hat{k}\in[n]$ be such that
$\delta=d_{\hat{k}}/\hat{k}$ and let $i$ be such that
$\lceil\alpha^{i-1}\rceil\leq\hat{k}\leq\lceil\alpha^{i}\rceil$. Notice that
obviously $\lceil\alpha^{i-1}\rceil\in A$, but also $\lceil\alpha^{i}\rceil\in
A$ as $\alpha^{i}\leq\alpha\hat{k}\leq\alpha n/2\leq n$ by Observation A.2. We
now distinguish two cases: (1)
$\lceil\alpha^{i}\rceil=\lceil\alpha^{i-1}\rceil+1$ and (2)
$\lceil\alpha^{i}\rceil\geq\lceil\alpha^{i-1}\rceil+2$. In case (1), we get
that $\hat{k}\in\\{\lceil\alpha^{i-1}\rceil,\lceil\alpha^{i}\rceil\\}\subseteq
A$ and consequently
$\tilde{\delta}\geq\max\\{(1-\varepsilon^{\prime})d_{k}/k:k\in
A\\}=(1-\varepsilon^{\prime})\cdot\delta\geq(1-\varepsilon)\cdot\delta$. In
case (2), it holds that

$\displaystyle\alpha^{i-1}\cdot\varepsilon^{\prime}=\alpha^{i}-\alpha^{i-1}\geq\lceil\alpha^{i}\rceil-1-\alpha^{i-1}\geq\lceil\alpha^{i-1}\rceil+1-\alpha^{i-1}\geq
1.$ (1)

Now let $\beta:=\lceil\alpha^{i}\rceil-\hat{k}$. We note that $d_{j+1}\geq
d_{j}-1$ for every $j$, as every distinct length-$j$ substring other than
possibly $S[n-j+1..n]$ gives at least one distinct length-$j+1$ substring.
Applying the same observation iteratively yields
$d_{\lceil\alpha^{i}\rceil}\geq d_{\hat{k}}-\beta$. Hence

$\displaystyle\tilde{\delta}\geq(1-\varepsilon^{\prime})\cdot\frac{d_{\lceil\alpha^{i}\rceil}}{\lceil\alpha^{i}\rceil}\geq(1-\varepsilon^{\prime})\cdot\frac{d_{\hat{k}}-\beta}{\hat{k}+\beta}=(1-\varepsilon^{\prime})\delta\cdot\frac{1-\frac{\beta}{d_{\hat{k}}}}{1+\frac{\beta}{{\hat{k}}}}\geq(1-\varepsilon^{\prime})\delta\cdot\frac{1-\frac{\beta}{2{\hat{k}}}}{1+\frac{\beta}{{\hat{k}}}},$

where we used that $\delta(S)=d_{\hat{k}}/\hat{k}\geq 2$ in the last step. We
can now upper bound $\beta$ by
$\lceil\alpha^{i}\rceil-\lceil\alpha^{i-1}\rceil\leq\alpha^{i}+1-\alpha^{i-1}=\alpha^{i-1}\cdot\varepsilon^{\prime}+1\leq\hat{k}\varepsilon^{\prime}+1$.
This yields

$\hat{\delta}\geq(1-\varepsilon^{\prime})\delta\cdot\frac{1-\frac{\varepsilon^{\prime}}{2}-\frac{1}{2\hat{k}}}{1+\varepsilon^{\prime}+\frac{1}{\hat{k}}}\geq(1-\varepsilon^{\prime})\delta\cdot\frac{1-\varepsilon^{\prime}}{1+2\varepsilon^{\prime}}\geq(1-\varepsilon)\cdot\delta,$

where the second inequality uses that $\hat{k}\geq 1/\varepsilon^{\prime}$
following from (1) and the last inequality uses the definition of
$\varepsilon^{\prime}$. ∎

##  Appendix B Details on Bookmarking the RLBWT

We show how to extract $S[|S|-k+1]$ from the RLBWT, for any of the sampled
lengths $k$. See also the example in Figure 1. We show how to initialize and
update (upon character extensions of the stream) an index (bookmark) $j$ such
that $BWT[j]=S[|S|-k+1]$. This allows us retrieving $S[|S|-k+1]$ in
$O(\log|S|)\subseteq O(\log n)$ time with a random access operation $BWT[j]$
on the RLBWT data structure.

We first discuss how to initialize the bookmark $j$ as soon as the stream’s
length becomes $S=k+1$ (before that, the window of the last $k$ characters is
not completely filled). The initialization works by setting $j=BWT.LF(i)$,
where $i$ is the position such that $BWT[i]=\$$. Since the LF mapping on the
BWT of the reversed stream corresponds to advancing one position in the
stream, it is easy to see that, after this operation, it holds
$BWT[j]=S[|S|-k+1]$. See Figure 1 for an example.

Suppose we are storing the bookmark $j$ such that $BWT[j]=S[|S|-k+1]$. We now
show how to update $j$ when a new character $a$ arrives; let $S^{\prime}=Sa$
be the updated stream. Our goal is to modify $j$ so that
$BWT[j]=S^{\prime}[|S^{\prime}|-k+1]$ holds. The observation is that, upon the
extension of the stream by one character $a$, the algorithm of [19] modifies
the $BWT$ as follows: letting $i$ being the index such that $BWT[i]=\$$, the
algorithm (1) replaces $BWT[i]\leftarrow a$, and (2) inserts $ in the position
$i^{\prime}$ corresponding to the lexicographic rank of the new reversed
stream $(Sa)^{R}$ (position $i^{\prime}$ is computed in $O(\log|S|)$ time
using basic operations on the RLBWT, see [19] and Example 1): the new BWT
becomes $BWT\leftarrow BWT[1,i^{\prime}-1]\cdot\$\cdot BWT[i^{\prime}+1,|S|]$.
If $j<i^{\prime}$ (i.e. $ is inserted after position $j$), then after these
modification we have that $BWT[j]=S^{\prime}[|S^{\prime}|-k]$; if, on the
other hand, $j\geq i^{\prime}$ (i.e. $ is inserted before position $j$), then
we increment $j$ as $j\leftarrow j+1$, and $BWT[j]=S^{\prime}[|S^{\prime}|-k]$
holds also in this case. Finally, we need to “advance” $j$ by one position on
the stream; this operation corresponds to one LF mapping step on the BWT:
$j\leftarrow BWT.LF(j)$ ($O(\log|S|)$ time). After these operations, we
finally have that $BWT[j]=S^{\prime}[|S^{\prime}|-k+1]$.

##  Appendix C Detailed Experimental Results

###  C.1 Estimation of $d_{k}$.

As mentioned above, there are two types of errors in the computation of the
approximation $\tilde{\delta}$ of $\delta$: (1) the error obtained when
approximating the values $d_{k}$ by $\tilde{d}_{k}$ and (2) the error due to
the restriction of the string’s offsets $[n]$ to the “sampled set” $A$. The
error of type (1) itself has two sources, namely (1.1) errors due to
collisions when computing fingerprints with Rabin’s hash function and (1.2)
errors due to the count-distinct sketch when applied to the fingerprints. We
experimentally evaluated the error of type (1.1) and (1.2) as follows. For the
Pizza&Chili repetitive corpus, we compute the exact values of $d_{k}$ and
their estimated values $\tilde{d}_{k}$ for $k\in\\{2^{i}:0\leq i\leq 7\\}$. We
observe that the error (1.1) caused by collisions in Rabin’s fingerprint were
negligible; the error in the ratio of the distinct number of fingerprints and
the actual number of distinct substrings was less than 0.01%. The error due to
the count-distinct sketch (1.2) was dependent on its parameter: the number of
registers used for estimation. It is worth noting that the number of registers
does not affect the time complexity when updating sketches, but only affects
the space usage by a constant factor (and the time to compute the actual
estimation at the end, which is negligible). When more than $2^{14}$ registers
were used for count-distinct sketches, the maximum relative error was observed
to be below $2\%$, and the average error on $d_{k}$ was below $0.5\%$; see
Table 2.

The number of registers | $2^{10}$ | $2^{12}$ | $2^{14}$ | $2^{16}$
---|---|---|---|---
Maximum | 0.0641 | 0.0313 | 0.0186 | 0.0077
Average | 0.0204 | 0.0089 | 0.0041 | 0.0017
Table 2: The relative error measured in estimating $d_{k}$

###  C.2 Experiments on phylogenetic tree reconstruction.

To show similar behavior of $\operatorname{NCD}_{\mathtt{xz}}$,
$\operatorname{NCD}_{\delta}$, and $\operatorname{NCD}_{\tilde{\delta}}$, we
conducted experiments on phylogenetic tree reconstruction using Gene-Trees
dataset from AF project777https://afproject.org/. It contains 11 groups of
sequences (651 sequences in total) where each group yields a tree. We
constructed 11 phylogenetic trees (i.e. one tree for each group) for each of
the NCD measures, and compare the constructed trees by measuring the
normalized Robinson-Foulds (nRF) distance, a widely-used distance measure for
this purpose. The distance tends to 0 as trees become similar, and tends to 1
when comparing with a random tree. The average nRF between
$\operatorname{NCD}_{\mathtt{xz}}$ and $\operatorname{NCD}_{\delta}$ was
measured as 0.2, indicating similar trees were reconstructed. The average nRF
between $\operatorname{NCD}_{\delta}$ and
$\operatorname{NCD}_{\tilde{\delta}}$ ranges from 0.115 to 0.250 depending on
the parameters. For ease of interpretation of these values, we depict two
similar phylogenetic trees with nRF=0.194 in Figure 3, which is an actual
example of reconstructed trees using NCD with $\tilde{\delta}$ and
$\mathtt{xz}$.

forked edges, /tikz/every label/.append style=font=, before typesetting nodes=
delay= where content=coordinate, , where n children=0tier=terminus,
label/.process=Owcontentright:#1, content=, , for tree= grow’=0, s sep’+=2pt,
l=1pt, l sep=9pt, ,
[[CITE1_MACMU][[CITE1_PANTR][[[CITE1_BOVIN][CITE1_CANFA]][[[CITE1_MOUSE][CITE1_RAT]][[CITE1_ORNAN][[CITE1_DANRE][[[[[[[[[[CITE2_HUMAN][CITE2_PANTR]][CITE2_BOVIN]][[CITE2_MOUSE][CITE2_RAT]]][CITE2_MONDO]][CITE2_MACMU]][CITE2_CHICK]][CITE2_XENTR]][CITE2_DANRE]][[[[CITE4_CHICK][[CITE4a_DANRE][[CITE4a_TAKRU][[CITE4b_DANRE][CITE4b_TAKRU]]]]][CITE4_XENTR]][[[[[[[CITE4_HUMAN][CITE4_PANTR]][CITE4_MACMU]][CITE4_BOVIN]][[CITE4_MOUSE][CITE4_RAT]]][CITE4_MONDO]][[CITED_BRAFL][CITED_NEMVE]]]]]]]]]][CITE1_HUMAN]]
forked edges, /tikz/every label/.append style=font=, before typesetting nodes=
delay= where content=coordinate, , where n children=0tier=terminus,
label/.process=Owcontentright:#1, content=, , for tree= grow’=0, s sep’+=2pt,
l=1pt, l sep=9pt, ,
[[CITE1_PANTR][[CITE1_MACMU][[[CITE1_BOVIN][CITE1_CANFA]][[[CITE1_MOUSE][CITE1_RAT]][[[CITE1_ORNAN][CITE1_DANRE]]
[
[[[[[[[[[[[CITE2_HUMAN][CITE2_PANTR]][CITE2_BOVIN]][[CITE2_MOUSE][CITE2_RAT]]][CITE2_MONDO]][CITE2_MACMU]][CITE2_CHICK]][CITE2_XENTR]][CITE2_DANRE]][CITED_BRAFL]][[[CITE4_CHICK][[[CITE4a_DANRE][[CITE4b_DANRE][CITE4b_TAKRU]]][CITE4a_TAKRU]]][CITE4_XENTR]]]
[[[[[[[CITE4_HUMAN][CITE4_PANTR]][CITE4_MACMU]][CITE4_BOVIN]][[CITE4_MOUSE][CITE4_RAT]]][CITE4_MONDO]][CITED_NEMVE]]
] ]]]][CITE1_HUMAN]]

Figure 3: Two similar phylogenetic trees constructed using normalized
compression distance (NCD) with estimated normalized substring complexity
$\tilde{\delta}$ (left) and a popular compression software $\mathtt{xz}$
(right). Normalized Robinson-Foulds distance is 0.194.

#### Running Time.

To construct a phylogenetic tree from a sequence set, we usually need to
compute all-pair distances. When sequences are long, computing NCDs can be
quite costly because we need to compress concatenated sequences for all pairs
of sequence in the input set. On the other hand, our sketching can be more
efficient because we only need to compute sketches for each sequence, then
merging sketches can be done very quickly compared to processing the entire
sequences all over again. To demonstrate this, we measure the running time for
computing all-pair NCDs on E.coli dataset from AF Project that consists of 29
sequences of average length 81,588. Computing all-pair NCDs with xz and the
exact $\delta$ took about 42 and 24 minutes. On the other hand, our sketching
method took only 3 minutes.


Precision measurements of jet and photon production at the ATLAS experiment
111Copyright 2022 CERN for the benefit of the ATLAS Collaboration. CC-BY-4.0
license.

Giuseppe Callea on behalf of the ATLAS collaboration1$\star$

1 University of Glasgow

*<EMAIL_ADDRESS>
51st International Symposium on Multiparticle Dynamics (ISMD2022)
Pitlochry, Scottish Highlands, 1-5 August 2022 10.21468/SciPostPhysProc.?

## Abstract

The production of jets and prompt isolated photons at hadron colliders
provides stringent tests of perturbative QCD. We present the latest
measurements of photon+jets and diphoton production using proton-proton
collision data collected by the ATLAS experiment at the LHC. The measurements
are compared to state-of-the-art NLO and NNLO predictions.

##  1 Introduction

The production of prompt photons in association with hadronic jets in proton-
proton collisions, $pp\to\gamma+\mathrm{jet}+X$, provides a testing ground for
perturbative Quantum Chromodynamics (pQCD) in a cleaner enviroment than in jet
production, since the colourless photon originates directly from the hard
interaction. The measurement of the angular correlations between the photon
and the jet can be used to probe the dynamics of the hard-scattering process.
Diphoton production offers a clean final state for the study of the properties
of the Higgs boson and a possible window into new physics phenomena. Prompt
photon production can proceed via ‘direct’ production, where the photon
originates from the hard scatter, or the ‘fragmentation’ process where it is
produced in the fragmentation of a parton with high $p_{\mathrm{T}}$ [1]. This
contribution summarises the studies on the photon + two jet [2] and photon
pair productions [3] with the ATLAS detector [4] at the Large Hadron Collider
(LHC).

##  2 Photon + 2 jets production

The production of prompt photons with two associated jets provides a solid
testing ground for pQCD by looking at the angular correlations between the
photon and each of the jets and between the two jets. In addition,
measurements of the invariant-mass distributions of the dijet system and the
$\gamma$-jet-jet system are sensitive to the dynamics of the hard interaction.
A comprehensive study of the observables describing this final state is of
relevance for the development of pQCD calculations as well as for the tuning
of Monte Carlo (MC) models. These studies were performed using a dataset with
an integrated luminosity of 36.1 fb-1, collected with the ATLAS detector
during Run 2.
Photons are required to pass isolation and identification requirements and
have $E_{\mathrm{T}}^{\gamma}$ > 150 GeV. Jets are reconstructed with the
anti-$k_{t}$ algorithm with radius parameter R=0.4 and are required to have
$p_{\mathrm{T}}$ > 100 GeV. A ‘fragmentation-enriched’ sample is selected by
requiring $E_{\mathrm{T}}^{\gamma}$ < $p_{\mathrm{T}}^{\mathrm{jet}2}$, where
$p_{\mathrm{T}}^{\mathrm{jet}2}$ is the transverse momentum of the sub-leading
jet. A ‘direct-enriched’ sample is selected with the requirement
$E_{\mathrm{T}}^{\gamma}$ > $p_{\mathrm{T}}^{\mathrm{jet}1}$, where
$p_{\mathrm{T}}^{\mathrm{jet}1}$ is the $p_{\mathrm{T}}$ of the leading jet.
The tree-level plus parton-shower predictions from leading order (LO) Sherpa
[5] and Pythia8 [6] as well as the next-to-leading-order (NLO) QCD predictions
from Sherpa are compared with the measurements. Figure 1 shows the
differential cross-section of the photon + 2 jet production as a function of
$E_{\mathrm{T}}^{\gamma}$ and the invariant mass of the $\gamma$-jet-jet
system, $m_{\gamma\textrm{-}\mathrm{jet}\textrm{-}\mathrm{jet}}$ . The NLO QCD
predictions of Sherpa describe the data adequately in shape and normalisation
except for regions of phase space such as those with high values of the
invariant mass, where the predictions overestimate the data. The Sherpa LO
predictions are in good agreement with the data.

(a)

(b)

Figure 1:  Differential cross-sections of the photon + 2 jet production
measured as a function of Figure $E_{\mathrm{T}}^{\gamma}$ (a) and
$m_{\gamma\textrm{-}\mathrm{jet}\textrm{-}\mathrm{jet}}$ (b) for both direct
and fragmentation process and compared to the theoretical MC predictions [2].

##  3 Diphoton production

Measurements of a photon pair production were studied with the full Run 2
dataset collected with the ATLAS detector with an integrated luminosity of 139
fb-1. At least two photon candidates are required per event with transverse
momentum above 40 GeV and 30 GeV, respectively. The photon pseudorapidity is
required to be $|\eta^{\gamma}|$ < 2.37, excluding the transition region
between the barrel and the endcap (1.37 < $|\eta^{\gamma}|$ < 1.52). In this
acceptance region, the high granularity of the calorimeter system allows the
efficient identification of photons. The unfolded data were compared to
several state-of-the-art predictions: Fixed-order NNLO with NNLOJet framework
[7] to obtain the diphoton predictions up to NNLO QCD accuracy; Fixed-order
NLO with Diphox [8]; Multi-leg Sherpa 2.2 generator with matrix elements for
the $pp\to\gamma\gamma$ \+ 0,1 jet process at NLO and $pp\to\gamma\gamma$ \+
2,3 jet process at LO, matched and merged with the parton shower using the
MEPS@NLO prescription [9]. Differential cross sections are measured as
functions of the main variables of the photon pair system, together with the
transverse component of $p_{\mathrm{T,\gamma\gamma}}$ with respect to the
thrust axis, $a_{\mathrm{T}}$. Figure 2a and 2b show the differential cross
section for the photon pair production as a function of
$p_{\mathrm{T}}^{\gamma,1}$ and $m_{\gamma\gamma}$, respectively. Only the
NNLO, as implemented by NNLOJET, and Sherpa NNLO predictions provide a
satisfactory description of the data. The fixed-order NNLO calculation has a
higher formal precision than Sherpa and leads to lower theoretical
uncertainties in the perturbative QCD governed regions. Fixed-order
predictions of DIPHOX and NNLOJET are not expected to be valid in regions
dominated by multiple collinear and soft QCD emissions. The $m_{\gamma\gamma}$
distribution is governed by the $p_{\mathrm{T}}$ requirements on the
individual photons of 40 GeV and 30 GeV, respectively. The
$\gamma\gamma$+multi-jet configurations are dominant in the $m_{\gamma\gamma}$
< 70 GeV region. The integrated fiducial cross sections are shown in Figure
2c. NNLOJet and Sherpa MEPS@NLO are in better agreement with respect to the
other predictions considered.

(a)

(b)

(c)

Figure 2:  Differential cross-sections of the photon pair production measured
as a function of $p_{\mathrm{T}}^{\gamma,1}$ (a) and $m_{\gamma\gamma}$ (b)
for both direct and fragmentation processes and compared to the theoretical MC
preditions. Figure (c) shows the integrated fiducial cross section comparisons
between the ATLAS data and the theoretical predictions [3].

##  4 Conclusion

In this contribution, high-precision measurements involving the production of
photons and jets at the centre-of-mass energy of 13 TeV are presented. All the
measurements are in good agreement with the pQCD predictions within the
theoretical uncertainties. The Sherpa multi-leg merged predictions lead to
larger uncertainties, but they are in better agreement with the data in all
regions. These studies will provide valuable physics inputs to PDF fits and MC
generators.

## References

  * [1] T. Pietrycki and A. Szczurek, Photon-jet correlations in $pp$ and $p\bar{p}$ collisions, PRD 76 (2007) 034003, 10.1103/PhysRevD.76.034003.
  * [2] ATLAS collaboration, Measurement of isolated-photon plus two-jet production in pp collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector, JHEP 03 (2020) 179, 10.1007/JHEP03(2020)179.
  * [3] ATLAS collaboration, Measurement of the production cross section of pairs of isolated photons in $pp$ collisions at 13 TeV with the ATLAS detector, JHEP 11 (2021) 169, 10.1007/JHEP11(2021)169.
  * [4] ATLAS collaboration, The ATLAS Experiment at the CERN Large Hadron Collider, JINST 3 (2008) S08033, 10.1088/1748-0221/3/08/S08003.
  * [5] E. Bothmann et al., Event generation with Sherpa 2.2, SciPost Phys. 7 (2019) 034, 10.21468/SciPostPhys.7.3.034.
  * [6] T. Sjöstrand, S. Mrenna and P. Skands, A brief introduction to PYTHIA 8.1, Comput. Phys. Commun. 178 (2008) 852, 10.1016/j.cpc.2008.01.036.
  * [7] T. Gehrmann et al., jet cross sections and transverse momentum distributions with NNLOJET, PoS RADCOR2017 (2018) 074, 10.48550/arXiv.1801.06415.
  * [8] T. Binoth, J. P. Guillet, E. Pilon and M. Werlen, A full next-to-leading order study of direct photon pair production in hadronic collisions, EPJC 16 (2000) 311, 10.1007/s100520050024.
  * [9] S. Hoeche et al, A critical appraisal of NLO+PS matching methods, JHEP09(2012)049, 10.1007/JHEP09%282012%29049.


# Dynamics of temperature from a relational approach

Sebastian Gemsheim<EMAIL_ADDRESS>Jan M. Rost  Max Planck Institute
for the Physics of Complex Systems, Nöthnitzer Straße 38, D-01187 Dresden,
Germany

# Statistical mechanics from relational complex time with a pure state

Sebastian Gemsheim<EMAIL_ADDRESS>Jan M. Rost  Max Planck Institute
for the Physics of Complex Systems, Nöthnitzer Straße 38, D-01187 Dresden,
Germany

###### Abstract

Thermodynamics and its quantum counterpart are traditionally described with
statistical ensembles. Canonical typicality has related statistical mechanics
for a system to ensembles of global energy eigenstates of system and its
environment analyzing their cardinality. We show that the canonical density
for a system emerges from a maximally entangled global state of system and
environment through relational complex time evolution between system and
environment without the need to maximize the entropy or to count states.

Canonical typicality has proven to be a powerful concept to understand how
statistical mechanics emerges for a (small) system $\mathscr{S}$ singled out
from its environment $\mathscr{E}$ where environment and system together are
in a pure energy eigenstate of the global Hilbert space [1, 2]. This is
achieved by analyzing what is "typical" about the ensemble of excessively many
possibilities to realize such an eigenstate in a small energy interval
$[E,E+\delta E]$ and assessing the exponentially growing number of typical
realizations in ensembles of increasing energy $E$. The approach carries the
spirit of statistical mechanics which deals with many degrees of freedom and
therefore naturally with ensembles. Yet, quantum mechanics introduces already
for a single degree of freedom a probabilistic perspective with the
wavefunction associated to its dynamics. Hence, one could hope to capitalize
on this fundamental quantum property to derive statistical mechanics from a
single pure state. The fact that, compared to classical mechanics, quantum
mechanics offers with entanglement another resource for variety already
encoded in a single state renders the existence of such a state as the source
of statistical mechanics more likely. We will show that statistical mechanics
is contained in a single, maximally entangled global eigenstate (MES) of
system and environment without the need to consider ensembles of global states
at different energies [1, 2] or to count degenerate states of the system [3].
Instead, we extend the relational time approach originally conceived by Page &
Wootters [4] to propagate density matrices in imaginary relational time. This
will allow us to conceive the fundamental thermodynamic relation from the
density matrix of a maximally entangled global state for any desired imaginary
time $\gamma$, which will turn out to represent temperature in the well known
way, $(k_{B}T)^{-1}=\gamma$, if infinite temperature $\gamma=0$ is linked to
the MES. Hence, relational complex time evolution turns out to be the element
which completes the quest for deriving statistical mechanics from a pure
quantum state without the need to resort to additional statistical
considerations. The result confirms the conjecture that temperature, similarly
as time, has a relational root [5] and provides a physical justification for
the striking connection between dynamics and thermal physics [6].

Relational time emerges from a splitting of the global Hamiltonian
$\hat{H}_{\mathrm{tot}}{=}\hat{H}\otimes\hat{\mathds{1}}_{\mathscr{E}}+\hat{\mathds{1}}\otimes\hat{H}_{\mathscr{E}}+\hat{V}$
into system $\hat{H}$ and environment $\hat{H}_{\mathscr{E}}$ interacting with
each other through $\hat{V}$. Since the global state $\Psi$ fullfills
$(\hat{H}_{\mathrm{tot}}-E_{\mathrm{tot}})\ket*{\Psi\rangle\\!}=0$ one can
formulate the _invariance principle_ [7]

$e^{\Lambda(\hat{H}_{\mathrm{tot}}-E_{\mathrm{tot}})}\ket*{\Psi\rangle\\!}=\ket*{\Psi\rangle\\!}$
(1)

which does not only hold for all imaginary-valued symmetry parameters
$\Lambda=i\lambda$, but also for real ones $\Lambda=\gamma$, or in general for
complex $\Lambda$.

With (1) we can write for the density matrix
$P_{\Psi}=\outerproduct*{\Psi\rangle\\!}{\\!\langle\Psi}$ of the global state

$\hat{P}_{\Psi}=e^{\Lambda(\hat{H}_{\mathrm{tot}}-E_{\mathrm{tot}})}\hat{P}_{\Psi}e^{\Lambda^{*}(\hat{H}_{\mathrm{tot}}-E_{\mathrm{tot}})}\,.$
(2)

From now on we will assume in accordance with the original proposal for
relational time [4] that the interaction $V$ is small enough to be neglected.
In fact, this condition can be softened without loosing the thrust of the
relational time concept by only demanding that
$[\hat{H}_{\mathscr{E}},\hat{V}]=0$, implying that the environment’s
eigenstates are not changed by the interaction, as we will see later.

For now with $\hat{V}{=}0$, the evolution according to the Hamiltonians for
the environment
$\hat{h}_{\mathscr{E}}\equiv\hat{H}_{\mathscr{E}}-E_{\mathrm{tot}}$ and system
$\hat{H}$ in (2) factorizes,

$e^{\Lambda(\hat{H}_{\mathrm{tot}}-E_{\mathrm{tot}})}=e^{\Lambda\hat{h}_{\mathscr{E}}}{\color[rgb]{0,0,0}\otimes}e^{\Lambda\hat{H}}\,,$
(3)

and we can obtain from (2) the density matrix of the system evolved in complex
time $\Lambda$,

$\hat{\rho}(\Lambda)=\frac{1}{{Z}(\Lambda)}e^{-\Lambda\hat{H}}\hat{\rho}_{0}e^{-\Lambda^{*}\hat{H}}\,,$
(4a) where
$e^{-\Lambda\hat{H}}\hat{\rho}_{0}e^{-\Lambda^{*}\hat{H}}=\tr_{\mathscr{E}}[\hat{\rho}_{\mathscr{E}}(-\Lambda^{*})\hat{P}_{\Psi}]$
(4b) with
$\hat{\rho}_{\mathscr{E}}(\Lambda)=\,e^{-\Lambda\hat{h}_{\mathscr{E}}}\hat{\rho}_{\mathscr{E}}(0)e^{-\Lambda^{*}\hat{h}_{\mathscr{E}}}\,,$
(4c)

which follows from (2) after tracing over the environment. We have left room
for a conditioning density $\hat{\rho}_{\mathscr{E}}(0)$ of the environment,
such that
$\hat{\rho}(0)={Z}(0)^{-1}\tr_{\mathscr{E}}[\hat{\rho}_{\mathscr{E}}(0)P_{\Psi}]$.
The normalization
${Z}(\Lambda)=\langle\\!\langle\Psi|\hat{\rho}_{\mathscr{E}}({\color[rgb]{0,0,0}-\Lambda^{*}})|\Psi\rangle\\!\rangle$
ensures that $\tr\hat{\rho}(\Lambda)=1$.

Note that according to (4b) the evolution of the system density $\hat{\rho}$
with $\Lambda$ means evolution of the environment density with $-\Lambda^{*}$
illustrating the relational character of the evolution in system and
environment, respectively. For real and continuous $\Lambda=\gamma$,
differentiation of (4a) leads to

$d\hat{\rho}(\gamma)/d\gamma=-\Bigl{\\{}\hat{H}-\langle\hat{H}\rangle(\gamma),\hat{\rho}(\gamma)\Bigl{\\}}\,,$
(5)

the von Neumann equation for imaginary time $\gamma$ emerging from complex
relational time. The mean energy
$\langle\hat{H}\rangle=\tr(\hat{H}\hat{\rho})$ of the system in the anti-
commutator $\\{\cdot,\cdot\\}$ ensures the normalization. For real time
$\Lambda=i\lambda$, $\lambda\in\mathds{R}$, the standard von Neumann equation

$d\hat{\rho}({\color[rgb]{0,0,0}i\lambda})/d\lambda=-i[\hat{H},\hat{\rho}({\color[rgb]{0,0,0}i\lambda})]\,,$
(6)

for the time-dependence of the density matrix follows from (4a) by
differentiation. The conditioning environmental density
$\hat{\rho}_{\mathscr{E}}$ in
$\hat{\rho}_{0}=\tr_{\mathscr{E}}[\hat{\rho}_{\mathscr{E}}(0)P_{\Psi}]$ allows
one to specify different initial conditions with the global state
$\ket*{\Psi\rangle\\!}$ as discussed for wavefunctions in [7]. Both
differential equations hold of course for all global states satisfying (1).

Statistical mechanics is encoded in (4). Standard tracing over the environment
(implying $\hat{\rho}_{\mathscr{E}}(0)\propto\hat{\mathds{1}}_{\mathscr{E}}$)
with eigenstates $\ket*{J}_{\mathscr{E}}$ and eigenenergies $E_{J}$ of
$\hat{H}_{\mathscr{E}}$ in (4b) renders the system density diagonal. Denoting
$\ket{\phi_{J}}=\innerproduct{J}{\Psi\rangle\\!}$ and
$\Lambda=\gamma/2+i\lambda$ we get from (4b) using (4c)

$\displaystyle\hat{\rho}(\Lambda)$
$\displaystyle=\frac{1}{{Z}(\Lambda)}\sum_{J}\matrixelement{J}{e^{\Lambda\hat{h}_{\mathscr{E}}}\hat{P}_{\Psi}e^{\Lambda^{*}\hat{h}_{\mathscr{E}}}}{J}_{\mathscr{E}}$
$\displaystyle=\frac{1}{{Z}(\Lambda)}\sum_{{\color[rgb]{0,0,0}J}}e^{-\gamma\epsilon_{J}}\outerproduct{\phi_{J}}{\phi_{J}}$
(7) $\displaystyle=\frac{1}{{Z}(\Lambda)}e^{-\gamma\hat{H}}\hat{\rho}_{0}\,,$

where $\epsilon_{J}=E-E_{J}$ and the normalization ${Z}(\Lambda)$ absorbs the
constant of proportionality in
$\hat{\rho}_{\mathscr{E}}(0)\propto\hat{\mathds{1}}_{\mathscr{E}}$. The last
identity follows from the invariance principle and (3).

Obviously, (Dynamics of temperature from a relational approach) represents the
canonical density for all temperatures if one identifies $\gamma=1/(kT)$ and
if $\innerproduct{\phi_{J}}{\phi_{J^{\prime}}}=c\delta_{JJ^{\prime}}$, where
the constant $c$ once again can be absorbed in $Z(\Lambda)$ and therefore may
be set to unity, $c=1$. Then, evaluating the entropy
$S=-\tr(\hat{\rho}\ln\hat{\rho})$ with (Dynamics of temperature from a
relational approach) gives
$S({\color[rgb]{0,0,0}\gamma/2})=\gamma\langle\hat{H}\rangle({\color[rgb]{0,0,0}\gamma/2})+\ln{Z}({\color[rgb]{0,0,0}\gamma/2})$
implying the thermodynamic relation $dS=\gamma d\langle\hat{H}\rangle$.

Note that we have derived the canonical density for the system solely by the
relational connection to the global state apart from requiring the
$\ket{\phi_{J}}$ to be orthonormal. This requirement has a physical reason:
With a system at infinite temperature $\gamma{=}0$ we associate minimal
structure implying in the relational context that the system states
$\ket{\phi_{J}}$ with their different energies $\epsilon_{J}$ should be
contained in the global state with equal probability
$\innerproduct{\phi_{J}}{\phi_{J}}=\langle\\!\langle\Psi|J\rangle\langle
J|\Psi\rangle\\!\rangle=1$ for all $J$. The global state with this property is
maximally entangled,

$\ket*{\Psi_{M}\rangle\\!}=\sum_{J}{\color[rgb]{0,0,0}\ket*{J}_{\mathscr{E}}\otimes\innerproduct{J}{\Psi_{M}\rangle\\!}}=\sum_{J}\ket*{J}_{\mathscr{E}}\otimes\ket{\phi_{J}}\,,$
(8)

where the $\\{\ket{\phi_{J}}\\}$ form an orthonormal basis for the system
Hamiltonian with eigenenergies $\epsilon_{J}$. Note, that thanks to the MES
and with
$\langle\hat{H}_{\mathscr{E}}\rangle_{\mathscr{E}}=\tr_{\mathscr{E}}(\hat{H}_{\mathscr{E}}\hat{\rho}_{\mathscr{E}})$,

$E=\langle\hat{H}\rangle(\Lambda)+\langle\hat{H}_{\mathscr{E}}\rangle_{\mathscr{E}}(\Lambda)$
(9)

holds for all complex times $\Lambda$, in particular for all temperatures. The
link of the MES to infinite temperature through complex relational time
motivates its role for statistical mechanics which is given in [3] on very
general grounds of symmetry considerations for a global pure state describing
system and environment.

This remarkably direct quantum route to the fundamental relation of
thermodynamics through relational imaginary time propagation of a maximally
entangled global state bears some interesting properties of relational complex
time, we would like to point out in the following.

  * (A)

The energy constraint (1) has the consequence, that

$[\hat{\rho}_{\mathscr{E}},\hat{h}_{\mathscr{E}}]=0\Leftrightarrow[\hat{\rho},\hat{H}]=0\,.$
(10)

In other words, if the conditioning density $\hat{\rho}_{\mathscr{E}}$ is
(block) diagonal in the energy eigenstates of the environment, the system’s
density matrix $\hat{\rho}$ commutes with the system Hamiltonian, rendering
$\hat{\rho}$ stationary in real time, see (6). As the imaginary time evolution
is ruled by the anti-commutator, the system’s density changes in imaginary
time, but only through the exponential energy factors, the coefficients in the
density remain constant. Since simply tracing over the environment implies
$\hat{\rho}_{\mathscr{E}}(0)\propto\hat{\mathds{1}}_{\mathscr{E}}$ diagonal,
this is also true for the imaginary time evolution giving rise to statistical
mechanics as derived: The evolution relates via the von Neumann equation in
imaginary time the system at different temperatures, maintaining equilibrium
(that is stationarity in real time), thereby explicating the "peaceful
coexistence of thermal equilibrium and emergence of time" as put forward in
[8].

  * (B)

As a pure state the MES carries vanishing entropy $S_{\mathrm{tot}}{=}0$. A
density matrix of the system singled out from the MES by tracing over the
environment,
$\hat{\rho}={\color[rgb]{0,0,0}D^{-1}}\tr_{\mathscr{E}}[{\color[rgb]{0,0,0}\hat{P}_{\Psi}}]={\color[rgb]{0,0,0}D^{-1}}\sum_{J=1}^{D}\outerproduct{{\color[rgb]{0,0,0}\phi_{J}}}{{\color[rgb]{0,0,0}\phi_{J}}}$,
is maximally mixed with (maximal) entropy, $S=\ln{D}$, where $D$ is the
dimension of the system’s Hilbert space. This property characterizes the
maximal deficit of knowledge or degree of in-determinism, which is typically
associated with a system’s state at infinite temperature.

The mathematical connection of a maximally determined global MES with a
maximally mixed (in-determined) state of the system originates from the way,
how information about the system is encoded in the MES, namely by exclusively
and evenly linking system wave functions to those of the environment. Hence,
giving up the link through tracing over the environment leaves minimal
information about the system. On the other hand, a global state which is not
an MES is in general a Schmidt state
$\ket*{\Psi\rangle\\!}=\sum_{J}a_{J}\ket*{J}_{\mathscr{E}}\otimes\ket{\phi_{J}}$,
$a_{J}\in\mathds{R}$. Links within (sub-)systems, represented by coefficients
$a_{J}\neq 1$ in the global state are preserved upon cutting the links to the
environment by tracing over it, and the resulting system state is not
maximally mixed.

In an experimental situation the global system designed may not be in a
maximally entangled or close by state, but only in a Schmidt state. In this
case it is still possible to mimic thermodynamics for the system via a
suitable conditioning density $\hat{\rho}_{\mathscr{E}}$ of the environment.
With
$\hat{\rho}_{\mathscr{E}}(0)=\sum_{J}{\color[rgb]{0,0,0}p_{J}}\outerproduct{J}{J}_{\mathscr{E}}$
one gets explicitly

$\displaystyle\hat{\rho}(\Lambda)=$
$\displaystyle\frac{1}{{Z}({\color[rgb]{0,0,0}\Lambda})}e^{-\gamma\hat{H}}\hat{\rho}(0)$
$\displaystyle=$
$\displaystyle\frac{1}{{Z}({\color[rgb]{0,0,0}\Lambda})}\sum_{J}e^{-\gamma\epsilon_{J}}{\color[rgb]{0,0,0}p_{J}}a_{J}^{2}\outerproduct{\phi_{J}}{\phi_{J}}\,,$
(11)

where $\Lambda=\gamma/2+i\lambda$. A conditioning density
$\hat{\rho}_{\mathscr{E}}$ with ${\color[rgb]{0,0,0}p_{J}}=a_{J}^{-2}$ gives
rise to a canonical density for the system and therefore statistical
mechanics. Such a setup is conceivable using cavities and quantum control. It
may be even more interesting to study with this kind of setup deviations from
statistical mechanics with densities of the form of (11). Another step of
complication would be to include real time dynamics if $\hat{\rho}$ does not
commute with $\hat{H}$, directly relevant for eigenstate thermalization [9,
10]. Of course this complicates, the expression for entropy which is also true
for explicit interaction $V\neq 0$ in $\hat{H}_{\mathrm{tot}}$.

Partially, the effect of interaction is included through the entanglement in
the global state: Without interaction the global state could be a simple
product state. Preserving the general outcome, an interaction $\hat{V}$ in (1)
can be included which commutes with the environmental Hamiltonian,
$[\hat{H}_{\mathscr{E}},\hat{V}]=0$, or due to the Jacobi identity
111$[\hat{H}_{\mathscr{E}},[\hat{V},\hat{P}_{J}]]+[\hat{V},[\hat{P}_{J},\hat{H}_{\mathscr{E}}]]+[\hat{P}_{J},[\hat{H}_{\mathscr{E}},\hat{V}]]=0$,
where $\hat{P}_{J}=\outerproduct{J}{J}_{\mathscr{E}}$. equivalently
$[V,{\color[rgb]{0,0,0}\outerproduct{J}{J}_{\mathscr{E}}}]=0$. This means that
the environment is big enough such that its eigenfunctions are not modified by
$\hat{V}$. The result (Dynamics of temperature from a relational approach)
remains the same, only that now the eigenfunctions $\ket{\phi_{J}}$ and
eigenenergies $\epsilon_{J}$ fulfill the Schrödinger equation
$(\hat{H}+{\color[rgb]{0,0,0}\matrixelement{J}{\hat{V}}{J}_{\mathscr{E}}}-\epsilon_{J})\ket{\phi_{J}}=0$.

In summary, we have shown from a relational dynamics point of view how
statistical physics of a system is encoded in maximally entangled states of
the system and its environment. The observation underscores their pivotal role
for fundamental quantum principles, as these states are also key in the
derivation of Born’s rule with envariance [12]. Maximally entangled states
have the least structure and are evenly distributed over system and
environment, corroborating the situation at infinite temperature. Without the
need of ensemble considerations such as canonical typicality, but by
capitalizing on the probabilistic nature of a quantum state itself,
statistical physics for the system at finite temperature emerges from
propagation of the maximally entangled state in relational complex time.

We thank A. Eisfeld and F. Fritzsch for helpful discussions.

## References

  * Goldstein  _et al._ [2006] S. Goldstein, J. L. Lebowitz, R. Tumulka, and N. Zanghì, Canonical typicality, Phys. Rev. Lett. 96, 050403 (2006).
  * Popescu  _et al._ [2006] S. Popescu, A. J. Short, and A. Winter, Entanglement and the foundations of statistical mechanics, Nature Physics 2, 754–758 (2006)
  * Deffner and Zurek [2016] S. Deffner and W. H. Zurek, Foundations of statistical mechanics from symmetries of entanglement, New Journal of Physics 18, 063013 (2016).
  * Page and Wootters [1983] D. N. Page and W. K. Wootters, Evolution without evolution: Dynamics described by stationary observables, Physical Review D 27, 2885 (1983).
  * Vedral [2017] V. Vedral, Time, (inverse) temperature and cosmological inflation as entanglement, in  _Time in Physics_ (Springer International Publishing, 2017) pp. 27–42.
  * Feynman and Hibbs [2010] R. P. Feynman and A. R. Hibbs,  _Quantum mechanics and path integrals_ , Dover Books on Physics (Dover Publications, Mineola, NY, 2010).
  * Gemsheim and Rost [2023] S. Gemsheim and J. M. Rost, Emergence of time from quantum interaction with the environment, Phys. Rev. Lett. 131, 140202, (2023)
  * Favalli and Smerzi [2022] T. Favalli and A. Smerzi, Peaceful coexistence of thermal equilibrium and the emergence of time, Phys. Rev. D 105, 023525 (2022).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Advances in Physics 65, 239–362 (2016)
  * Deutsch [2018] J. M. Deutsch, Eigenstate thermalization hypothesis, Reports on Progress in Physics 81, 082001 (2018).
  * Note [1] $[\hat{H}_{\mathcal{E}},[\hat{V},\hat{P}_{J}]]+[\hat{V},[\hat{P}_{J},\hat{H}_{\mathcal{E}}]]+[\hat{P}_{J},[\hat{H}_{\mathcal{E}},\hat{V}]]=0$, where ${\color[rgb]{0,0,0}\hat{P}_{J}=\outerproduct{J}{J}_{\mathscr{E}}}$.
  * Zurek [2003] W. H. Zurek, Environment-assisted invariance, entanglement, and probabilities in quantum physics, Phys. Rev. Lett. 90, 120404, (2003)