# Hausdorff’s forgotten proof that
almost all numbers are normal

Edmund Weitz111E-mail<EMAIL_ADDRESS>
Hamburg University of Applied Sciences

(January 31, 2021)

###### Abstract

In 1914, Felix Hausdorff published an elegant proof that almost all numbers
are simply normal in base 2. We generalize this proof to show that almost all
numbers are normal. The result is arguably the most elementary proof for this
theorem so far and should be accessible to undergraduates in their first year.

Mathematics Subject Classification (2020): 11K16

Keywords: normal numbers

##  1 Introduction

In 1909, Émile Borel [1] introduced normal numbers and proved that almost all
numbers are normal.222According to [2], the original proof was “unmendably
faulty,” though. Today, several different proofs of Borel’s theorem exist.
This includes some that are generally considered to be elementary, e.g., [6],
[7], and [3].

One of the earliest proofs can be found in Felix Hausdorff’s Grundzüge der
Mengenlehre [5] from 1914. But Hausdorff only proved that almost all numbers
are simply normal in base 2 and then claimed it would be “evident” that the
statement was true for other bases as well. He didn’t define normal numbers
and gave no indication how to prove a stronger version of his result. As we
will show in this article, Hausdorff’s argument isn’t hard to generalize,
although the way to do it might not be totally obvious either.

Anyway, to the author’s knowledge, nobody has picked up Hausdorff’s elegant
idea so far. [4] and [8] contain proofs which argue along similar lines but
require more technical finesse and are less direct.

This article is intended to be accessible to undergraduates at the beginning
of their studies and we thus won’t presuppose a lot of previous knowledge
except for basic combinatorics, basic calculus, and a bit of set theory (up to
the definition of countable). Everything else will be defined and proved,
including enough (informal) measure theory to state and prove the main
theorem.

##  2 “Almost all”

The idea of a measure is to assign non-negative numbers to sets (of real
numbers) in such a way that these numbers can intuitively be interpreted as
the sizes of the sets. Two obviously meaningful requirements for a measure are
that the empty set is assigned the measure zero and that the measure is
additive: the measure of the union of two (or finitely many) sets which are
mutually disjoint must be the sum of their measures. In order to be useful in
analysis, measures are actually required to be $\sigma$-additive: the above
must also hold for countably many sets (in which case the sum of the measures
becomes a series).

The most important measure, and the one to be used in this article, is the
Lebesgue measure which we’ll denote with the letter $\lambda$. The Lebesgue
measure of an interval of real numbers is its length, e.g.,
$\lambda([1,5])=5-1=4$. And the open interval $(1,5)$ has the same measure.
This implies that the endpoints “don’t count:” finite sets like $\\{1,5\\}$
are null sets, their measure is zero. Generally, a set is a null set if, for
every $\varepsilon>0$, one can find countably many intervals such that the set
is a subset of the union of these intervals and the sum of the measures of the
intervals is at most $\varepsilon$. An important example of a null set is the
set $\mathbb{Q}$ of rational numbers. More generally:

######  Lemma 2.1.

Every countable set is a null set.

###### Proof.

Let $A$ be a countable set and let $(a_{n})_{n\in\mathbb{N}}$ be an
enumeration of $A$.333The enumeration doesn’t have to be one-to-one, so the
proof also applies to finite sets. Also note that in this article
$0\in\mathbb{N}$ and we write $\mathbb{N}^{+}$ for the set of positive
integers. Furthermore, let $\varepsilon$ be an arbitrary positive number. For
each $n\in\mathbb{N}$, let $I_{n}$ be an interval of length
$\varepsilon/2^{n+1}$ which includes $a_{n}$. $A$ is then covered by the
intervals $I_{n}$ and we have:

$\sum_{n=0}^{\infty}\lambda(I_{n})=\sum_{n=0}^{\infty}\frac{\varepsilon}{2^{n+1}}=\varepsilon\qed$

######  Lemma 2.2.

The countable union of null sets is again a null set.

###### Proof.

The main idea is that, for a given $\varepsilon$, we’ll cover the first null
set with intervals which have a total measure of at most $\varepsilon/2$, the
second one with intervals with a total measure of at most $\varepsilon/4$, and
so on, i.e., we’ll again use the geometric series. The details are left to the
reader. ∎

We also note in passing that a subset of a null set is a null set and that
more generally the additivity of $\lambda$ implies $\lambda(A)\leq\lambda(B)$
for $A\subseteq B$ because $B$ is the disjoint union of $A$ und $B\setminus A$
and thus $\lambda(B)=\lambda(A)+\lambda(B\setminus A)$.

If $A$ is a set with a positive measure and $P$ is a property the elements of
$A$ can have or not have, then we say that almost all elements of $A$ have
property $P$ if the set of numbers not having property $P$ is a null set.

##  3 Normal numbers

For $\alpha\in\mathbb{R}$ let

$\alpha=\sum_{j=-d}^{\infty}\alpha_{j,r}r^{-j}$

be its representation in base (or radix) $r\in\mathbb{N}\setminus\\{0,1\\}$.
This means that all $\alpha_{j,r}$ are digits in this base, i.e., elements of
the set $\Sigma_{r}=\\{0,\dots,r-1\\}$. In order to ensure uniqueness, we also
require:

  1. (i)

$\alpha_{-d,r}\neq 0$ for $\alpha\neq 0$ and $d=-1$ for $\alpha=0$

  2. (ii)

There is no $k$ with $\alpha_{j,r}=r-1$ for all $j\geq k$.

We sometimes write $\alpha_{j}$ instead of $\alpha_{j,r}$ if $r$ is implicit.

For example, for $\alpha=\pi$ and $r=10$ we have $d=0$, $\alpha_{0}=3$,
$\alpha_{1}=1$, $\alpha_{2}=4$, and so on. For $\beta=1/2$ and $r=2$ we have
$d=-1$, $\beta_{1}=1$ and $\beta_{j}=0$ for $j>1$. Note that the alternative
representation $d=-2$ and $\beta_{j}=1$ for all $j$ is forbidden by the second
requirement above.444We will soon see that this decision is irrelevant in the
context of normality.

To make up for the fact that we don’t have enough digits for bases greater
than 10, we will sometimes write $[a]$ for the $a$-th digit. For example,
$\pi$ in base $r=100$ will start like this:

$\pi=[3].[14][15][92]\dots$

This means that $\alpha_{2}$ is $[15]$, the 15th digit in base 100.555The 15th
digit is usually written as F in hexadecimal notation.

We now define for $b\in\\{0,\dots,r\\}$:

$p_{b,n,\alpha,r}=|\\{j\in\mathbb{N}^{+}:j\leq n\land\alpha_{j,r}=b\\}|$ (1)

This counts how often the digit $b$ occurs among the first $n$ digits after
the radix point.666We will usually only be concerned with digits after the
radix point and will from now on in most cases omit the phrase “after the
radix point.” Again, we might omit the $r$ (and even the $\alpha$).

As $\alpha=\pi$ starts like this

$3.14159265\bm{3}58979\bm{3}2\bm{3}846264\bm{3}\bm{3}8\bm{3}279502884197169\bm{3}99\bm{3}7510\dots$

in base $r=10$, we have $p_{3,50,\pi}=8$.

A number $\alpha\in\mathbb{R}$ is called simply normal in base $r$, if

$\lim_{n\to\infty}p_{b,n,\alpha,r}/n=1/r$ (2)

holds for all digits $b\in\Sigma_{r}$, i.e., if each digit occurs with the
same relative frequency “in the long run.”

An example that demonstrates how this property depends on the base is the
number $1/3$ which is simply normal in base 2, but obviously not simply normal
in the usual base 10.

$\alpha$ is called normal in base $r$ if $r^{m}\alpha$ is simply normal in
base $r^{n}$ for all $n\in\mathbb{N}^{+}$ and all $m\in\mathbb{N}$. So, the
number $\alpha=1/3$ from above is simply normal in base $r=2$, but it is not
normal in this base because it is not simply normal in base 4 (using $m=0$ and
$n=2$): in base 4 we have $\alpha_{j}=1$ for all $j\in\mathbb{N}^{+}$ and thus
the limit in (2) is 1 for the digit $b=1$ and 0 for the other three digits,
but never $1/4$.

######  Lemma 3.1.

If $r^{m}\alpha$ is simply normal in base $r^{n}$ for all $n\in\mathbb{N}^{+}$
and all $m<n$, then $\alpha$ is normal in base $r$.

###### Proof.

Let $n$ be fixed and use $r^{n}$ as the base. For $m\geq n$ we find non-
negative integers $c$ and $k$ with $c<n$ and $m=kn+c$, and we have:

$r^{m}\alpha=r^{c}\cdot(r^{n})^{k}\alpha$

The digits after the radix point of $(r^{n})^{k}\alpha$ are the same as those
of $\alpha$ beginning at the $k$-th digit. And thus the digits after the radix
point of $r^{m}\alpha$ are the same as those of $r^{c}\alpha$ beginning at the
$k$-the digit.

The reason for this is that the digits in base $r^{n}$ are obtained by
combining groups of $n$ digits in base $r$. Multiplication with $r^{c}$ thus
has the same effect on both sequences of digits. ∎

The following example demonstrates the “shift effect” described in the proof
above (for a number $\alpha$ which is obviously not simply normal in base 10
or 1000):777The horizontal line (vinculum) marks an infinitely repeating digit
sequence. The subscript denotes the base.

$\displaystyle r$ $\displaystyle=10$ $\displaystyle n$ $\displaystyle=3$
$\displaystyle m$ $\displaystyle=7=2\cdot n+1$ $\displaystyle\alpha$
$\displaystyle=0.[123]\overline{[345][42]}_{1000}=0.123\overline{345042}_{10}$
$\displaystyle r^{m}\alpha$ $\displaystyle=10^{7}\cdot
0.123\overline{345042}_{10}=1233450.\overline{423450}_{10}$
$\displaystyle=[1][233][450].\overline{[423][450]}_{1000}$ $\displaystyle
r^{c}\alpha$ $\displaystyle=10\cdot
0.123\overline{345042}_{10}=1.23\overline{345042}_{10}=[1].[233]\overline{[450][423]}_{1000}$

We will call a finite sequence $w=b_{1}b_{2}b_{3}\dots b_{n}$ of digits in
base $r$ an $r$-word (or simply a word) and write $|w|$ for its length $n$.
For the word $\alpha_{n_{1},r}\dots\alpha_{n_{2},r}$ consisting of the digits
of $\alpha$ beginning at position $n_{1}$ and ending at position $n_{2}$,
we’ll write $\alpha_{[n_{1},n_{2}],r}$. For an arbitrary word $w$ we define:

$p_{w,n,\alpha,r}=|\\{j\in\mathbb{N}^{+}:j+|w|-1\leq
n\land\alpha_{[j,j+|w|-1],r}=w\\}|$

This number counts how often the block $w$ of digits appears as a substring of
$\alpha_{[1,n],r}$, i.e., of the first $n$ digits of $\alpha$. Note that this
definition agrees with (1) for words consisting of just one digit.

As an example, consider $\alpha=0.11010111011$. We have $p_{101,11,\alpha}=3$
which means that the word 101 occurs three times among the first 11 digits.
Note that it doesn’t matter that the first two occurrences overlap.

######  Lemma 3.2.

If $\alpha$ is normal in base $r$, then we have for all $r$-words $w$:

$\lim_{n\to\infty}p_{w,n,\alpha,r}/n=1/r^{|w|}$

So, if $\alpha$ is normal in some base, then each finite sequence of digits of
this base, no matter how long, appears infinitely often in the representation
of $\alpha$ and with the same frequency “in the long run” as all other
sequences of the same length. Which is pretty fascinating if you think about
it. Imagine the text of your favorite book stored in a computer file and
viewed as a sequence of ones and zeros. If $\alpha$ is normal in base 2, then
your book will appear infinitely often in the binary representation of
$\alpha$, as will any other book—and your favorite songs as well!

One should think that such “magic” numbers are pretty rare or don’t exist at
all. But the whole purpose of this article is to prove that they are “normal”
in the sense that numbers which don’t have this strange property are extremely
scarce. On the other hand, we don’t know many normal numbers. The ones we do
know about were “breeded” for this purpose while the numbers we deal with on a
daily basis are either obviously not normal, like the rational numbers, or it
is unknown whether they are normal. It is for example an open question whether
$\sqrt{2}$, $\pi$, or $\mathrm{e}$ are normal.

By the way, the property described in lemma 3.2 is sometimes used to define
normality. And it is in fact equivalent to our definition. However, proving
the equivalence requires a lot of technical effort which we’ll forego. See [8]
if you’re interested.

Proving lemma 3.2 is not that hard, though. But instead of a formal proof
(which would probably be confusing because of the notation), we’ll go through
an example which is hopefully illuminating enough to illustrate the general
idea.888If you really need a formal proof, grab any normal number. A very
elegant proof, much better than anything I could write, can be found somewhere
in its digit sequence… Consider the base $r=2$, a number $\alpha$ normal in
this base, and the word $w=11$ consisting of two digits. Furthermore, let
$\varepsilon$ be some positive real number. Because $\alpha$ is normal in base
2, it is simply normal in base 4. That means we can find a number $n_{1}$ such
that for $n\geq n_{1}$ approximately $n/4$ of the first $n$ digits in base 4
are the digit 3. Approximately here is supposed to mean that the actual number
deviates from $n/4$ by no more than $\varepsilon n$. But that implies that
among the first $2n$ digits in the base 2 representation of $\alpha$ we will
have $n/4\pm\varepsilon n$ occurrences of the word 11 (which corresponds to
the digit 3 in base 4).

Base 4: 0 2 1 0 0 3 2 3 1 3 3 2 0 1 0 0 1 2 2 3 0 1 3 2 0 Base 2:
00100100001110110111111000010000011010110001111000

And because $\alpha$ is normal in base 2, $2\alpha$ is also simply normal in
base 4. Which entails that we can find a number $n_{2}$ which has the same
property for $2\alpha$ that $n_{1}$ has for $\alpha$. And we can certainly
arrange for $n_{2}$ to be at least as big as $n_{1}$. So, for $n\geq n_{2}$
we’ll again find $n/4\pm\varepsilon n$ occurrences of the word 11, this time
among the first $2n$ digits of the base 2 representation of $2\alpha$. But
that’s just the base 2 representation of $\alpha$ shifted by one digit and so
these are new occurrences we haven’t counted yet.

Base 4: 1 0 2 0 1 3 1 2 3 3 3 0 0 2 0 0 3 1 1 2 0 3 3 0 1 Base 2:
01001000011101101111110000100000110101100011110001

Combined with the ones we already had we now have $n/2\pm\varepsilon 2n$
places where 11 is a substring. Apart from the possible deviation by
$\varepsilon 2n$ that’s one quarter of $2n$ base 2 digits and that’s what we
needed to show.

The final definition is the following: $\alpha$ is called (absolutely) normal
if it is normal in any integer base greater than 1.

##  4 The main lemma

The proof that almost all numbers are normal relies on a technical lemma which
generalizes a computation from [5, p. 420 f]:

######  Lemma 4.1.

If $r$ is an integer greater than 1, then there’s a positive constant $D$ such
that the following inequality holds for all $n\in\mathbb{N}^{+}$:

$\sum_{p=0}^{n}\binom{n}{p}\frac{(r-1)^{n-p}}{r^{n}}\biggl{(}\frac{p}{n}-\frac{1}{r}\biggr{)}^{\\!4}\leq\frac{D}{n^{2}}$

###### Proof.

We fix positive integers $s$ and $n$ and define some functions recursively:

$\displaystyle f_{0}(x,y)$
$\displaystyle=\sum_{p=0}^{n}\binom{n}{p}x^{sp}y^{n-p}$ $\displaystyle
f_{k+1}(x,y)$ $\displaystyle=x\cdot\frac{\partial}{\partial
x}f_{k}(x,y)-y\cdot\frac{\partial}{\partial
y}f_{k}(x,y)\quad\quad\quad(k\in\mathbb{N})$ By working with individual
summands, it is easy to check that $\displaystyle f_{k}(x,y)$
$\displaystyle=\sum_{p=0}^{n}\binom{n}{p}\bigl{(}(s+1)p-n\bigr{)}^{k}x^{sp}y^{n-p}$
(3)

holds for all $k$.

By the binomial theorem, we know that $f_{0}$ can also be written like this:

$f_{0}(x,y)=(x^{s}+y)^{n}$

It is a tedious—but completely elementary—exercise to compute $f_{4}$ based on
this representation.999As we don’t live in Hausdorff’s times anymore, we can
use a computer algebra system. We get:

$\displaystyle f_{4}(x,y)$ $\displaystyle=n(x^{s}+y)^{n-4}Q$
$\displaystyle\text{with }Q$
$\displaystyle=6n^{2}(s+1)^{2}x^{s}y(sx^{s}-y)^{2}+n^{3}(sx^{s}-y)^{4}+{}$
$\displaystyle\mathrel{\phantom{=}}(s+1)^{4}x^{s}y(x^{2s}-4x^{s}y+y^{2})+{}$
$\displaystyle\mathrel{\phantom{=}}n(s+1)^{3}x^{s}y(7(s+1)x^{s}y-4sx^{2s}-4y^{2})$

We now set $r=s+1$, $x=1/\sqrt[s]{r}$, and $y=(r-1)/r$. This implies
$sx^{s}-y=0$ and two of the four summands in $Q$ vanish. The remaining terms
simplify to this:

$\displaystyle f_{4}(x,y)$
$\displaystyle=3(r-1)^{2}n^{2}+(r^{3}-7r^{2}+12r-6)n$ That’s a second degree
polynomial in $n$ and we thus know that $\displaystyle f_{4}(x,y)$
$\displaystyle\leq Cn^{2}$

for some constant $C$ independent of $n$.

If we now replace $f_{4}(x,y)$ with the term from (3), we get:

$\sum_{p=0}^{n}\binom{n}{p}\frac{(r-1)^{n-p}}{r^{n}}(rp-n)^{4}\leq Cn^{2}$

Dividing by $(rn)^{4}$ yields the inequality we’re after with $D=C/r^{4}$. ∎

##  5 Almost all numbers are normal.

For the rest of this text, we will concentrate on numbers in the interval
$[0,1)$. We fix some base $r\geq 2$. If we look at a specific sequence of $n$
digits, then the set of numbers starting with this sequence is an interval
with Lebesgue measure $1/r^{n}$. For example, in base 10, the set of numbers
starting with the sequence $141$ is the interval $[0.141,0.142)$ which
includes numbers like $\pi-3$. Its measure is $1/1000$.

We now also fix a specific digit $b$ of $\Sigma_{r}$. We want to know the
measure of the set of numbers that have exactly $p$ occurrences of this digit
among their first $n$ digits. That’s also easy to compute: There are
$\binom{n}{p}$ ways to pick $p$ of the available $n$ positions. For the
remaining $n-p$ positions we can pick any of the other $r-1$ digits and there
are $(r-1)^{n-p}$ ways to do that. And each of the sequences thus created
results in an interval of length $1/r^{n}$ disjoint from all other intervals
of the same type. The measure therefore is:

$\binom{n}{p}\cdot(r-1)^{n-p}\cdot\frac{1}{r^{n}}$ (4)

For a positive real number $\varepsilon$ we now look at the set of all numbers
where the relative frequency of $b$’s among the first $n$ digits deviates from
the “expected” value $1/r$ by at least $\varepsilon$:

$\displaystyle M_{b}(n,\varepsilon)$
$\displaystyle=\\{\alpha\in[0,1):|p_{b,n,\alpha,r}/n-1/r|\geq\varepsilon\\}$
With (4), we can compute the measure of this set:
$\displaystyle\lambda(M_{b}(n,\varepsilon))$
$\displaystyle=\sum_{\begin{subarray}{c}p=0\\\
|p/n-1/r|\geq\varepsilon\end{subarray}}^{n}\binom{n}{p}\frac{(r-1)^{n-p}}{r^{n}}$
Using the constant $D$ from lemma 4.1 we get
$\displaystyle\varepsilon^{4}\cdot\lambda(M_{b}(n,\varepsilon))$
$\displaystyle=\sum_{\begin{subarray}{c}p=0\\\
|p/n-1/r|\geq\varepsilon\end{subarray}}^{n}\binom{n}{p}\frac{(r-1)^{n-p}}{r^{n}}\cdot\varepsilon^{4}$
$\displaystyle\leq\sum_{\begin{subarray}{c}p=0\\\
|p/n-1/r|\geq\varepsilon\end{subarray}}^{n}\binom{n}{p}\frac{(r-1)^{n-p}}{r^{n}}\biggl{(}\frac{p}{n}-\frac{1}{r}\biggr{)}^{4}\leq\frac{D}{n^{2}}$
and thus: $\displaystyle\lambda(M_{b}(n,\varepsilon))$
$\displaystyle\leq\frac{D}{\varepsilon^{4}}\cdot\frac{1}{n^{2}}$ (5)

Let $M_{b}(\varepsilon)$ be the set of numbers $\alpha\in[0,1)$ where the
relative frequency $p_{b,n,\alpha,r}/n$ deviates from $1/r$ by at least
$\varepsilon$ for infinitely many $n$. In other words, $\alpha\in
M_{b}(\varepsilon)$ iff for each $m\in\mathbb{N}$ there’s an $n\geq m$ such
that $\alpha\in M_{b}(n,\varepsilon)$:

$\displaystyle M_{b}(\varepsilon)$
$\displaystyle=\bigcap_{m=1}^{\infty}S_{b}(m,\varepsilon)$ $\displaystyle
S_{b}(m,\varepsilon)$
$\displaystyle=\bigcup_{n=m}^{\infty}M_{b}(n,\varepsilon)$

By (5), we have

$\lambda(S_{b}(m,\varepsilon))\leq\sum_{n=m}^{\infty}\lambda(M_{b}(n,\varepsilon))\leq\frac{D}{\varepsilon^{4}}\sum_{n=m}^{\infty}\frac{1}{n^{2}}$

and because the series on the right converges, the measure of
$S_{b}(m,\varepsilon)$ will become arbitrarily small if $m$ is just big
enough. $M_{b}(\varepsilon)$ must therefore be a null set as it is contained
in all $S_{b}(m,\varepsilon)$.

Finally, let $M_{b}$ be the set of all numbers in $[0,1)$ that are not simply
normal in base $r$ because condition (2) is violated by at least the digit
$b$. By the definition of a limit, $M_{b}$ will look like this

$M_{b}=\bigcup_{k=1}^{\infty}M_{b}(1/k)$

and as a countable union of null sets it is itself a null set by lemma 2.2.
The set of the elements of $[0,1)$ which are not simply normal in base $r$ is
then also a null set as it is the union of the $r$ sets $M_{0}$ to $M_{r-1}$.
We just proved:

######  Theorem 5.1.

If $r\geq 2$ is an arbitrary base, then almost all numbers are simply normal
in this base.101010We can drop the restriction to the interval $[0,1)$ if we
want. The proof obviously works just as well for any interval $[m,m+1)$ where
$m$ is an integer and $\mathbb{R}$ is the countable union of such intervals.

We are not quite done yet, but the rest is fairly easy. Let’s again fix a base
$r$. If we multiply each element of $[0,1)$ with a factor $r^{m}$ for some
$m>0$, then the set of products is “spread” over the following intervals:

$r^{m}\cdot[0,1)=[0,1)\cup[1,2)\cup\dots\cup[r^{m}-1,r^{m})$

But as the set of numbers not simply normal in base $r$ in each of these
intervals is a null set (we just proved that), their union is also a null set,
again by lemma 2.2.

Another application of lemma 2.2 yields that the set of numbers
$\alpha\in[0,1)$ such that $r^{m}\alpha$ is not simply normal in base $r$ for
at least one $m\in\mathbb{N}$ is a null set. But the same argument also works
for the bases $r^{2}$, $r^{3}$, and so on. Invoking lemma 2.2 a third time we
get:

######  Corollary 5.2.

If $r\geq 2$ is an arbitrary base, then almost all numbers are normal in this
base.

Finally, you guessed it, we use lemma 2.2 for the last time, utilizing that
there are only countably many bases:

######  Corollary 5.3.

Almost all numbers are absolutely normal.

## References

  * [1] Borel, Émile.  Les probabilités dénombrables et leurs applications arithmétiques.  Suppl. Rend. Circ. mat. Palermo 27, 247–271 (1909)
  * [2] Doob, Joseph Leo.  The development of rigor in mathematical probability.  Am. Math. Monthly 103(4), 586–595 (1996)
  * [3] Filip, Ferdinánd and Šustek, Jan.  An elementary proof that almost all real numbers are normal.  Acta Univ. Sapientiae, Math. 2(1), 99–110 (2010)
  * [4] Hardy, Godfrey Harold and Wright, Edward Maitland.  An introduction to the theory of numbers. Clarendon Press, Oxford (1938)
  * [5] Hausdorff, Felix.  Grundzüge der Mengenlehre. Verlag von Veit & Comp., Leipzig (1914)
  * [6] Kac, Mark.  Statistical independence in probability, analysis and number theory. Carus Math. Monogr., no. 12, Wiley, New York (1959)
  * [7] Nillsen, Rodney.  Normal numbers without measure theory.  Am. Math. Monthly 107(7), 639–644 (2000)
  * [8] Niven, Ivan.  Irrational numbers. Carus Math. Monogr., no. 11, Wiley, New York (1956)


# Beam pointing stabilization of an acousto-optic modulator with thermal
control

Xiao Zhang  1 Yang Chen  1 Jianxiong Fang  1 Tishuo Wang  1 Jiaming Li  1,3
and Le Luo1,2,4  1School of Physics and Astronomy, Sun Yat-Sen University,
Zhuhai, Guangdong, China 519082
2Department of Physics, Indiana University-Purdue University Indianapolis
(IUPUI), Indianapolis, Indiana 46202, USA
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

Diffraction beams generated by an acousto-optic modulator (AOM) are widely
used in various optical experiments, some of which require high angular
stability with the temporal modulation of optical power. Usually, it is
difficult to realize both angular stability and high-power modulation in a
passive setup without a servo system of radio-frequency compensation. Here, we
present a method to suppress the angular drift and pointing noise only with
the thermal management of the AOM crystal. We analyze the dependence of the
angular drift on the refractive index variation, and find that the angular
drift is very sensitivity to the temperature gradient which could induce the
refractive index gradient inside the AOM crystal. It reminds us such angular
drift could be significantly suppressed by carefully overlapping the zero
temperature gradient area with the position of the acousto-optic interaction
zone. We implement a water-cooling setup, and find that the angular drift of
an AOM is reduced over 100 times during the thermal transient, and the angular
noise is also suppressed to 1/3 of the non-cooled case. It should be
emphasized that this thermal control method is a general to suppress the beam
drift in both the diffraction and the perpendicular-to-diffraction directions.
The refractive index thermal coefficient of tellurium dioxide crystal at 1064
nm determined by this angular drift-temperature model is 16$\times$10-6 K-1
consistent with previous studies. This thermal control technique provides
potential applications for optical trapping and remote sensoring that demand
for intensity ramps.

††journal: oe

##  1 Introduction

Acousto-optic modulators are widely used to precisely control light intensity
and frequency, in which a traveling acoustic wave, generated by a
piezoelectric transducer, creates a modulated refractive index and results in
the Bragg diffraction of the incoming light [1]. In ultracold gas experiments,
AOMs are often applied to tune the light intensity of a high power laser to
generate an optical dipole trap (ODT) with time-dependent trap depth for
evaporative cooling or parametric cooling [2, 3]. These applications demand
for minimizing the pointing noise of the diffracted beam so that the wide-band
pointing noise will not parametrically heat ultracold atoms [3, 4]. Meanwhile,
slow angular drifting should also be avoided, which causes the instability of
optical trapping, especially for a cross-beam ODT [5]. Previously, a large
thermal effect on the diffracted beam, induced by radio-frequency (RF)
driving, has been observed in many ODT experiments [6, 7, 8, 9]. It has also
been verified that the angular drift problem is more serious for the high-
power applications, companying with a larger AOM crystal as well as higher RF
driving power [6].

There have been several methods to reduce the angular drift, such as keeping
the total RF power constant with multiple-frequency components during thermal
transient [7, 4], passing AOM twice to cancel out the angular drift [6], RF
servo schemes [9, 10] et. al. However, there are still some drawbacks existing
in these methods, such as constant RF driving at high power causing the beam
shape into elliptical [11], double-passing setup decreasing 25$\%$ of the
output power. Especially, most of the above methods can not compensate the
angular drift along the perpendicular axis of the diffraction direction, for
example Piggott detected a 0.1 mrad angular drift of CrystalTech 3080-197 AOM
in the vertical direction [8].

In this paper, we theoretically modeled the temperature distribution in the
undiffractive direction and measured the refractive index thermal coefficient
of tellurium dioxide crystal. We then proposed and demonstrated a method to
reduce the angular drift by water-cooling. The angular drift of an AOM has
been reduced over 100 times during the thermal transient comparing with the
non-cooling case, and the angular noise is also suppressed to 1/3 accordingly.
Furthermore, we also find this thermal control method is general to the
angular drift suppression in the diffraction direction.

##  2 Experimental setup and theoretical model

Figure 1: Measurement of the angular drift of the diffracted beam of an AOM.
(a) Experimental setup. Light source is a 1064 nm CW fiber laser (IPG,
YLR-100-1064-LP) with an output power of 12 W. The Gaussian Beam waist is 1.1
mm. Beam sampler (Thorlabs, BSF10-B) picks out 2% of the laser power. A CMOS
beam profiler (Dataray, S-WCD-LCM4C) with a 2048(H)$\times$2048(V) pixels and
5.5 $\times$ 5.5 $\mu$m pixel size is put at $D$=9.5 meter away from the AOM.
(b) AOM structure. T${}_{\text{up}}$ and T${}_{\text{down}}$ are the
temperature of the top and bottom water-cooling plates respectively. They are
controlled by two chilled systems independently.

The schematic of the angular drift $\Delta\theta$ measurement is shown in Fig.
1a. The separation angle $2\theta_{B}$ between the zeroth and first order
diffracted beams is determined by $\lambda f/V$, where $\lambda$ is the
wavelength of the incoming light, $f$ is the radio-frequency (RF) of the
acoustic wave, and $V$ is the acoustic-velocity of the medium. The first order
diffraction beam is picked and measured by a beam profiler at the far end. The
exposure time of the beam profiler is limited to sub-millisecond with the
sampling rate of 10 Hz. We fit the beam profile in the horizontal and vertical
direction with a 2D Gaussian profile, and use the center of the Gaussian
fitting as the beam position $(X,Y)$. Therefore, it is easy to get
$\Delta\theta_{x(y)}=\Delta$X(Y)$/D$, where $\Delta$X(Y) is the position shift
related to the initial position at the far end. The AOM (IntraAction,
ATM-804DA6B, tellurium dioxide) has an aperture size of 3 mm. A cartoon
structure of the AOM is illustrated in Fig. 1b. AOM is driven by a 80 MHz RF
wave along the horizontal direction. The highest diffraction efficiency of the
first order is about 92% at 4.5 Watt RF power. The Bragg diffraction angle is
around 10 mrad at 1064 nm wavelength. T${}_{\text{up}}$ and
T${}_{\text{down}}$ are measured by two thermal couples from 15 to 25 ℃ with
0.1 ℃ precision.

We measured $\Delta\theta$ of the AOM caused by thermal transients, as shown
in Fig. 2. For a non-cooled AOM, the beam position in vertical direction is
slowly drifted to one side in about 10 s. The magnitude of $\Delta\theta$ is
almost linear proportion to the RF driving power. With the increase of the RF
power, T${}_{\text{down}}$ increases slower than T${}_{\text{up}}$ because a
copper mount is used in the down side to mount the AOM.

The temperature space-temporal distribution of the AOM crystal $\Delta T(x,t)$
can be approximated by a 1D heat diffusion equation [12]

$\frac{\partial^{2}T(x,t)}{\partial
x^{2}}+\frac{\dot{q}}{k}=\frac{1}{\alpha}\frac{\partial T(x,t)}{\partial t}$
(1)

The RF driving generates an uniform energy per unit volume $\dot{q}$ inside
the AOM. $\alpha$ is the thermal diffusivity, and $k$ is the thermal
conductivity. When $t\rightarrow\infty$, the system reaches to a steady state.
The temperature distribution along the $x$ direction is [13]

$T(x,\infty)=\frac{\dot{q}L^{2}}{2k}(1-\frac{x^{2}}{L^{2}})+\frac{T_{\text{up}}-T_{\text{down}}}{2}\frac{x}{L}+\frac{T_{\text{up}}+T_{\text{down}}}{2}$
(2)

where $L$ is the height of the crystal.

Figure 2: AOM angular drift due to thermal transient. (a) $\Delta\theta$ in
the vertical direction are measured with different input locations of the
input light beam. $\Delta$X is the off center distance of the laser injection
position at AOM in the perpendicular-to-diffraction direction. (b)
$\Delta\theta$ for different RF driving power with the incoming light location
at $\Delta x=0$ (middle point). Blue square, light blue cross, red circle, and
black diamond are measured under 4.5, 3.0, 1.8, and 0.8 W RF power,
respectively.

When the RF driving is turned on, the thermal transient inside the crystal can
be simplified as

$T(x,t)=\theta_{i}e^{-t/\tau}+T(x,\infty)$ (3)

where $\theta_{i}\equiv T(x,0)-T(x,\infty)$ is the initial temperature
difference. The temperature difference $\Delta T(x,t)=T(x,t)-T(x,0)$ affects
the refraction of the crystal. The refractive index as a function of radius
and time can be obtained from $\Delta T(x,t)$ as expressed

$n(x,t)=n_{0}+\frac{dn}{dT}\Delta T(x,t)$ (4)

where $n_{0}$ is the refractive index at the initial temperature. This
equation assumes a decrease in refractive index with increased temperature.

Figure 3: Schematic of the angular drift due to refraction distribution
$n(x,t)$.

Fig. 3 illustrates the angular drift model for the refraction gradient. We
simplify the practical case by assuming 1D temperature gradient, so that
$n_{\text{i}}\cos\theta_{\text{i}}=n_{\text{i+1}}\cos\theta_{\text{i+1}},i\in{1,2,3\ldots}$,
and results in the relation between the input and output angle

$\sin^{2}\theta_{\text{out}}-\sin^{2}\theta_{\text{in}}=n_{i}^{2}-n_{0}^{2}$
(5)

With small angle approximation, the angular shift is given by
$\Delta\theta=\theta_{\text{out}}-\theta_{\text{in}}$

$\Delta\theta=\frac{dn}{dT}\frac{dT}{dx}n_{0}W$ (6)

where $W$ is the length of the crystal.

It is notable that we ignore the beam displacement effect at the AOM in Eq. 6.
The reason is that, when $W\ll D$, the position shift inside the crystal is
much smaller than the position shift after the AOM. We confirm this by
measuring the angular shift at different values of $D$. The result turns out
the beam displacement effect is too small to be considered.

Eq. 2 and Eq. 6 quantitatively explain the drift effect shown in Fig. 2.
Without water cooling, $T_{\text{up}}>T_{\text{down}}$, results in a
temperature gradient along the vertical direction which causes $\Delta\theta$
with RF ramps. In such case, there is no zero temperature gradient area for
suppressing $\Delta\theta$.

##  3 Results

###  3.1 Refractive index thermal coefficient determination

We fix $T_{\text{bottom}}$ to 20 ℃ and change $T_{up}$ from 15 to 25 ℃. Then
the steady angular drift of the beam is measured with a constant RF power. The
result is presented in Fig. 4. The linear dependence of $\Delta\theta$ on the
temperature difference supports the theoretical model.

In this static measurement, $\dot{q}L^{2}/2k$ can be ignored. We simulate
$dT/dx$ with Eq. 2, and get $dT/dx\approx\Delta T/2L$, which means the
temperature gradient is constant. From Eq. 6, we can easily obtain
$dn/dT=16\times 10^{-6}/K$ with $n_{0}=2.2079$ at 20 ℃ [14], $W$=20 mm, $L$=6
mm, which is closed to previous results [15].

Figure 4: The dependence of $\Delta\theta$ on the temperature difference
between the top and bottom surfaces of an AOM. The linear fitting result is
$\Delta\theta/\Delta T=0.059$ mrad/℃.

###  3.2 Reduce the angular drift

Figure 5: AOM angular drift of the thermal transient with a water cooling
setup. (a), the position drift in the vertical direction of AOM, when the
laser position is in different value of the AOM. (b), steady position of (a)
and position noise comparison between water cooling and non-cooled cases. The
standard derivations of water-cooled experiments is 5.5 $\mu$rad, while the
non-cooled one is 17.5 $\mu$rad.

In order to minimize the angular drift of the AOM, we operate the experiments
with water cooling, and $T_{up}=T_{down}=20$ ℃. From Eq. 2, we can get a
parabolic-like temperature distribution with the highest temperature and zero
temperature gradient in the middle point, resulting in null $dn/dx$ in this
spot as shown in Fig. 5, the angular shift is reduced from 0.60 mrad to zero
mrad level from non-cooling case to cooling-case. From Fig. 5a and Fig. 2a, we
can estimate that the magnitude of the refractive index gradient in the middle
point is reduced about 10 times by water-cooling. Fig. 5b shows that the
water-cooled AOM has 3 times lower angular noise than the non-cooled one,
which could be explained by weakening the heating of the crystal [10]. Using
the water cooling scheme, our AOM system could be suitable for cooling in
micro-ODT [16] and lattice trap [17], both of which highly demand pointing
stability and noise reduction.

We test the evaporative cooling with the water-cooled AOM [18]. The result is
shown in Fig. 6, indicating that the cooling setup could solve the angular
drift problem fairly well.

Figure 6: The comparison of the angular drift during an evaporative cooling
process between water-cooled and non-cooled AOMs with exponential decay RF
power for evaporative cooling. The trap depth is lower to 1$\%$ in 20 seconds.

We also use the cooling scheme to suppress the angular drift along the
diffraction direction. Different with the perpendicular direction, there is an
attached piezoelectric transducer at one side of the diffracting direction,
which can not be water-cooled. The solution is that we can control the
temperature of the other side to the same temperature as the transducer side
by thermal control to make an symmetry temperature distribution along the
diffraction direction [6]. The result is shown in Fig. 7.

Figure 7: Measurement results of angular drift in the Bragg diffraction
direction with a water cooling setup in the opposite side of piezoelectric
transducer. $\Delta$Y is the off-the-center distance of the laser injection
position at AOM in the diffraction direction.

##  4 Conclusions

We studied the angular drift and noise of the AOM with RF power ramps. A
thermal-control setup is applied to reduce the angular drift significantly as
well as the pointing-noise. The results agree with the the refractive index-
temperature model of the AOM crystal very well. This technique opens a route
for high-precision beam control for laser cooling and trapping with micro-ODT
and lattice trap.

## Funding

National Natural Science Foundation of China (NSFC) under Grant No.11804406,
Fundamental Research Funds for Sun Yat-sen University 18lgpy78. Science and
Technology Program of Guangzhou 2019-03-01-05-3001-0035. NSFC No.11774436.
SunYat-sen University Discipline Construction Fund, Guangdong ProvinceYouth
Talent Program under Grant No.2017GC010656.

## References

  * [1] I. C. Chang, “I. acoustooptic devices and applications,” IEEE Transactions on Sonics and Ultrasonics 23, 2–21 (1976).
  * [2] R. Grimm, M. Weidemüller, and Y. B. Ovchinnikov, “Optical dipole traps for neutral atoms,” Advances In Atomic, Molecular, and Optical Physics 42, 95 – 170 (2000).
  * [3] J. Li, J. Liu, W. Xu, L. de Melo, and L. Luo, “Parametric cooling of a degenerate fermi gas in an optical trap,” Phys. Rev. A 93, 041401 (2016).
  * [4] L. Luo, “Entropy and superfluid critical parameters of a strongly interacting fermi gas,” Ph.D. thesis, Duke University (2008).
  * [5] C. S. Adams, H. J. Lee, N. Davidson, M. Kasevich, and S. Chu, “Evaporative cooling in a crossed dipole trap,” Phys. Rev. Lett. 74, 3577–3580 (1995).
  * [6] J. Kobayashi, Y. Izumi, M. Kumakura, and Y. Takahashi, “Stable all–optical formation of bose–einstein condensate using pointing–stabilized optical trapping beams,” App. Phys. B 83, 21–25 (2006).
  * [7] B. Fröhlich, T. Lahaye, B. Kaltenhäuser, H. Kübler, S. Müller, T. Koch, M. Fattori, and T. Pfau, “Two-frequency acousto-optic modulator driver to improve the beam pointing stability during intensity ramps,” Review of Scientific Instruments 78, 043101 (2007).
  * [8] A. Piggot, “Angular drift of crystaltech 3080–197 aoms due to thermal transients,” (2010).
  * [9] W. R. Mcgehee, “Transport and disorder-induced localization of unltracold fermi gases,” Ph.D. thesis, University of Illinois at Urbana-Champaign (2015).
  * [10] V. I. Balakshy, V. B. Voloshinov, V. A. Karasev, V. Y. Molchanov, and V. Semenkov, “Compensation of thermal effects in acousto-optic deflector,” Proc. SPIE 2713, 164–171 (1996).
  * [11] “We measured three aoms, such as intraaction model atm-804da6b, isomet model m1306-t80l-4 and model m1135-t80l-3. all of their diffraction beams become elliptical when the high driving power reaches half of the maximum diffraction efficient,” .
  * [12] S. J. Sheldon, L. V. Knight, and J. M. Thorne, “Laser-induced thermal lens effect: a new theoretical model,” Appl. Opt. 21, 1663–1669 (1982).
  * [13] F. P. Incropera, D. P. Dewitt, and T. L. Bergman, _Fundamentals of heat and mass transfer_ (John Wiley Sons, 2007).
  * [14] N. Uchida, “Optical properties of single-crystal paratellurite (te${\mathrm{o}}_{2}$),” Phys. Rev. B 4, 3736–3745 (1971).
  * [15] H. Li, J. Lousteau, W. N. MacPherson, X. Jiang, H. T. Bookey, J. S. Barton, A. Jha, and A. K. Kar, “Thermal sensitivity of tellurite and germanate optical fibers,” Opt. Express 15, 8857–8863 (2007).
  * [16] A. N. Wenz, G. Zürn, S. Murmann, I. Brouzos, T. Lompe, and S. Jochim, “From few to many: Observing the formation of a fermi sea one atom at a time,” Science 342, 457–460 (2013).
  * [17] S. Blatt, A. Mazurenko, M. F. Parsons, C. S. Chiu, F. Huber, and M. Greiner, “Low-noise optical lattices for ultracold ${}^{6}\mathrm{Li}$,” Phys. Rev. A 92, 021402 (2015).
  * [18] K. M. O’Hara, M. E. Gehm, S. R. Granade, and J. E. Thomas, “Scaling laws for evaporative cooling in time-dependent optical traps,” Phys. Rev. A 64, 051403 (2001).


pagerefparensbackrefparens

# Topological Graph Signal Compression

Guillermo Bernárdez∗, Lev Telyatnikov†, Eduard Alarcón∗, Albert Cabellos-
Aparicio∗,
Pere Barlet-Ros∗ and Pietro Liò§
∗Universitat Politècnica de Catalunya, †Sapienza Università di Roma,
§University of Cambridge

###### Abstract

Recently emerged Topological Deep Learning (TDL) methods aim to extend current
Graph Neural Networks (GNN) by naturally processing higher-order interactions,
going beyond the pairwise relations and local neighborhoods defined by graph
representations. In this paper we propose a novel TDL-based method for
compressing signals over graphs, consisting in two main steps: first, disjoint
sets of higher-order structures are inferred based on the original signal –by
clustering $N$ datapoints into $K\ll N$ collections; then, a topological-
inspired message passing gets a compressed representation of the signal within
those multi-element sets. Our results show that our framework improves both
standard GNN and feed-forward architectures in compressing temporal link
signals from two real-world Internet Service Provider Networks’ datasets –from
$30\%$ up to $90\%$ better reconstruction errors across all evaluation
scenarios–, suggesting that it better captures and exploits spatial and
temporal correlations over the whole graph-based network structure.

##  1 Motivation

Graph Neural Networks (GNNs)[27] have demonstrated remarkable performance in
modelling and processing relational data on the graph domain, which naturally
encodes binary interactions. Topological Deep Learning (TDL)[16] methods take
this a step further by working on domains that can feature higher-order
relations. By leveraging (algebraic) topology concepts to encode multi-element
relationships (e.g. simplicial[7], cell[6] and combinatorial complexes[11]),
Topological Neural Networks (TNNs) allows for a more expressive representation
of the complex relational structure at the core of the data. In fact, despite
its recent emergence, TDL is already postulated to become a relevant tool in
many research areas and applications[11], including complex physical
systems[5], signal processing[2, 3], molecular analysis[7, 19] and social
interactions[18].

We argue that the task of data compression can hugely benefit from TDL by
enabling to exploit multi-way correlations between elements beyond pre-defined
local neighborhoods to get the desired lower-dimensional representations. To
the best of our knowledge, current Machine Learning (ML) compression
approaches mainly rely on Information Theory (IT) and are narrowed to Computer
Vision applications[13, 24]. In contrast to that, and inspired by zfp[9] –the
current state-of-the-art lossy compression method for floating-point data,
more details in A.2–, we propose in this paper a novel TDL framework to (a)
first detect higher-order correlated structures over a given data, and (b)
then directly apply TNNs to obtain compressed representations within those
multi-element sets.

This work provides evidence supporting that TDL could have great potential in
compressing relational data. With the long-term objective of outperforming
zfp, our current goal is to assess if the proposed framework naturally
exploits multi-datapoint interactions –between possibly distant elements– in a
way that makes it more suitable for compression than other ML architectures
(even if data comes from the graph domain). To do so, we consider the critical
problem of traffic storage in today’s Internet Service Providers (ISP)
networks[1], and set the target to compress the temporal per-link traffic
evolution –Figure 4– for two real-world datasets extracted from [15] (more
details and motivation of this use case provided in A.1). Once the original
link-based signal is divided into processable temporal windows, we benchmark
our method against a curated set of GNN-based architectures –and a Multi-Layer
Perceptron (MLP)– properly designed for compression as well. Obtained results
clearly suggest that our topological framework defines the best baseline for
lossy neural compression.

##  2 Methodology

This section describes the proposed Topological (Graph) Signal Compression
framework, which is divided into the following three primary modules:

##### 1) Topology Inference Module.

The first stage of the proposed model infers the computational topological
structure –both pairwise and higher-order relationships– from the data
measurements. In general, the framework assumes to have a set $\mathcal{S}$ of
$M$ signals, $\mathcal{S}=\\{S_{i}\\}_{i=1}^{M}$, where $S_{i}$ consists of
$N$ vector-valued measurements $x_{j}$ of a pre-defined dimension $d$, i.e.
$S_{i}=\\{x^{i}_{j}\\}_{j=1}^{N}$, $x^{i}_{j}\in\mathbb{R}^{d}$. Thus, the
pipeline that we describe as follows (see also Figure 1) is independently
applied to every signal $S_{i}$.

Figure 1: Topology Inference Module. For each subsignal $S_{i}\in\mathcal{S}$,
it outputs a topological object
$\mathcal{T}=\left(\mathcal{V},\mathcal{E},\mathcal{W}\right)$ determined by
$K$ disjoint hyperedges.

Similarity Matrix: The initial signal $\\{x_{j}\\}_{j=1}^{N}$ 111For the sake
of simplicity, and as abuse of notation, we will avoid writing the superscript
$i$ when referring to the measurements of a generic signal
$S_{i}\in\mathcal{S}$. is encoded with a MLP into an embedding space
$h_{j}^{0}=\psi_{\theta_{\mathcal{V}}}(x_{j})\in\mathbb{R}^{d^{\prime}}$,
$\forall j\in\\{1,\dots,N\\}$. Next, we compute the pairwise similarity matrix
$M_{S}=(m_{uv})\in\mathbb{R}^{d^{\prime}\times d^{\prime}}$ where
$m_{uv}:=f_{S}\left(h^{0}_{u},h^{0}_{v}\right)$ and
$f_{S}:\mathbb{R}^{d^{\prime}}\times\mathbb{R}^{d^{\prime}}\rightarrow\mathbb{R}$
is a similarity function.

Higher-order Relationships: We use clustering techniques on the similarity
matrix $M_{S}$ to deduce $K$ higher-order structures, over which a Topological
Message Passing pipeline –see next module 2)– performs the compression. In
fact, the idea is to compress the signal within the inferred multi-element
sets and encode compressed representations of the data into the final hidden
states of these hyperedges. Therefore, the number of higher-order structures
$K$ is desired to be considerably lower than the number $N$ of datapoints
($K\ll N$); we design the following clustering scheme for this purpose:

  1. 1.

The number of hyperedges are defined as $K=\lfloor N/p\rceil$, where $p$ is a
hyperparameter that identifies the maximum hyperedge length.

  2. 2.

For every row in the similarity matrix $M_{S}$, we extract the top $p-1$
highest entries and calculate their sum. We then select the row that
corresponds to the highest summation value. This chosen row becomes the basis
for forming a hyperedge as we gather the indices of the $p-1$ selected columns
along with the index of the row itself. Then the gathered indices are removed
from the rows and columns of the similarity matrix $M_{S}$, obtaining a
reduced $\hat{M}_{S}\in\mathbb{R}^{(d^{\prime}-p)\times(d^{\prime}-p)}$.

  3. 3.

Previous step 2 is repeated with subsequent $\hat{M}_{S}$ until $K$ disjoint
hyperedges are obtained.222When $N/p$ is not an even division, at some point
of the process the ranking starts considering the row-wise $p-2$ higher
entries to form $p-1$-length hyperedges, so that at the end a total of
$K=\lfloor N/p\rceil$ hyperedges of lengths $p$ and $p-1$ are obtained; see
A.4.2 for further details.

On the other hand, the choice of the similarity function becomes a crucial
aspect for the compression task. Supported by our early experiments (see
Section 3), our framework makes use of the Signal to Noise Ratio (SNR)
distance metric presented in [25], proposed in the context of deep metric
learning as it jointly preserves the semantic similarity and the correlations
in learned features[25].

Pairwise relationships: Besides higher-order structures, our framework can
optionally leverage graph-based relational interactions, either (i) by
considering the original graph connectivity if it is known, or (ii) by
inferring the edges via the similarity matrix as well –by connecting each
element with a subset of top $k$ row-based entries in $M_{S}$. In our
experiments only intra-hyperedge edges have been considered to keep the
inferred higher-order structures completely disjoint from each other.

##### 2) Compression Module via Topological Message Passing

We implemented two topological Message Passing (MP) compression pipelines,
named SetMP and CombMP. SetMP is a purely set-based architecture that operates
only over hyperedges and nodes; more details in A.3. In this section we
describe CombMP, our most general architecture that leverages the three
different structures (nodes, edges, hyperedges) in a hierarchical way,333Edges
and hyperedges are distinguished because, analogously to recent Combinatorial
Complexes (CCC) models[11], edges can hierarchically communicate with
hyperedges if they are contained in them; in fact, the name CombMP relates to
these general topological constructions (more details in Appendix A.2). and
can be seen as a generalisation of SetMP.

Figure 2: Compression Module workflow for the CombMP architecture; it is
independently applied to each hyperedge $w\in\mathcal{W}$ of the inferred
topological object
$\mathcal{T}=\left(\mathcal{V},\mathcal{E},\mathcal{W}\right)$. The arguments
of the functions visually represent either initial node features or the
corresponding element embedding.

For a given signal $S_{i}$ and its corresponding initial embeddings
$\\{h_{j}^{0}\\}_{j=1}^{N}$, our model operates over a topological object
$\mathcal{T}=\left(\mathcal{V},\mathcal{E},\mathcal{W}\right)$ where
$\mathcal{V}$ denotes the set of elements or nodes, $|\mathcal{V}|=N$;
$\mathcal{E}\in\mathcal{V}\times\mathcal{V}$ represent the set of edges; and
$\mathcal{W}\in(\mathcal{V}\times\dots\times\mathcal{V})$ the set of
hyperedges. The compression pipeline (visualized in Figure 2) can be described
as follows:

Initial embeddings: First, we generate initial embeddings for the three
considered topological structures. For nodes, we use the previously computed
embeddings $\\{h_{v}^{0}\\}_{v=1}^{N}$. For edges and hyperedges, (learnable)
permutation invariant functions are applied over the initial embeddings of the
nodes they contain; respectively,
$h_{e}^{0}=\phi_{\theta_{\mathcal{E}}}\left(\oplus_{v\in e}h_{v}^{0}\right)$
for each $e\in\mathcal{E}$, and
$h_{w}^{0}=\phi_{\theta_{\mathcal{W}}}\left(\oplus_{v\in w}h_{v}^{0}\right)$
for each $w\in\mathcal{W}$. The same dimension $d^{\prime}$ is used for all
initial and intermediate hidden representations.

Edge-Hyperedge Message Passing: We define a hierarchical propagation of
messages between edges and hyperedges. First, neighboring edges communicate to
each other to update their representations; denoting the edge neighbors of an
edge $e\in\mathcal{E}$ by
$\mathcal{N}_{e}^{\mathcal{E}}:=\\{e^{\prime}=(u,v)\in\mathcal{E}|e^{\prime}\neq
e,u\in e\lor v\in e\\}$, its new hidden state becomes
$h_{e}^{1}=\phi_{\theta_{\mathcal{E}\rightarrow\mathcal{E}}}\left(\oplus_{e^{\prime}\in\mathcal{N}_{e}^{\mathcal{E}}}\psi_{\theta_{\mathcal{E}\rightarrow\mathcal{E}}}\left(h_{e}^{0},h^{0}_{e^{\prime}}\right)\right)$.
Next, hyperedges also update their hidden states based on the updated edge
representations according to
$h_{w}^{1}=\phi_{\theta_{\mathcal{E}\rightarrow\mathcal{W}}}\left(\oplus_{e\in\mathcal{E},e\subset
w}\psi_{\theta_{\mathcal{E}\rightarrow\mathcal{W}}}\left(h_{w}^{0},h^{1}_{e}\right)\right)$,
for each $w\in\mathcal{W}$. Then the idea is to propagate downwards towards
the edges, i.e. from hyperedges to edges,
$h_{e}^{2}=\phi_{\theta_{\mathcal{W}\rightarrow\mathcal{E}}}\left(\oplus_{w\in\mathcal{W},e\subset
w}\psi_{\theta_{\mathcal{W}\rightarrow\mathcal{E}}}\left(h_{e}^{1},h^{1}_{w}\right)\right)$;
and only between edges again. This whole communication process can be iterated
$T$ times.

Edge-to-Node Compression: At this point, we perform a first compression step
over the nodes by leveraging the updated edge hidden representations, the
initial node embeddings, as well as the original node data as a residual
connection. Formally, for each node $v\in\mathcal{V}$ we get a compressed
hidden representation
$h_{v}^{c}=\phi_{\theta_{\mathcal{E}\rightarrow\mathcal{V}}}\left(\oplus_{e\in\mathcal{E},v\in
e}\psi_{\theta_{\mathcal{E}\rightarrow\mathcal{V}}}\left(x_{v},h_{v}^{0},h^{t}_{e}\right)\right)\in\mathbb{R}^{d_{\mathcal{V}}^{c}}$.

Node-to-Hyperedge Compression: Finally, a second and last compression step is
performed over the hypergraph representations, in this case leveraging a
residual connection to the original measurements, the previously computed
compressed representations of nodes, as well as the updated hidden
representations of hyperedges. More in detail, each hyperedge
$w\in\mathcal{W}$ obtains its final compressed hidden representation as
$h_{w}^{c}=\phi_{\theta_{\mathcal{V}\rightarrow\mathcal{W}}}\left(\oplus_{v\in\mathcal{V},v\in
w}\psi_{\theta_{\mathcal{V}\rightarrow\mathcal{W}}}\left(x_{v},h_{v}^{c},h^{t}_{w}\right)\right)\in\mathbb{R}^{d_{\mathcal{W}}^{c}}$.

The final node and hyperedge states,
$\left\\{\\{h_{v}^{c}\\}_{v\in\mathcal{V}},\\{h_{w}^{c}\\}_{w\in\mathcal{W}}\right\\}$,
encode the compressed representation of a signal
$S_{i}=\\{x_{j}\\}_{j=1}^{N}$. Consequently, the compression factor $r_{c}$
can be expressed as: $r_{c}=\frac{N\cdot d_{\mathcal{V}}^{c}+K\cdot
d_{\mathcal{W}}^{c}}{N\cdot d}$ (1)

Figure 3: Decompression Module. It is applied over each hyperedge-dependent
compressed representation set generated by the Compression Module.

##### 3) Decompression Module.

This last module learns to reconstruct the original signal of every node
through its compressed representation and the final hidden state of the
hyperedge it belongs to. More formally, for each $v\in\mathcal{V}$ and its
corresponding hyperedge $v\in w\in\mathcal{W}$, the reconstructed signal
$\hat{x}_{v}$ is obtained as
$\hat{x}_{v}=\phi_{dec}\left(h_{v}^{c},h_{w}^{c}\right)$, where $\phi_{dec}$
is implemented as a MLP in our framework. The whole model is trained end-to-
end to minimize the (mean squared) reconstruction error.

##  3 Evaluation & Discussion

Table 1: Reconstruction Mean Squared Error (MSE) and Mean Absolute Error (MAE)
over the test set of the considered datasets for two different compression
factors. Top: Methods that leverage the original graph-based network
structure. Middle: Methods with no inductive biases. Bottom: Methods that
leverage higher-order structures (ours).

| Abilene | Geant
---|---|---
$r_{c}=1/3$ | $r_{c}=2/3$ | $r_{c}=1/3$ | $r_{c}=2/3$
MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE
GNN | $1.95\cdot 10^{-2}$ | $1.08\cdot 10^{-1}$ | $1.95\cdot 10^{-2}$ | $1.08\cdot 10^{-1}$ | $2.33\cdot 10^{-2}$ | $1.21\cdot 10^{-1}$ | $2.32\cdot 10^{-2}$ | $1.20\cdot 10^{-1}$
MPNN | $7.88\cdot 10^{-4}$ | $1.24\cdot 10^{-2}$ | $7.92\cdot 10^{-4}$ | $1.24\cdot 10^{-2}$ | $8.45\cdot 10^{-3}$ | $4.13\cdot 10^{-2}$ | $1.82\cdot 10^{-3}$ | $2.39\cdot 10^{-2}$
MLP | $1.04\cdot 10^{-3}$ | $1.88\cdot 10^{-2}$ | $9.71\cdot 10^{-4}$ | $1.80\cdot 10^{-2}$ | $3.76\cdot 10^{-3}$ | $3.96\cdot 10^{-2}$ | $3.62\cdot 10^{-3}$ | $3.89\cdot 10^{-2}$
SetMP | $3.22\cdot 10^{-4}$ | $8.75\cdot 10^{-3}$ | $2.03\cdot 10^{-4}$ | $6.80\cdot 10^{-3}$ | $6.93\cdot 10^{-4}$ | $1.52\cdot 10^{-2}$ | $2.90\cdot 10^{-4}$ | $1.05\cdot 10^{-2}$
CombMP | $5.81\cdot 10^{-4}$ | $1.12\cdot 10^{-2}$ | $3.76\cdot 10^{-4}$ | $1.06\cdot 10^{-2}$ | $1.07\cdot 10^{-3}$ | $1.88\cdot 10^{-2}$ | $7.04\cdot 10^{-4}$ | $1.61\cdot 10^{-2}$

Experimental Setup: For the evaluation, we use two public real-world datasets
–Abilene, Geant– from [15]. They are pre-processed to generate subsignals
$S_{i}$ of network link-level traffic measurements in temporal windows of
length $d=10$, to which then a random $60/20/20$ split is performed for
training, validation and test, respectively. In this context, our method is
compared against:

  * •

w we implemented several standard GNN architectures (GCN[14], GAT[21],
GATv2[8], GraphSAGE[12]) to perform signal compression over the network graph
topology; we take the best result among them in each evaluation scenario.

  * •

MPNN: a custom MP-based GNN scheme –over the original network graph structure
as well– whose modules and pipeline are similar to our proposed topological
MPs.

  * •

MLP: a feed-forward auto-encoder architecture with no inductive biases over
subsignals $S_{i}$.

GNN and MPNN baselines implement a decompression module similar to that of our
TDL-based methods. More details about the evaluation and all model
implementations are provided in A.4.

Experimental Results: Table 1 shows the reconstruction error (MSE and MAE)
obtained by our framework and the baselines in both datasets for two
compression factors, $1/3$ and $2/3$. SetMP, our topological edge-less
architecture, clearly performs the best in all scenarios –improving on average
by $75\%$ and $48\%$, respectively, the best MSE and MAE obtained by
baselines–, followed by our most general CombMP method. As for the baselines,
in overall MPNN slightly outperforms MLP, and GNN performs the worst. In
addition, a comparison against state-of-the-art zfp can be found in A.4.5.

Discussion: These results support our hypothesis that taking into account
higher-order interactions could help in designing more expressive ML-based
models for (graph) signal compression tasks, specially due to the fact that
these higher-order structures can go beyond the (graph) local neighborhood and
connect possibly distant datapoints whose signals may be strongly correlated
(e.g. generator and sink nodes in ISP Networks). In that regard, TDL can
provide us with novel methodologies that naturally encompass and exploit those
multi-element relations. Moreover, it is interesting to see how our set-based
architecture outperforms the combinatorial-based one in every scenario,
suggesting that intermediate binary connections might add noise in the process
of distilling compressed representations. Further discussion on future work,
focusing on the current limitations of our method and how to possibly address
them, can be found in A.5.

## Acknowledgements

This publication is part of the Spanish I+D+i project TRAINER-A (ref.
PID2020-118011GB-C21), funded by MCIN/AEI/10.13039/ 501100011033. This work is
also partially funded by the Catalan Institution for Research and Advanced
Studies (ICREA), the Secretariat for Universities and Research of the Ministry
of Business and Knowledge of the Government of Catalonia, and the European
Social Fund.

## References

  * [1] Paul Almasan et al.  “Leveraging Spatial and Temporal Correlations for Network Traffic Compression”  In _arXiv preprint arXiv:2301.08962_ , 2023
  * [2] Sergio Barbarossa and Stefania Sardellitti  “Topological signal processing over simplicial complexes”  In _IEEE Transactions on Signal Processing_ 68 IEEE, 2020, pp. 2992–3007
  * [3] Claudio Battiloro, Paolo Di Lorenzo and Sergio Barbarossa  “Topological Slepians: Maximally Localized Representations of Signals Over Simplicial Complexes”  In _ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2023, pp. 1–5  IEEE
  * [4] Claudio Battiloro et al.  “From Latent Graph to Latent Topology Inference: Differentiable Cell Complex Module”  In _arXiv preprint arXiv:2305.16174_ , 2023
  * [5] Federico Battiston et al.  “The physics of higher-order interactions in complex systems”  In _Nature Physics_ 17.10 Nature Publishing Group UK London, 2021, pp. 1093–1098
  * [6] Cristian Bodnar et al.  “Weisfeiler and lehman go cellular: Cw networks”  In _Advances in Neural Information Processing Systems_ 34, 2021, pp. 2625–2640
  * [7] Cristian Bodnar et al.  “Weisfeiler and lehman go topological: Message passing simplicial networks”  In _International Conference on Machine Learning_ , 2021, pp. 1026–1037  PMLR
  * [8] Shaked Brody, Uri Alon and Eran Yahav  “How attentive are graph attention networks?”  In _ICLR_ , 2022
  * [9] James Diffenderfer et al.  “Error analysis of zfp compression for floating-point data”  In _SIAM Journal on Scientific Computing_ 41.3 SIAM, 2019, pp. A1867–A1898
  * [10] Matthias Fey and Jan Eric Lenssen  “Fast graph representation learning with PyTorch Geometric”  In _arXiv preprint arXiv:1903.02428_ , 2019
  * [11] Mustafa Hajij et al.  “Topological Deep Learning: Going Beyond Graph Data”, 2023  arXiv:2206.00606 [cs.LG]
  * [12] Will Hamilton, Zhitao Ying and Jure Leskovec  “Inductive representation learning on large graphs”  In _Advances in Neural Information Processing Systems_ 30, 2017
  * [13] Marton Havasi  “Advances in Compression using Probabilistic Models”  Apollo - University of Cambridge Repository, 2021  DOI: 10.17863/CAM.79008
  * [14] Thomas N. Kipf and Max Welling  “Semi-Supervised Classification with Graph Convolutional Networks”  In _ICLR_ , 2017
  * [15] Sebastian Orlowski, Roland Wessäly, Michal Pióro and Artur Tomaszewski  “SNDlib 1.0—Survivable network design library”  In _Networks: An International Journal_ 55.3 Wiley Online Library, 2010, pp. 276–286
  * [16] Mathilde Papillon, Sophia Sanborn, Mustafa Hajij and Nina Miolane  “Architectures of topological deep learning: A survey on topological neural networks”  In _arXiv preprint arXiv:2304.10031_ , 2023
  * [17] Arjun Roy et al.  “Inside the social network’s (datacenter) network”  In _Proceedings of the 2015 ACM Conference on Special Interest Group on Data Communication_ , 2015, pp. 123–137
  * [18] Michael T Schaub et al.  “Random walks on simplicial complexes and the normalized Hodge 1-Laplacian”  In _SIAM Review_ 62.2 SIAM, 2020, pp. 353–391
  * [19] Yair Schiff, Vijil Chenthamarakshan, Karthikeyan Natesan Ramamurthy and Payel Das  “Characterizing the latent space of molecular deep generative models with persistent homology metrics”  In _arXiv preprint arXiv:2010.08548_ , 2020
  * [20] Paul Tune, Matthew Roughan, H Haddadi and O Bonaventure  “Internet traffic matrices: A primer”  In _Recent Advances in Networking_ 1 ACM SIGCOMM, 2013, pp. 1–56
  * [21] Petar Veličković et al.  “Graph Attention Networks”  In _ICLR_ , 2018
  * [22] Chenxin Xu et al.  “Groupnet: Multiscale hypergraph neural networks for trajectory prediction with relational reasoning”  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 6498–6507
  * [23] Fengli Xu et al.  “Understanding mobile traffic patterns of large scale cellular towers in urban environment”  In _IEEE/ACM transactions on networking_ 25.2 IEEE, 2016, pp. 1147–1161
  * [24] Yibo Yang, Stephan Mandt and Lucas Theis  “An introduction to neural data compression”  In _Foundations and Trends® in Computer Graphics and Vision_ 15.2 Now Publishers, Inc., 2023, pp. 113–200
  * [25] Tongtong Yuan et al.  “Signal-to-noise ratio: A robust distance metric for deep metric learning”  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2019, pp. 4815–4824
  * [26] Manzil Zaheer et al.  “Deep sets”  In _Advances in neural information processing systems_ 30, 2017
  * [27] Jie Zhou et al.  “Graph neural networks: A review of methods and applications”  In _AI open_ 1 Elsevier, 2020, pp. 57–81

##  Appendix A Appendix

###  A.1 Network Traffic Compression

Computer Network’s traffic has significantly increased in recent years[20],
specially driven by the development of new applications –such as vehicular
networks, Internet of Things, virtual reality, video streaming– and the
advancement of network technology –e.g. the fast improvements in link speed.
In fact, current Internet Service Providers (ISP) Networks can easily produce
hundreds of terabytes of traffic traces per day[17], and this keeps growing.

However, network operators continuously need to store and analyze network
traffic data for various network management purposes, including network
planning, traffic engineering, traffic classification, anomaly detection or
network forensics. With those huge amounts of generated data, the efficient
storage of all this information is then becoming a crucial aspect for
them[23].

Figure 4: The goal is to compress a signal $\mathcal{S}$ over graph-based ISP
Networks representing the temporal evolution of each link utilization.

This is what motivated us to choose ISP Traffic Compression for testing our
proposed method. Not only it provides with complex data naturally represented
in the graph domain, but also represents a relevant use case for the
networking community. In addition to that, for network management tasks there
is a reasonable loss tolerance of the compression, so it would make sense to
consider lossy methods such as zfp[9] or ours. Finally, it is also important
to note that there exists public real-world datasets of ISP backbone
networks[15] that, despite of their limited network sizes, already reflect
complex traffic patterns that may go beyond the provided graph structure (e.g.
with distant elements possibly having strong correlations, such as links that
are adjacent to generator and sink nodes). Thus, we argue that this use case
perfectly serves for our first testing purposes. In particular, in our
experiments we consider Abilene and Geant datasets from [15], which are the
ones with the higher number of traffic traces; more details about them in
Appendix A.4.1.

###  A.2 Related Work

As already mentioned in Section 1, there already exist ML-based models in the
literature that target compression tasks[13, 24], but they are mainly
entangled to Information Theory concepts used in classical compression
algorithms, and applied to Computer Vision domains. Our approach differs from
them in these two basic aspects, and has been inspired by zfp[9], the state-
of-the-art lossy method for floating-point data compression. As detailed in
[9], the first step of zfp consists in dividing floating matrices or tensors
in disjoint blocks of a fixed dimension, which then are independently
processed to extract compressed representations. This has obvious similarities
with our search of higher-order structures, with the difference that in zfp
the divisions are totally determined by the input elements’ order.

Precisely our proposed topology inference module has some resemblances to that
of [22], which is also based on grouping entries of an affinity (i.e.
similarity) matrix computed over a set of element’s neural embeddings.
Nevertheless, due to the constraints imposed by the compression task there
exists relevant differences between that inference procedure and ours: whereas
we design a clustering methodology to get disjoint and uniform-length sets, in
[22] they implement a multi-scale hyperedge forming pipeline where each node
ends up belonging to an arbitrary number of higher-order structures –and not
only among different hyperedge degrees, but also within the same scale.

Regarding our topological-inspired MP architectures, we highlight the paper of
[11] on Combinatorial Complexes (CCC) and the survey [16] on TDL
architectures. Our most general method –CombMP– was conceptualized before the
publication of these works, but we note that it can be formalized in terms of
CCCs’ notation by considering nodes as 0-cells, edges as 1-cells, and
hyperedges as 2-cells; it is due to this fact that we called it CombMP, which
stands for Combinatorial Message Passing. SetMP, on the other hand, belongs to
the hypergraph family of TDL architectures according to the classification of
[16]. More precisely, it can be linked to DeepSets architectures[26], which is
why it received that name.

Finally, a special mention should be given to [1], which also leverages a ML
model to specifically perform traffic compression on ISP networks. However,
authors of this work implement a spatio-temporal GNN that acts as a predictor
of a traditional lossless compression method (in this case, Arithmetic
Coding), which defines a totally different conceptual approach.

###  A.3 Compression via Topological Message Passing (SetMP)

Figure 5: Compression Module workflow for the SetMP architecture; it is
independently applied to each inferred hyperedge $w\in\mathcal{W}$. The
arguments of the functions visually represent either initial node features or
the corresponding element embedding.

In this section we describe the topological MP pipeline of SetMP (Figure 5).
In contrast to CombMP, this architecture disregards binary connections and
consequently operates over a topological object
$\mathcal{T}=\left(\mathcal{V},\mathcal{W}\right)$, where again $\mathcal{V}$
denotes the set of $N$ nodes and
$\mathcal{W}\in(\mathcal{V}\times\dots\times\mathcal{V})$ the set of $K$
inferred hyperedges. We describe the differences in the compression pipeline
in this scenario:

##### Initial embeddings:

Initial embeddings for nodes and hyperedges are generated in the same way:
$\\{h_{v}^{0}\\}_{v=1}^{N}$ for the nodes, and
$h_{w}^{0}=\phi_{\theta_{\mathcal{W}}}\left(\oplus_{v\in w}h_{v}^{0}\right)$
for each $w\in\mathcal{W}$ for the hyperedges (both of them with dimension
$d^{\prime}$).

##### Hyperedge-to-Node Compression:

Without edges as intermediaries, we directly perform the node compression, in
this case based on both the initial node and hyperedge embeddings and the
original signal $S_{i}$. Formally, for each node $v\in\mathcal{V}$ we get a
compressed hidden representation
$h_{v}^{c}=\phi_{\theta_{\mathcal{W}\rightarrow\mathcal{V}}}\left(\oplus_{w\in\mathcal{W},v\in
e}\psi_{\theta_{\mathcal{W}\rightarrow\mathcal{V}}}\left(x_{v},h_{v}^{0},h^{0}_{w}\right)\right)\in\mathbb{R}^{d_{\mathcal{V}}^{c}}$.

##### Node-to-Hyperedge Compression:

The second and last compression step over the hyperedge representations is
exactly the same as in the CombMP –except for the fact that now the considered
hyperedge hidden states have not been updated. Therefore, each hyperedge
$w\in\mathcal{W}$ obtains its final compressed hidden representation as
$h_{w}^{c}=\phi_{\theta_{\mathcal{V}\rightarrow\mathcal{W}}}\left(\oplus_{v\in\mathcal{V},v\in
w}\psi_{\theta_{\mathcal{V}\rightarrow\mathcal{W}}}\left(x_{v},h_{v}^{c},h^{0}_{w}\right)\right)\in\mathbb{R}^{d_{\mathcal{W}}^{c}}$.

Thus, again the final node and hyperedge representations
–$\left\\{\\{h_{v}^{c}\\}_{v\in\mathcal{V}},\\{h_{w}^{c}\\}_{w\in\mathcal{W}}\right\\}$–
encode the compressed representation of a signal
$S_{i}=\\{x_{j}\\}_{j=1}^{N}$, and the same expression of the compression
factor $r_{c}$ applies to SetMP. Finally, we note that the decompression
module does not vary either.

###  A.4 Further Details of Evaluation

In this section we provide more technical details about the datasets and model
implementations. Lastly, we also show the performance of zfp in the considered
evaluation scenarios.

####  A.4.1 Datasets

As stated in Section 3, the traffic traces of both datasets are publicly
available at [15]. After distributing them into links, Abilene dataset
contains link-level traffic utilization measurements over 6 months –in
intervals of 5 minutes– for a topology with 12 nodes and 30 directional links.
Considering a temporal window of length $d=10$, this results in $N=4,809$
subsignal samples after data cleaning. On the other hand, Geant dataset
contains analogous measurements for a period of 4 months and a time interval
of 15 minutes, in this case for a topology with 22 nodes and 72 directional
links. After data cleaning, it has a total of $1,075$ final samples with the
same window size.

The topology structure of both networks is also provided, and we use it in our
GNN-based baselines. The aforementioned $60/20/20$ split is performed over
these resulting link-based subsignals.

####  A.4.2 Implementation of our Proposed Models

Regarding the Topology Inference module, we recall the relevance of the hyper-
parameter $p$ that defines the maximum allowed hyperedge length, and which
also defines the number of those hyperedges by $K=\lfloor N/p\rceil$, being
$N$ the total number of datapoints. In particular, in our implementation we
try to form $K$ $p$-uniform disjoint hyperedges if possible, but otherwise
build a combination of $p-1$ and $p$-uniform disjoint hyperedges –so that
every datapoint is contained in one of them. Moreover, we always consider
$p>4$ so that $K\ll N$ holds.

As shown in Equation 1, this parameter $p$ together with the dimensions of the
final node and hyperedge compressed representations, $d_{\mathcal{V}}^{c}$ and
$d_{\mathcal{W}}^{c}$, define the compression factor of our method. After some
hyperparameter tuning, in our experiments we used $p=8$,
$d_{\mathcal{V}}^{c}=2$ and $d_{\mathcal{W}}^{c}=10$ for achieving
$r_{c}=1/3$, and $p=6$, $d_{\mathcal{V}}^{c}=5$ and $d_{\mathcal{W}}^{c}=10$
for $r_{c}=2/3$; respectively, this resulted in $4$ and $5$ inferred
hyperedges for the Abilene dataset, and $9$ and $12$ for Geant.

These parameters apply to both SetMP and CombMP architectures, and in both
scenarios we also consider node, edge and hyperedge hidden representations of
dimension $d^{\prime}=20$, double the the dimension $d$ of node signals. All
message and update functions, $\psi_{\theta_{\cdot}}$ and
$\phi_{\theta_{\cdot}}$, are implemented as MLPs, and permutation-invariant
aggregator functions $\oplus$ consist of the concatenation three element-wise
operations: mean, max and min. In the case of CombMP, only 1 iteration of
topological message passing (Figure 2.2) is performed. Finally, we note that
the Decompression Module has the same MLP structure in both compression
pipelines.

We will publicly release the code in the future extended version of this work.

####  A.4.3 Baseline Implementations

Analogously to what we do with our proposed architectures, we have fine-tuned
the involved hyperparameters of all implemented baselines to perform the
compression task. Moreover, they follow a similar compression scheme than our
proposed TDL methods. We provide more details below:

##### Graph-based

These methods leverage the original graph-like network structure present in
Abilene and Geant datasets. In this case, since the signal is over the edges,
we compute the (dual) line graph of the network, so that edges become nodes
and are connected between them if they share origin/destination. The idea is
then to perform several iterations of message passing over this line graph to
get a compressed representation of the original signal in the hidden state of
these link-based nodes. The difference between our two graph-based baselines
precisely relies on the nature of this message passing:

  * •

GNN: Under this name we gather the results of implementing several standard
GNN architectures (GCN[14], GAT[21], GATv2[8], GraphSAGE[12]). In all of these
cases, two consecutive message interchanges (with relatively high dimensional
hidden states, $64$ and $32$ in our experiments) are performed before a third
one gets the desired compressed representation. We have used the available
implementations of PyTorch Geometric [10] for the convolutional GNN layers,
and performed an exhaustive hyperparameter-tuning for each of them (testing
different hidden dimensions, aggregations, normalizations, dropout values,
number of heads, etc.). As stated in Section 3, in each dataset/compression
ratio scenario we select the best performing model among this set of GNN
architectures to perform the evaluation (shown in Table 1). However, we note
that there is not a significant difference in performance among the different
GNN models, and within a single model different hyperparameter settings do not
result in huge performance variations either; we will further expand this
analysis in future work.

  * •

MPNN: In this case we implement a custom Message Passing GNN whose pipeline
resembles that of the CombMP architecture: edges to edges, edges to nodes,
nodes to edges, and a final edge to node communication with a residual
connection to the original node signal that performs the compression. In this
case the intermediate node and edge hidden states’ dimension is set to
$d^{\prime}=20$, just as in our topological-inspired methods. Notably, the
implementation of this baseline follows directly from the topological
architectures, restricting everything to the graph domain. As a result, it
underwent hyperparameter tuning in exactly the same way as CombMP and SetMP.

In both cases, the final hidden states obtained from the last MP step
represent the node signal compressed representations. Since the original
window-based signals have length $d=10$, we set this final dimension to $4$
and $7$ when benchmarking our method against them for getting compression
factors $r_{c}$ of $1/3$ and $2/3$, respectively. Finally, we note that both
graph-based baselines implement a MLP for the decompression task totally
analogous to the one of our TDL-based architectures (in this case simply
having as input the final node compressed representation).

##### MLP

We also implemented a MLP auto-encoder architecture that considers all
possible connections between all elements of each subsignal
$S_{i}=\\{x_{j}\\}_{j=1}^{N}$. In particular, each of the subsignals is
flattened –i.e. $d_{input}=N\cdot d$– and passed to a feed forward encoder
network (with $1024$ and $512$ hidden dimensions and ReLu activation function
after our architecture search) that outputs a final compressed representation
of the full subsignal with dimension $d_{output}=\lceil r_{c}\cdot
d_{input}\rceil$, being $r_{c}$ the considered compression factor. A symmetric
decoder network reconstructs the signal from that representation.

####  A.4.4 Training and Validation Pipeline

All models are trained for a maximum of $200$ epochs in Abilene dataset and
$500$ in Geant using Adam optimizer (with learning rate $0.003$ and epsilon
$0.001$) and using batches of $25$ samples if required. The model state
corresponding to the best performing iteration over the validation set is
selected for the test evaluation, whose results are represented in Table 1 of
the main body of the paper.

####  A.4.5 Comparison against zfp

Table 2: Reconstruction Mean Squared Error (MSE) and Mean Absolute Error (MAE)
over the test set obtained by our best performing architecture (SetMP) and the
state-of-the-art method zfp.

| Abilene | Geant
---|---|---
$r_{c}=1/3$ | $r_{c}=2/3$ | $r_{c}=1/3$ | $r_{c}=2/3$
MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE
SetMP | $3.22\cdot 10^{-4}$ | $8.75\cdot 10^{-3}$ | $2.03\cdot 10^{-4}$ | $6.80\cdot 10^{-3}$ | $6.93\cdot 10^{-4}$ | $1.52\cdot 10^{-2}$ | $2.90\cdot 10^{-4}$ | $1.05\cdot 10^{-2}$
zfp | $9.19\cdot 10^{-5}$ | $7.34\cdot 10^{-3}$ | $4.02\cdot 10^{-7}$ | $4.84\cdot 10^{-4}$ | $1.04\cdot 10^{-4}$ | $7.83\cdot 10^{-3}$ | $4.18\cdot 10^{-7}$ | $4.95\cdot 10^{-4}$

Finally, we extend the evaluation by showing the relative performance of our
best behaved TDL-based methodology –SeMP– with respect to the state-of-the-art
lossy compression method zfp for the desired precision; see Table 2. As it can
be seen, zfp gets better reconstruction errors than our model in all
scenarios, although we note that differences are considerably lower for
smaller (i.e. more challenging) compression factors. Overall, we reckon that
these are promising results; our topological models are postulated as strong
ML-based baselines for lossy compression, and by addressing some of their
current limitations (see Section A.5) there may be room for for shortening the
gap w.r.t zfp.

In particular, we hypothesize that the main limitation of our current method
revolves around its topology inference procedure, as it can potentially gather
elements that do not necessarily share any correlation –especially at the end
of the clustering process, where groups are built in a greedy manner
regardless of the actual similarities between their elements. Whereas zfp
implements a sophisticated method to deal with such uncorrelated elements, our
current methodology mainly relies on existing correlations to perform
compression; this can lead to a systematic poor reconstruction error of the
involved signals. We will delve deeper into this in the extension of this
work.

###  A.5 Future Work

Despite the promising preliminary results obtained, and because of them as
well, there are many aspects and limitations of our proposed methodology that
are being investigated and will be addressed in future work. The following
list summarizes some of the main lines of research:

##### Topology Inference

Apart from exploring other metrics beyond SNR, it would be interesting to
consider more flexible clustering approaches that could dynamically adapt the
length/number of the inferred higher-order structures (e.g. by defining a
Reinforcement Learning pipeline to perform the division, or adapting a
solution like the Differentiable Cell Complex Module[4] to our scenario).
Moreover, we would also like to explore non-greedy clustering techniques such
as Affinity Propagation.

##### Implementation Issues

Our current proposal does not scale well with the original (graph) signal
dimension, as it mainly relies on an iterative pair-wise similarity
computation between all (graph) elements. More research about how to relax
this aspect is required in order to make its deployment feasible, and in this
regard some ideas we want to explore are:

  * •

To test whether static higher-order structures generated from training samples
can perform well in testing time; not only this would reduce the execution
time, but also the necessity of storing the cluster sequence at each
iteration.

  * •

To check the feasibility of storing simultaneously sparse matrices indicating
when the compression loss exceeds a certain threshold; apart from becoming a
potential indicator of anomalies in the signal, if the distribution of big
reconstruction errors is really sparse, combining them with the compressed
representations will provide with loss tolerance guarantees, and could be key
for matching zfp performance.

  * •

For large graphs, to divide the original graph into independent subgraphs, to
which then our method is applied.

##### Topological MP

Another goal is to shed more light on the performance difference observed
between our two architectures, SetMP and CombMP. Owing to current results, it
seems that the lower complexity of SetMP is a clear advantage, and suggests
that intermediate edge communications noisily interfere in the compression
process. This should be further validated with other datasets, and possibly by
testing some variations of the edge-based communication pipeline in the
general CombMP architecture.