# Does causal dynamics imply local interactions?

Zoltán Zimborás  Wigner Research Centre for Physics, H-1121, Budapest, Hungary
MTA-BME Lendület Quantum Information Theory Research Group, Budapest, Hungary
Mathematical Institute, Budapest University of Technology and Economics,
H-1111, Budapest, Hungary  Terry Farrelly  Institut für Theoretische Physik,
Leibniz Universität, Appelstraße 2, 30167 Hannover, Germany  ARC Centre for
Engineered Quantum Systems, School of Mathematics and Physics, University of
Queensland, Brisbane, QLD 4072, Australia  Szilárd Farkas  Wigner Research
Centre for Physics, H-1121, Budapest, Hungary  Lluis Masanes  Computer Science
Department, University College London, United Kingdom.  London Centre for
Nanotechnology, University College London, United Kingdom.

###### Abstract

We consider quantum systems with causal dynamics in discrete spacetimes, also
known as quantum cellular automata (QCA). Due to time-discreteness this type
of dynamics is not characterized by a Hamiltonian but by a one-time-step
unitary. This can be written as the exponential of a Hamiltonian but in a
highly non-unique way. We ask if any of the Hamiltonians generating a QCA
unitary is local in some sense, and we obtain two very different answers. On
one hand, we present an example of QCA for which all generating Hamiltonians
are fully non-local, in the sense that interactions do not decay with the
distance. On the other hand, we show that all one-dimensional quasi-free
fermionic QCAs have quasi-local generating Hamiltonians, with interactions
decaying exponentially in the massive case and algebraically in the critical
case. We also prove that some integrable systems do not have local, quasi-
local nor low-weight constants of motion; a result that challenges the
standard classification of integrability.

##  I Introduction

Quantum cellular automata (QCAs) originally arose in the context of quantum
computation as the generalization of classical cellular automata Feynman82
and were proven to be universal quantum computers Watrous95 . QCAs can also be
understood as the many-body generalization or “second quantization” of quantum
walks 19Farrelly . From a physics perspective, QCAs are quantum field theories
in discrete spacetimes obeying strict causality SW04 ; SFW06 ; 19Farrelly ;
19Arrighi ; DP14 ; BDPT15 ; GUWZ10 . This means that after one time-step
information only propagates a finite distance. Hence, QCAs provide a rigorous
regularization of (continuous) quantum field theories which simultaneous
preserve causality and unitarity, something impossible in Hamiltonian lattice
field theory 19Farrelly . In Lagrangian lattice field theory, the path
integral is equivalent to a QCA for some field theories FS20 . On a more
speculative level, some arguments suggest that spacetime might be discrete at
the Planck scale, and that all of the more familiar continuous spacetime
physics emerges as an effective description at larger scales. This opens the
possibility of considering QCAs as a Planck-scale theories.

The mathematical formulation discrete-time quantum dynamics is different from
that of continuous time. In the discrete case, dynamics is represented by a
one-time-step unitary evolution operator $W$ and in the continuous case by a
Hamiltonian $H$. The eigenstates of a Hamiltonian $H$ can be ordered with
increasing energies, but the eigenstates of a unitary $W$ cannot be ordered
because the corresponding quasi-energies are defined modulo $2\pi$. This also
makes unclear what should be the Gibbs states associated to $W$. An exception
to this are the unitaries that are close to the identity
$W\approx\mathbbm{1}-\mathrm{i}\epsilon H$, which arise when continuous-time
dynamics is Trotterized for simulations NielsenChuang . But despite the above-
mentioned differences, it is reasonable to expect that, at large time scales,
discrete-time models converge to continuous-time models. The results presented
in this work suggest that this convergence is not straightforward.

In this work we address the following question. If $W$ is the evolution
operator of a QCA, we consider all Hamiltonians $H$ which generate it via
$W=\mathrm{e}^{-\mathrm{i}H}$, and ask whether one of these Hamiltonians is in
some sense local. In general, the Hamiltonian $H$ cannot have finite-range
interactions, because the exponential of a finite-range $H$ is only
approximately causal, as constrained by the Lieb-Robinson bound LR72 . But $H$
can be local in the weaker sense of having interactions that decay with the
distance. In this work we present two extreme examples of QCAs in one spatial
dimension, with opposite decaying behavior.

The first model that we analyze (Section II) is a so called “fractal QCA”
introduced in GUWZ10 . We prove that any of its generating Hamiltonians has
interactions which do not decay with the distance, and that the weight (number
of qubits acted on) of the interaction terms is unbounded. The evolution
operator $W$ of the fractal QCA is a Clifford unitary SW04 , and these share
some features with quasi-free bosonic unitaries, like the fact that dynamics
can be represented in a symplectic phase space of dimension linear in the
number of modes (qubits), which allows for efficiently simulating the dynamics
with a classical computer. But despite sharing these features with integrable
systems, we prove in Section II.2 that all conserved quantities of the fractal
QCA (i.e. operators that commute with $W$) are non-local and have unbounded
weight (like the Hamiltonians). This is a very interesting fact because it
challenges one of the standard characterizations of integrable systems in
terms of local (or low-weight) conserved quantities GME11 . And suggest that,
in the discrete-time scenario, integrability should be characterized
differently.

The second family of QCAs that we analyze (Section III) have general quasi-
free fermion dynamics in one spatial dimension. In this case we show that
there always exists a Hamiltonian with decaying interactions. This decay is
exponential in the gapped case and inversely proportional to the distance in
the critical case. However, we need to define what do we mean by gapped and
critical when (quasi-)energy is defined modulo $2\pi$. We prove that the whole
algebra of operators corresponding to a type of quasi-particle drifts to the
right at a constant speed equal to the winding number of the quasi-energy band
associated to this quasi-particle. Hence, when this winding number is not
zero, the quasi-particle behaves like massless particles in quantum field
theory. For this reason we say that a quasi-free fermionic QCA is critical
when some quasi-energy bands have non-zero winding number. In contrast, when
all winding numbers are zero, we say that the QCA is gapped.

##  II The fractal QCA

###  II.1 Description of the model

Clifford QCAs schlingemann2008structure ; GUWZ10 ; haah2018nontrivial ;
haah2019clifford  are QCAs on lattices of qubits with the property that
products of Pauli operators are mapped to products of Pauli operators. In Ref.
GUWZ10  one-dimensional Clifford QCAs were dived in two classes depending on
the spacetime graph of the evolution of a single-site Pauli operator: those
with a periodic structure (periodic Clifford QCAs) and those with a spacetime
graph that is self-similar over long timescales (fractal Clifford QCAs). This
classification has since then been turned out to have importance in schemes of
measurement-based quantum computation built on Clifford QCAs
stephen2019subsystem  and the fractal property of various QCAs (Clifford and
non-Clifford) has been studied intensely recently GNW10 ; yoshida2013exotic ;
hillberry2020entangled . In what follows we define a particular fractal
Clifford QCA that was studied in GUWZ10 .

Consider a spin chain with $L$ qubits labelled by $r\in\\{0,1,\ldots,L-1\\}$
and periodic boundary conditions. We denote by
$\sigma_{x}^{r},\sigma_{y}^{r},\sigma_{z}^{r}$ the Pauli sigma matrices acting
on qubit $r$. The evolution operator $W$ is determined by conditions

$\begin{split}W^{\dagger}\sigma_{z}^{r}W&=\sigma_{x}^{r}\ ,\\\
W^{\dagger}\sigma_{x}^{r}W&=\sigma_{x}^{r-1}\sigma_{y}^{r}\sigma_{x}^{r+1}\
,\end{split}$ (1)

for all $r$. To see this, recall that
$\sigma_{y}^{r}=\mathrm{i}\,\sigma_{x}^{r}\sigma_{z}^{r}$ and use
$W^{\dagger}\sigma_{y}^{r}W=\mathrm{i}\,W^{\dagger}\sigma_{x}^{r}WW^{\dagger}\sigma_{z}^{r}W$.

To keep track of the evolution of a general $n$-qubit Pauli operator it is
more convenient to use the phase-space description. Then each $n$-qubit Pauli
operator $\sigma_{\mathbf{u}}$ is represented by a phase-space vector
$\mathbf{u}=(q_{0},p_{0},\ldots,q_{L-1},p_{L-1})\in\\{0,1\\}^{2L}$ so that

$\displaystyle\sigma_{\mathbf{u}}=\bigotimes_{r=0}^{L-1}(\sigma_{x}^{r})^{q_{r}}(\sigma_{z}^{r})^{p_{r}}\
.$ (2)

With this notation we have that
$\sigma_{\mathbf{u}}\sigma_{\mathbf{u}^{\prime}}\propto\sigma_{\mathbf{u}+\mathbf{u}^{\prime}}$
where addition in phase space is defined modulo 2. We can write the phase
space of the system as

$\mathcal{V}=\bigoplus_{r=0}^{L-1}\mathcal{V}_{r}\ ,$ (3)

where subspace $\mathcal{V}_{r}\cong\\{0,1\\}^{2}$ is associated to qubit $r$.
For each region of the lattice $\mathcal{R}\subseteq\\{0,\ldots,L-1\\}$ we
define the corresponding subspace of Pauli operators

$\mathcal{V}_{\mathcal{R}}=\bigoplus_{r\in\mathcal{R}}\mathcal{V}_{r}\ .$ (4)

Since Paulis form an operator basis we can express any of the Hamiltonians of
the fractal QCA $W=\mathrm{e}^{-\mathrm{i}H}$ as

$H=\sum_{\mathbf{u}\in\mathcal{V}}h_{\mathbf{u}}\,\sigma_{\mathbf{u}}\ ,$ (5)

where $h_{\mathbf{u}}$ are the real coefficients. We say that a Hamiltonian is
local if the coefficients $h_{\mathbf{u}}$ decay with the size of the support
of $\mathbf{u}$. That is, there is a monotonically decreasing function
$f:\mathbb{N}\to\mathbb{R}^{+}$ such that

$|h_{\mathbf{u}}|\leq f(D(\mathbf{u}))\ ,$ (6)

where the diameter of $\mathbf{u}$ is

$D(\mathbf{u})=\min\\{d(r_{2},r_{1}):\mathbf{u}\in\mathcal{V}_{[r_{1},r_{2}]}\\}\
,$ (7)

with

$d(r_{2},r_{1})=\begin{cases}r_{2}-r_{1}\ \mathrm{if}\ r_{2}-r_{1}\geq 0\\\
r_{2}-r_{1}+L\ \mathrm{otherwise.}\end{cases}$ (8)

###  II.2 No local Hamiltonians

The following lemma tells us that, for the fractal QCA defined in equation
(1), _none_ of the Hamiltonians satisfying $W=\mathrm{e}^{-\mathrm{i}H}$ is
local in the sense of (6).

######  Lemma 1.

Let $W$ be the QCA defined in equation (1), and let
$[r_{1},r_{2}]\subseteq\\{0,\ldots,L-1\\}$ be any interval of the spin chain.
For each Pauli operator in the interval
$\mathbf{u}\in\mathcal{V}_{[r_{1},r_{2}]}$ there is another Pauli operator in
the larger interval $\mathbf{u}^{\prime}\in\mathcal{V}_{[r_{1}-1,r_{2}+1]}$
which is not in any smaller interval and has the same coefficient
$h_{\mathbf{u}^{\prime}}=h_{\mathbf{u}}$ in _any_ of the Hamiltonians $H$
satisfying $W=\mathrm{e}^{-\mathrm{i}H}$, .

###### Proof.

Take _any_ Hamiltonian $H$, such that $W=e^{-iH}$. Then, $H$ satisfies
$W^{\dagger n}HW^{n}=H$ for any $n\in\mathbb{Z}$. We can expand the
Hamiltonian in terms of Pauli operators, rewriting equation (5), as

$H=\sum_{\mathbf{u}\in\mathcal{O}}h_{\mathbf{u}}\sum_{n=0}^{m(\mathbf{u})-1}W^{\dagger
n}\sigma_{\mathbf{u}}W^{n},$ (9)

where now $\mathcal{O}\subset\mathcal{V}$ is a set of labels of Pauli strings
with one vector for _each closed orbit_ under applying $W$ up to
$m(\mathbf{u})$ times. As a convention, we can choose the
$\mathbf{u}\in\mathcal{O}$ that labels each orbit to be the label
corresponding to the Pauli string with the smallest support in the orbit
(i.e., the Pauli string with the smallest $D(\mathbf{u})$). Here
$m(\mathbf{u})$ is the length of the orbit, meaning that $m(\mathbf{u})$ is
the smallest positive integer such that

$W^{\dagger
m(\mathbf{u})}\sigma_{\mathbf{u}}W^{m(\mathbf{u})}=\sigma_{\mathbf{u}}.$ (10)

Note that $m(\mathbf{u})$ must exist for each orbit because we have a finite
quantum system, and there are only a finite number of Pauli strings, so the
orbits must be closed.

The next step is to show that terms like $W^{\dagger
n}\sigma_{\mathbf{v}}W^{n}$ will generally spread over larger and larger
regions for fractal QCAs. Note that the region $\sigma_{\mathbf{v}}$ is
supported on is not guaranteed to grow for _only a single application_ of $W$.
Take, for example,
$\sigma_{\mathbf{v}}=\sigma_{z}^{r}\otimes\sigma_{y}^{r+1}\otimes\sigma_{y}^{r+2}\otimes\sigma^{r+3}_{z}$.
Then we have
$W^{\dagger}\sigma_{\mathbf{v}}W=\sigma_{y}^{r+1}\otimes\sigma_{y}^{r+2}$.

We first note the fact that this QCA $W$ has no gliders, which is proved in
appendix A. A glider is an operator $\sigma_{\mathbf{v}}$ with the property
that $W^{\dagger}\sigma_{\mathbf{v}}W=S^{\dagger k}\sigma_{\mathbf{v}}S^{k}$,
for some $k\in\mathbb{Z}$, where $S$ is just the unitary that shifts qubits
one step to the left. Next, we use lemma II.15 in GUWZ10 . This shows that, if
a Clifford QCA $W$ has no gliders then $W^{n}$ has no gliders for any
$n\in\mathbb{Z}$.

Consider $\sigma_{\mathbf{v}}$ with support on an interval with $l$ sites with
$l<L-2$. If we consider how $\sigma_{\mathbf{v}}$ evolves over time, it
follows that $W^{\dagger n}\sigma_{\mathbf{v}}W^{n}$ must eventually spread
over an interval larger than $l$ sites. To see this, we can argue by
contradiction: suppose that that $W^{\dagger n}\sigma_{\mathbf{v}}W^{n}$
always remains localised on at most $l$ sites (the region may shift left or
right but the range of sites is at most $l$). Now there are only a finite
number of Pauli strings on $l$ sites ($4^{l}$ many). So at some value of
$n\neq 0$ we must have

$W^{\dagger n}\sigma_{\mathbf{v}}W^{n}=S^{\dagger k}\sigma_{\mathbf{v}}S^{k},$
(11)

for some $k\in\mathbb{Z}$. But this means that $W^{n}$ has a glider, which is
impossible because $W$ has no gliders.

Then we can apply this logic to the orbits in the sum in equation (9).
Therefore, we see that $H$ contains interactions between arbitrarily far
regions _with no decay of interaction strength_ with distance. More precisely,
for any Pauli operator $\sigma_{\mathbf{u}}$ in the interval
$\mathbf{u}\in\mathcal{V}_{[r_{1},r_{2}]}$, there is another operator in the
larger interval $\mathbf{u}^{\prime}\in\mathcal{V}_{[r_{1}-1,r_{2}+1]}$ (and
not in any smaller interval), which has the same coefficient
$h_{\mathbf{u}^{\prime}}=h_{\mathbf{u}}$ in any Hamiltonians $H$ satisfying
$W=e^{-iH}$. ∎

###  II.3 No local integrals of motion

The proof of the above lemma not only applies to operators $H$ such that
$W=\mathrm{e}^{-\mathrm{i}H}$ but to any operator which commutes with $W$.
This shows that any constant of motion of the fractal QCA is non-local in the
same sense that the Hamiltonians.

We finish this section by commenting about the continuous-time dynamics of any
Hamiltonian $H$ that generates the fractal QCA $W=\mathrm{e}^{-\mathrm{i}H}$.
As the lemma tells us, the continuous-time dynamics
$\mathrm{e}^{-\mathrm{i}Ht}$ for $t\in\mathbb{R}$ is fully non-local but it
has a particular type of destructive interference that cancels out all non-
causal effects every time that $t$ reaches an integer value $t\in\mathbb{Z}$.

##  III Quasi-free fermionic QCAs

###  III.1 Fermionic systems

In this section we introduce some formalism for working with general quasi-
free fermionic systems. Consider $N$ fermionic modes with associated Majorana
operators $a_{i}$ with $i\in\\{1,2,\ldots,2N\\}$. These operators are
Hermitian $a_{i}^{\dagger}=a_{i}$ and satisfy the canonical anti-commutation
relations

$\left\\{a_{i},a_{j}\right\\}=2\,\mathbbm{1}\delta_{ij}\ ,$ (12)

where $\mathbbm{1}$ is the identity and $\delta_{ij}$ the Kronecker-delta
function. Instead of Majorana operators we could use the creation and
annihilation ones

$\displaystyle f_{r}^{\dagger}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(a_{2r}-\mathrm{i}\,a_{2r+1}\right),$
$\displaystyle f_{r}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(a_{2r}+\mathrm{i}\,a_{2r+1}\right),$
(13)

for $r=1,\ldots,N$, but Majoranas simplify our expressions with no loss of
generality.

A unitary operator $W$ is quasi-free if it maps each Majorana operator onto a
linear combinations of them

$W^{\dagger}a_{i}W=\sum_{j}O_{ij}\,a_{j}\ .$ (14)

The Hermiticity of $a_{i}$ together with the anti-commutation relations (12)
imply that the matrix $O$ is orthogonal. Imposing that $W$ commutes with the
fermionic parity operator

$Q=\bigotimes_{r}(2f_{r}^{\dagger}f_{r}-1)\ ,$ (15)

implies (Lemma 4 in appendix C) that $O$ has unit determinant $O\in{\rm
SO}(2N)$. Note that particle-number conserving models (where $W$ commutes with
$\sum_{r}f_{r}^{\dagger}f_{r}$) are a small subset of the quasi-free models.

A system of $N$ fermion modes can be mapped to $N$ qubits via the Jordan-
Wigner transformation JW28

$\displaystyle
a_{2r}=(\mbox{$\prod_{s=1}^{r-1}$}\sigma_{z}^{s})\,\sigma_{x}^{r}\ ,$
$\displaystyle
a_{2r+1}=(\mbox{$\prod_{s=1}^{r-1}$}\sigma_{z}^{s})\,\sigma_{y}^{r}\ ,$ (16)

where $\sigma_{x,y,z}^{r}$ denote the Pauli sigma matrices acting on qubit
$r$. Although this representation is not local, the product of an even number
of operators from $\\{f_{r},f^{\dagger}_{r},f_{r+1},f^{\dagger}_{r+1}\\}$ only
acts on qubits $r$ and $r+1$ for systems on an infinite line.

For each $O\in{\rm SO}(2N)$ there is a (non-unique) real antisymmetric matrix
$Z$ such that $O=\mathrm{e}^{Z}$. For any anti-symmetric matrix $A$ we define
$\alpha(A)=\frac{1}{4}\sum_{ij}A_{ij}a_{i}a_{j}$. From the anti-commutation
relations (12) it follows that $[\alpha(A),a_{i}]=-\sum_{j}A_{ij}a_{j}$. Then
using the identity

$\displaystyle\mathrm{e}^{A}B\mathrm{e}^{-A}$
$\displaystyle=\mathrm{e}^{[A,\,\cdot\,]}B$
$\displaystyle=B+[A,B]+\frac{1}{2}[A,[A,B]]+\cdots$ (17)

we arrive at

$\mathrm{e}^{-\alpha(Z)}a_{i}\,\mathrm{e}^{\alpha(Z)}=\sum_{j}(\mathrm{e}^{Z})_{ij}\,a_{j}\
.$ (18)

Therefore, up to a phase, any quasi-free unitary $W$ can be written as

$W=\mathrm{e}^{\alpha(Z)}\ .$ (19)

The Hamiltonian $H=\mathrm{i}\alpha(Z)$ is a possible generator for $W$.

###  III.2 QCAs

Consider a spin chain with $L$ sites labelled by
$r\in\mathbb{Z}_{L}=\\{0,1,\ldots,L-1\\}$ and periodic boundary conditions.
Each site $r$ contains $n$ fermionic modes represented by $2n$ Majorana
operators $a_{r}^{l}$ with $l\in\\{1,\ldots,2n\\}$. Complex linear combination
of Majorana operators can be represented by vectors in
$\mathbb{C}^{L}\otimes\mathbb{C}^{2n}$, where we separate the spatial $r$ and
internal $l$ degrees of freedom. The orthogonal matrix $O$ associated to the
QCA’s evolution operator $W$ via (14) acts on the space
$\mathbb{C}^{L}\otimes\mathbb{C}^{2n}$.

Let $|r\rangle$ be the orthonormal basis for $\mathbb{C}^{L}$ corresponding to
the position. Define the translation (or shift) operator $S$ acting on
$\mathbb{C}^{L}$ via $S\,|r\rangle=|r+1\bmod L\rangle$. The properties of
translation invariance and causality imply that $O$ can be understood as the
dynamics of a discrete-time quantum walk with coin space $\mathbb{C}^{2n}$. It
is well known (see, e.g., FS14 ) that translation-invariance implies the
structure

$O=\sum_{q}S^{q}\otimes A_{q}\ ,$ (20)

where $q\in\mathbb{Z}_{L}$ and the operators $A_{q}$ act on the coin space
$\mathbb{C}^{2n}$. In particular, the operator $A_{q}$ specifies how the
information that is translated $q$ sites (to the right) is processed.
Additionally, causality enforces the existence of a neighborhood radius (also
known as interaction range) $R$ beyond which information does not flow after
only one timestep:

$A_{q}=0\mbox{ for all }q\notin[-R,R]\ .$ (21)

###  III.3 The Hamiltonian

Next we obtain the spectral decomposition of (20). The eigenvectors of the
translation operator

$S^{q}|k\rangle=\mathrm{e}^{-iqk}|k\rangle\ ,$ (22)

are the quasi-momentum states

$|k\rangle=\frac{1}{\sqrt{L}}\sum_{r=1}^{L}\mathrm{e}^{ikr}|r\rangle,$ (23)

with $k\in\frac{2\pi}{L}\,\\{0,\ldots,L-1\\}$. If we consider the ansatz
$|k\rangle\otimes|v\rangle$ as an eigenvector of (20) then we obtain

$O\,|k\rangle\\!\otimes\\!|v\rangle=|k\rangle\\!\otimes\\!M_{k}|v\rangle\ ,$
(24)

where we define

$M_{k}=\sum_{q}A_{q}\,\mathrm{e}^{-iqk}\ .$ (25)

The ansatz $|k\rangle\otimes|v\rangle$ is an eigenvector of $O$ if $|v\rangle$
is an eigenvector of $M_{k}$.

The orthogonality of (20) implies the unitarity of $M_{k}$ for all $k$. Hence,
the spectral decomposition

$\displaystyle M_{k}$ $\displaystyle=$
$\displaystyle\sum_{s}\theta_{k}^{s}\,P_{k}^{s}\ ,$ (26)

has complex-phase eigenvectors $\theta_{k}^{s}$ and orthogonal spectral
projectors $P^{s}_{k}$ labelled by $s$.

In all what follows on fermionic QCAs we work in the thermodynamic limit
$L\to\infty$ to maximize the clarity of the results. For finite $L$ the
results are essentially the same but more cumbersome to express. We recall
that in the thermodynamic limit locations are labelled by $r\in\mathbb{Z}$ and
the momentum becomes continuous $k\in[0,2\pi)$.

In this limit, the spectral decomposition of (20) can be written as

$O=\int_{0}^{2\pi}\\!\frac{dk}{2\pi}\ \sum_{s}\theta_{k}^{s}\
|k\rangle\\!\langle k|\otimes P_{k}^{s}\ .$ (27)

Causality (21) implies that the matrix (25) is holomorphic in $k$ (as a
complex variable $k\in\mathbb{C}$). Therefore we can choose the eigenvalues
$\theta^{s}_{k}$ and spectral projections $P^{s}_{k}$ to be holomorphic in $k$
too (see theorem $1.10$ in chapter $2$ of Kato  and also Ahlbrecht12 ;
Agaltsov18 ). Of course, we are primarily interested in the range
$k\in[0,2\pi)$ as opposed to $k\in\mathbb{C}$.

The matrix $M_{k}$ is periodic in $k\in[0,2\pi)$, but this may not be true for
an individual eigenvalue band $\theta_{k}^{s}$. However, we can group bands
together to form continuous and periodic energy bands in the following way.
Suppose the $n$ bands $s_{1},s_{2},\ldots,s_{n}$ form a closed curve

$\displaystyle\lim_{k\to 2\pi}\theta_{k}^{s_{1}}$
$\displaystyle=\theta_{0}^{s_{2}}\ ,$ (28) $\displaystyle\lim_{k\to
2\pi}\theta_{k}^{s_{2}}$ $\displaystyle=\theta_{0}^{s_{3}}\ ,$ (29)
$\displaystyle\vdots$ $\displaystyle\lim_{k\to 2\pi}\theta_{k}^{s_{n}}$
$\displaystyle=\theta_{0}^{s_{1}}\ .$ (30)

Then we can define the periodic holomorphic function $\Theta:[0,2\pi
n)\to\mathbb{C}$ as

$\displaystyle\Theta(\mathbf{k})=\begin{cases}\theta_{\mathbf{k}}^{s_{1}}&\textrm{for}\
\mathbf{k}\in[0,2\pi)\\\ \theta_{\bf k-2\pi}^{s_{2}}&\textrm{for}\
\mathbf{k}\in[2\pi,4\pi)\\\ \ \vdots&\vdots\\\
\theta_{\mathbf{k}-2\pi(n-1)}^{s_{n}}&\textrm{for}\
\mathbf{k}\in[2\pi(n-1),2\pi n)\end{cases},$ (31)

and the periodic holomorphic projector

$\displaystyle\Pi(\mathbf{k})=\begin{cases}P_{\mathbf{k}}^{s_{1}}&\textrm{for}\
\mathbf{k}\in[0,2\pi)\\\ P_{\mathbf{k}-2\pi}^{s_{2}}&\textrm{for}\
\mathbf{k}\in[2\pi,4\pi)\\\ \ \vdots&\vdots\\\
P_{\mathbf{k}-2\pi(n-1)}^{s_{n}}&\textrm{for}\ \mathbf{k}\in[2\pi(n-1),2\pi
n)\end{cases}.$ (32)

Each closed curve is associated to a type of quasi-particle labelled by $\nu$.
Quasi-particle $\nu$ is characterized by the objects
$n^{\nu},\Theta^{\nu}(\mathbf{k}),\Pi^{\nu}(\mathbf{k})$, which also carry the
label $\nu$. This allows to write (27) as

$O=\sum_{\nu}\int_{0}^{2\pi n^{\nu}}\\!\frac{d\mathbf{k}}{2\pi}\
\Theta^{\nu}(\mathbf{k})\
|\mathbf{k}\rangle\\!\langle\mathbf{k}|\otimes\Pi^{\nu}(\mathbf{k})\ .$ (33)

The periodicity of the momentum eigenstates $|k+2\pi\rangle=|k\rangle$ allows
to label them with the extended momentum $\mathbf{k}\in[0,2\pi n^{\nu})$. The
physical interpretation of the extended momentum is the following. Conditions
(28)-(30) for quasi-particle $\nu$ suggest that each lattice site
$r\in\mathbb{Z}$ contains $n^{\nu}$ internal sites and that the dynamics $W$
enjoys a finer translational symmetry (for quasi-particle $\nu$) involving
translations of fractional length $1/n^{\nu}$ which take into account the
extra internal sites. If the number $n^{\nu}$ is the same for all $\nu$ then
we could fine-grain the lattice so that the new system has $n^{\nu}=1$ for all
$\nu$, and the local number of modes is equal to the different types of quasi-
particle (and no more).

Let $w^{\nu}$ denote the winding number of the periodic function
$\Theta^{\nu}:[0,2\pi n^{\nu})\to\mathbb{C}$. This integer is the number of
net loops around the unit circle in $\mathbb{C}$ that the function does across
the interval $[0,2\pi n^{\nu})$. Because $\Theta^{\nu}(\mathbf{k})$ is a
holomorphic function, there is another holomorphic function
$E^{\nu}(\mathbf{k})$ such that

$\Theta^{\nu}(\mathbf{k})=\exp\\!\left(-\mathrm{i}E^{\nu}(\mathbf{k})+\mathrm{i}\frac{w^{\nu}}{n^{\nu}}\mathbf{k}\right)\
.$ (34)

Furthermore, in the range $\mathbf{k}\in[0,2\pi n^{\nu})$ the function
$E^{\nu}(\mathbf{k})$ is periodic and takes real values. Note that these real
values are not restricted to $[0,2\pi)$ due to the continuity imposed in
definition (31).

Let us construct a Hamiltonian whose single-particle energy bands are the
functions $E^{\nu}(\mathbf{k})$. Our choice of matrix $Z$ satisfying
$O=\mathrm{e}^{Z}$ is

$\mathrm{i}Z=\sum_{\nu}\int_{0}^{2\pi
n^{\nu}}\\!\\!\frac{d\mathbf{k}}{2\pi}\left(E^{\nu}(\mathbf{k})-\frac{w^{\nu}}{n^{\nu}}\mathbf{k}\right)|\mathbf{k}\rangle\\!\langle\mathbf{k}|\otimes\Pi^{\nu}(\mathbf{k})\
.$

The reason for writing the quasi-energies as the sum of two terms
$(E^{\nu}(\mathbf{k})-\frac{w^{\nu}}{n^{\nu}}\mathbf{k})$ will be clear below.

Our choice of Hamiltonian $H$ satisfying $W=\mathrm{e}^{-\mathrm{i}H}$ is
$H=\mathrm{i}\alpha(Z)$, which can be written as

$\displaystyle H=$
$\displaystyle\sum_{\nu,r,r^{\prime},l,l^{\prime}}\int_{0}^{2\pi
n^{\nu}}\\!\\!\frac{d\mathbf{k}}{2\pi}\left(E^{\nu}(\mathbf{k})-\frac{w^{\nu}}{n^{\nu}}\mathbf{k}\right)\times$
$\displaystyle\times\,\mathrm{e}^{\mathrm{i}(r-r^{\prime})\mathbf{k}}\,\langle
l|\Pi^{\nu}(\mathbf{k})|l^{\prime}\rangle\,a_{r}^{l}\,a_{r^{\prime}}^{l^{\prime}}\
.$ (35)

In the next section we analyze the locality of this Hamiltonian.

###  III.4 Zero winding implies locality

In this section we consider the case where all energy bands have zero winding
number $w^{\nu}=0$.

In this case the coupling between lattice sites $r,r^{\prime}\in\mathbb{Z}$
specified by Hamiltonian (35) is

$\displaystyle\langle r,l|Z|r^{\prime},l^{\prime}\rangle$ $\displaystyle=$
$\displaystyle-\mathrm{i}\sum_{\nu}\int_{0}^{2\pi
n^{\nu}}\\!\\!\frac{d\mathbf{k}}{2\pi}\,E^{\nu}(\mathbf{k})\,\mathrm{e}^{\mathrm{i}(r-r^{\prime})\mathbf{k}}\,\langle
l|\Pi^{\nu}(\mathbf{k})|l^{\prime}\rangle\,\ ,$ (36)

for any pair $l,l^{\prime}$. Recall that the functions $E^{\nu}(\mathbf{k})$
and $\langle l|\Pi^{\nu}(\mathbf{k})|l^{\prime}\rangle$ are analytic and
periodic in the integration range $\mathbf{k}\in[0,2\pi n^{\nu})$. Therefore,
expression (36) is the Fourier transform of a periodic analytic function. This
is the premise of Lemma 2 from Appendix B, which tells us that

$|\langle r,l|Z|r^{\prime},l^{\prime}\rangle|\leq
C_{1}\,\mathrm{e}^{-\beta_{1}|r-r^{\prime}|}\ ,$ (37)

for some constants $C_{1},\beta_{1}>0$. That is, in the non-critical (gapped
$H$) case interactions decay exponentially with the distance.

###  III.5 Non-zero winding implies weak locality

Suppose the band $\nu$ has non-zero winding number $w^{\nu}\neq 0$. In the
next subsection we see that this can be interpreted as the critical case,
because any quasi-particle of type $\nu$ moves at constant speed
irrespectively of its initial state. Mimicking the behavior of massless
particles in quantum field theory.

The contribution of quasi-particle $\nu$ to the interaction between sites
$r,r^{\prime}\in\mathbb{Z}$ in Hamiltonian (35) is

$\displaystyle\int_{0}^{2\pi
n^{\nu}}\\!\\!\frac{d\mathbf{k}}{2\pi}\left(E^{\nu}(\mathbf{k})-\frac{w^{\nu}}{n^{\nu}}\mathbf{k}\right)\mathrm{e}^{\mathrm{i}(r-r^{\prime})\mathbf{k}}\,\langle
l|\Pi^{\nu}(\mathbf{k})|l^{\prime}\rangle\,\ .$ (38)

The part proportional to $E^{\nu}(\mathbf{k})$ gives an exponential decay as
in (37). The part proportional to $w^{\nu}$ is the Fourier transform of the
product of the analytic function $\langle
l|\Pi^{\nu}(\mathbf{k})|l^{\prime}\rangle$ times the discontinuous (on the
$[0,2\pi n^{\nu})$ torus) function $\mathbf{k}$. The Fourier transform of a
product of two functions is the convolution of their Fourier transforms. The
Fourier transform of the analytic part can be upper-bounded as

$\displaystyle\left|\int_{0}^{2\pi
n^{\nu}}\\!\\!\frac{d\mathbf{k}}{2\pi}\,\mathrm{e}^{\mathrm{i}(r-r^{\prime})\mathbf{k}}\,\langle
l|\Pi^{\nu}(\mathbf{k})|l^{\prime}\rangle\right|\leq
C_{2}\,\mathrm{e}^{-\beta_{2}|r-r^{\prime}|}\ .$ (39)

And the Fourier transform of the discontinuous part is

$\int_{0}^{2\pi
n^{\nu}}\\!\\!\frac{d\mathbf{k}}{2\pi}\,\frac{\mathbf{k}}{n^{\nu}}\,\mathrm{e}^{\mathrm{i}(r-r^{\prime})\mathbf{k}}=\begin{cases}\frac{\mathrm{i}}{r^{\prime}-r}&\mbox{if
}r\neq r^{\prime}\\\ \pi n^{\nu}&\mbox{if }r=r^{\prime}\end{cases}.$ (40)

Hence, their convolution can be upper-bounded by bounding the absolute value
of each term

$\displaystyle\left|\int_{0}^{2\pi
n^{\nu}}\\!\\!\frac{d\mathbf{k}}{2\pi}\frac{w^{\nu}}{n^{\nu}}\mathbf{k}\mathrm{e}^{\mathrm{i}(r-r^{\prime})\mathbf{k}}\,\langle
l|\Pi^{\nu}(\mathbf{k})|l^{\prime}\rangle\right|$ $\displaystyle\leq\ $
$\displaystyle C_{3}\,\sum_{q\neq
0}\,\frac{1}{|q|}\,\mathrm{e}^{-\beta_{2}|r-r^{\prime}-q|}+C_{4}\,\mathrm{e}^{-\beta_{2}|r-r^{\prime}|}$
$\displaystyle\leq\ $ $\displaystyle\frac{C_{5}}{|r-r^{\prime}|}\ ,$ (41)

where $\beta_{2},C_{3},C_{4},C_{5}$ are some constants. Hence, we conclude
that when at least one of the functions $\Theta^{\nu}(\mathbf{k})$ has non-
zero winding number, the Hamiltonian $H$ involves interactions that decay no
slower than the inverse of the distance. This is a much weaker form of
locality than the exponential decay (37). As an example, the massless Dirac
QCA D'Ariano12a ; D'Ariano12b  has non-zero winding numbers and its
Hamiltonian decays as $1/|r-r^{\prime}|$ exactly.

For (single-particle) quantum walks with gapped spectra, effective quasi-local
Hamiltonians were constructed in Osborne08 . In contrast, we get bounds on the
locality of quantum walk Hamiltonians with or without a gap. Note, however,
that Osborne08  had no assumption of translational invariance.

###  III.6 Criticality as drift dynamics

Let us consider a quasi-particle $\nu$ with non-zero winding number
$w^{\nu}\neq 0$. The corresponding sub-algebra of operators is generated by

$\sum_{r,l}\mathrm{e}^{\mathrm{i}\mathbf{k}r}\,\langle
l|\Pi^{\nu}(\mathbf{k})|v\rangle\,a_{r}^{l}\ ,$ (42)

for all $|v\rangle\in\mathbb{C}^{2n}$. If the projector has rank one
$\Pi^{\nu}(\mathbf{k})=|v^{\nu}(\mathbf{k})\rangle\\!\langle
v^{\nu}(\mathbf{k})|$ then we can write the simpler expression

$b^{\nu}(\mathbf{k})=\sum_{r,l}\mathrm{e}^{\mathrm{i}\mathbf{k}r}\,\langle
l|v^{\nu}(\mathbf{k})\rangle\,a_{r}^{l}\ ,$ (43)

for the generators of the sub-algebra of quasi-particle $\nu$. By
construction, the time evolution of these generators is

$W^{\dagger}b^{\nu}(\mathbf{k})W=\mathrm{e}^{-\mathrm{i}E^{\nu}(\mathbf{k})}\,\mathrm{e}^{\mathrm{i}\frac{w^{\nu}}{n^{\nu}}\mathbf{k}}\,b^{\nu}(\mathbf{k})\
.$ (44)

The first phase $\mathrm{e}^{-\mathrm{i}E^{\nu}(\mathbf{k})}$ corresponds the
dynamics generated by an (exponentially) local Hamiltonian. The second term
$\mathrm{e}^{\mathrm{i}\frac{w^{\nu}}{n^{\nu}}\mathbf{k}}$ corresponds to a
spatial translation $r\mapsto r-\frac{w^{\nu}}{n^{\nu}}$ for all the algebra
of operators of quasi-particle $\nu$. This implies that, irrespective of its
initial state, the quasi-particle $\nu$ drifts at a constant speed
$\frac{w^{\nu}}{n^{\nu}}$. Mimicking the behavior of massless particles in
quantum field theory.

###  III.7 Discussion

Our approach for fermions was different but superficially similar to taking a
quantum walk and applying fermionic second quantization 19Farrelly . While our
results include such models (which preserve particle number), they are more
general.

On another topic, the sum of all winding numbers is the index of the
corresponding quantum walk

$\mathcal{I}=\sum_{\nu}w^{\nu}\ ,$ (45)

defined in GNVW12 . This can be interpreted as the net amount of information
flow along the chain. A quantum walk with non-zero index has gapless spectrum
CGGSVWW16 ; asboth2012symmetries .

##  IV Outlook

This work gives rise to the following important open questions. Do QCAs have a
continuous-time limit? How should we describe this limit? Our result suggest
that this is a complex problem which may not have a universal solution. On
another topic, how should we define integrable dynamics (as opposed to chaotic
dynamics) in QCAs? Our results suggest that the existence of local or low-
weight constants of motion is not a good criterion for integrability.

### Acknowledgment

The authors would like to thank Tobias J. Osborne for useful discussions. LM
acknowledges financial support by the UK’s Engineering and Physical Sciences
Research Council (grant number EP/R012393/1). TF was supported by the ERC
grants QFTCMPS and SIQS, the cluster of excellence EXC201 Quantum Engineering
and Space-Time Research, the DFG through SFB 1227 (DQ-mat), and the Australian
Research Council Centres of Excellence for Engineered Quantum Systems (EQUS,
CE170100009). ZZ acknowledges support from the János Bolyai Research
Scholarship, the UKNP Bolyai+ Grant, and the NKFIH Grants No. K124152, K124176
KH129601, K120569, and from the Hungarian Quantum Technology National
Excellence Program, Project No. 2017-1.2.1-NKP-2017-00001.

## References

  * [1] R. P. Feynman.  Simulating physics with computers.  International Journal of Theoretical Physics, 21:467–488, 1982\.
  * [2] J. Watrous.  On one-dimensional quantum cellular automata.  Proceedings of the 36th Annual IEEE Symposium on Foundations of Computer Science, pages 528–537, 1995.
  * [3] T. Farrelly.  A review of quantum cellular automata.  arXiv:1904.13318, 2019.
  * [4] B. Schumacher and R. F. Werner.  Reversible Quantum Cellular Automata.  arXiv:0405174, 2004.
  * [5] R. F. Werner D. J. Shepherd, T. Franz.  Universally programmable quantum cellular automaton.  Phys. Rev. Lett., 97:020502, 2006.
  * [6] P. Arrighi.  An overview of quantum cellular automata.  Natural Computing, 2019.  https://doi.org/10.1007/s11047-019-09762-6.
  * [7] G. M. D’Ariano and P. Perinotti.  Derivation of the Dirac equation from principles of information processing.  Phys. Rev. A, 90:062106, 2014.
  * [8] A. Bisio, G. M. D’Ariano, P. Perinotti, and A. Tosini.  Weyl, Dirac and Maxwell quantum cellular automata.  Foundations of Physics, 45(10):1203–1221, 2015.
  * [9] J. Gütschow, S. Uphoff, R. F. Werner, and Z. Zimborás.  Time asymptotics and entanglement generation of clifford quantum cellular automata.  Journal of Mathematical Physics, 51(1):015203, 2010.
  * [10] T. Farrelly and J. Streich.  Discretizing quantum field theories for quantum simulation.  arXiv:2002.02643, 2020.
  * [11] M. A. Nielsen and I. L. Chuang.  Quantum Computation and Quantum Information.  Cambridge University Press, 2010.
  * [12] E. H. Lieb and D. W. Robinson.  The finite group velocity of quantum spin systems.  Communications in Mathematical Physics, 28(3):251–257, 1972.
  * [13] C. Gogolin, M. P. Müller, and J. Eisert.  Absence of thermalization in nonintegrable systems.  Phys. Rev. Lett., 106:040401, 2011.
  * [14] D. M. Schlingemann, H. Vogts, and R. F. Werner.  On the structure of clifford quantum cellular automata.  Journal of Mathematical Physics, 49(11):112104, 2008.
  * [15] J. Haah, L. Fidkowski, and M. B. Hastings.  Nontrivial quantum cellular automata in higher dimensions.  arXiv preprint arXiv:1812.01625, 2018.
  * [16] J. Haah.  Clifford quantum cellular automata: Trivial group in 2d and witt group in 3d.  arXiv preprint arXiv:1907.02075, 2019.
  * [17] D. T. Stephen, H. P. Nautrup, J. Bermejo-Vega, J. Eisert, and R. Raussendorf.  Subsystem symmetries, quantum cellular automata, and computational phases of quantum matter.  Quantum, 3:142, 2019.
  * [18] J. Gütschow, V. Nesme, and R. F. Werner.  The fractal structure of cellular automata on abelian groups.  Discrete Mathematics & Theoretical Computer Science, DMTCS Proceedings vol. AL, Automata 2010 - 16th Intl. Workshop on CA and DCS, 2010\.
  * [19] B Yoshida.  Exotic topological order in fractal spin liquids.  Physical Review B, 88(12):125122, 2013.
  * [20] L. E. Hillberry, M. T. Jones, D. L. Vargas, P. Rall, N. Y. Halpern, N. Bao, S. Notarnicola, S. Montangero, and L. D. Carr.  Entangled quantum cellular automata, physical complexity, and goldilocks rules.  arXiv preprint arXiv:2005.01763, 2020.
  * [21] P. Jordan and E. Wigner.  Über das Paulische Äquivalenzverbot.  Z. Physik, 47:631–651, 1928.
  * [22] T. C. Farrelly and A. J. Short.  Discrete spacetime and relativistic quantum particles.  Phys. Rev. A, 89:062109, 2014.
  * [23] T. Kato.  Perturbation theory for linear operators.  Springer, 2nd edition, 1995.
  * [24] A. Ahlbrecht.  Asymptotic behavior of decoherent and interacting quantum walks.  PhD thesis, Leibniz Universität Hannover, 2012.
  * [25] A. Agaltsov.  Eigenvalues of analytic families of operators, 2018.  https://www2.mps.mpg.de/homes/ agaltsov/notes/holomkato76.html.
  * [26] G. M. D’Ariano.  The quantum field as a quantum computer.  Physics Letters A, 376(5):697–702, 2012.
  * [27] G. M. D’Ariano.  Physics as quantum information processing: quantum fields as quantum automata.  Foundations of Probability and Physics - 6, AIP Conf. Proc., page 1424 371, 2012.
  * [28] T. J. Osborne.  Approximate locality for quantum systems on graphs.  Phys. Rev. Lett., 101:140503, 2008.
  * [29] D. Gross, V. Nesme, H. Vogts, and R. F. Werner.  Index theory of one dimensional quantum walks and cellular automata.  Communications in Mathematical Physics, 310:419–454, 2012.
  * [30] C. Cedzich, T. Geib, F. A. Grünbaum, C. Stahl, L. Velázquez, A. H. Werner, and R. F. Werner.  The topological classification of one-dimensional symmetric quantum walks.  Annales Henri Poincaré, 19(2):325–383, 2018.
  * [31] B. Tarasinski, J. K. Asbóth, and J. P. Dahlhaus.  Scattering theory of topological phases in discrete-time quantum walks.  Physical Review A, 89(4):042327, 2014.
  * [32] S. G. Krantz and H. R. Parks.  A Primer of Real Analytic Functions.  A Primer of Real Analytic Functions. Birkhäuser Boston, 2002.
  * [33] C. D. Offner.  A little harmonic analysis.  http://www.cs.umb.edu/~offner/files/harm_an.pdf.
  * [34] H. Brezis and L. Nirenberg.  Degree theory and BMO; part I: Compact manifolds without boundaries.  Selecta Mathematica, 1(2):197–263, 1995.

##  Appendix A The fractal QCA has no gliders

To simplify notations, we shall visualize the action of $W$ on monomials of
Pauli-matrices as follows. The matrices $\sigma_{x}$, $\sigma_{y}$, and
$\sigma_{z}$ themselves will be denoted by x, y, and z, respectively.
Monomials of the Pauli matrices obtained by successive applications of $W$
will be written under each other. For example, the action of $W$ on the sigma
matrices can be written schematically as

| x |
---|---|---
x | y | x
| y |
---|---|---
x | z | x
z,
---
x,

and the inverse of $W$ is

| y |
---|---|---
z | x | z
| z |
---|---|---
z | y | z
x.
---
z.

Since they are inessential for the argument, signs are ignored by our
notations, so a monomial is denoted by the same string as minus the same
monomial.

Our proof for the absence of gliders consists of three steps. First we show
that if there were gliders, their length could not change under the time
evolution. Then we demonstrate that this property makes it possible to define
even simpler gliders, which we will call rigid gliders. Finally, a case by
case study rules out the existence of rigid gliders.

(1) The neighborhood of the left end of a string evolves in the following way:

z | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$
---|---|---|---|---|---
| z | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$
|  | z | $\bullet$ | $\bullet$ | $\bullet$
|  |  | y | $\bullet$ | $\bullet$
|  | x | $\bullet$ | $\bullet$ | $\bullet$
| x | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$
x | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$
---|---|---|---|---|---
| $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$
|  | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$
|  |  | x | $\bullet$ | $\bullet$
|  | x | $\bullet$ | $\bullet$ | $\bullet$
| x | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$
x | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$
z | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ |
---|---|---|---|---|---|---
| z | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ |
|  | z | $\bullet$ | $\bullet$ | $\bullet$ |
|  |  | z | $\bullet$ | $\bullet$ | $\Longleftarrow$
|  | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ |
| $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ |
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ |

The arrow indicates the time at which we specify the leftmost operator, and
the three patterns correspond to the three possible choices. The bullets
($\bullet$) indicates that any of the operators x, y, or z may be assigned to
the particular position, but it is also possible that the identity is assigned
to it. So where the frontier consists of bullets, we actually do not know
exactly where the frontier lies. Similarly, the neighborhood of the right
frontier looks like this:

$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | z
---|---|---|---|---|---
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | z |
$\bullet$ | $\bullet$ | $\bullet$ | z |  |
$\bullet$ | $\bullet$ | y |  |  |
$\bullet$ | $\bullet$ | $\bullet$ | x |  |
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | x |
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | x
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$
---|---|---|---|---|---
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ |
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ |  |
$\bullet$ | $\bullet$ | x |  |  |
$\bullet$ | $\bullet$ | $\bullet$ | x |  |
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | x |
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | x
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | z |
---|---|---|---|---|---|---
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | z |  |
$\bullet$ | $\bullet$ | $\bullet$ | z |  |  |
$\bullet$ | $\bullet$ | z |  |  |  | $\Longleftarrow$
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ |  |  |
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ |  |
$\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ | $\bullet$ |

We would like to match the left and right frontiers so that they enclose the
evolution of a glider. For example, take a glider whose left frontier contains
a z. The past of this site is represented by a z-sequence extending
indefinitely in the left upward direction. On the right hand side, this has to
be matched with a parallel x-sequence, which in turn implies that the
z-sequence on the left hand side has to extend ad infinitum in the future as
well. This way we obtain a right-moving glider. Assuming that the left
frontier contains an x, we get a left-moving glider. Clearly, no glider can
have a y in the frontier. So the only possible gliders evolve as

x $\mathrm{S}_{1}$ z z $\mathrm{S}_{1}$ x

x $\mathrm{S}_{2}$ z z $\mathrm{S}_{2}$ x

x $\mathrm{S}_{3}$ z z $\mathrm{S}_{3}$ x,

where $S_{1},S_{2},\dots$ are potentially different strings of equal length.

(2) Suppose that there is a glider for which
$\mathrm{S}_{1}\neq\mathrm{S}_{2}$. Then there is also a shorter glider.
Consider for example a left-moving glider. Evolve it by one time step,
translate it by one lattice site to the right, and multiply the result by the
original operator. If $*$ denotes the product of the corresponding operators,
so that x$*$x, y$*$y, and z$*$z are the empty strings, then the new glider
evolves in the following way:

(x$*$x) $(\mathrm{S}_{1}\\!*\\!\mathrm{S}_{2})$ (z$*$z)=
$(\mathrm{S}_{1}\\!*\\!\mathrm{S}_{2})$

(x$*$x) $(\mathrm{S}_{2}\\!*\\!\mathrm{S}_{3})$ (z$*$z) =
$(\mathrm{S}_{2}\\!*\\!\mathrm{S}_{3})$

(x$*$x) $(\mathrm{S}_{3}\\!*\\!\mathrm{S}_{4})$ (z$*$z)=
$(\mathrm{S}_{3}\\!*\\!\mathrm{S}_{4})$.

So the new glider is shorter at least by two than the original one. Since the
evolution of this glider must follow the previously derived pattern, the left
and right frontiers are again x and z, respectively. We can redefine
$\mathrm{S}_{i}$ as the strings representing the inner part of the new glider.
If these strings change in time, we can repeat the procedure, thereby
obtaining an even shorter glider. After some iterations, we eventually get a
rigid glider, for which $\mathrm{S}_{1}=\mathrm{S}_{2}$, so it is simply
translated by the time evolution.

(3) To rule out rigid gliders, it is enough to look at the operator to the
right of the frontier and see what we get after one time step:

| x | x | $\bullet$
---|---|---|---
x | z | $\bullet$ | $\bullet$
| x | y | $\bullet$
---|---|---|---
x | z | $\bullet$ | $\bullet$
| x | z | $\bullet$
---|---|---|---
x | y | $\bullet$ | $\bullet$
| x |  | $\bullet$
---|---|---|---
x | y | $\bullet$ | $\bullet$

None of these are rigid. An argument analogous to (2) and (3) shows that there
are no right-moving gliders either. (This already follows from the non-
existence of left-moving gliders because the time evolution is invariant under
reflection.)

##  Appendix B Fourier lemmas

The following two lemmas are proved in References [32, 33] and [34]
respectively.

######  Lemma 2.

Let $f:\mathcal{T\to T}$ be an analytic function on the torus with Fourier
decomposition

$f(k)=\sum_{r\in\mathbb{Z}}\hat{f}(r)\,\mathrm{e}^{\mathrm{i}rk}\ .$ (46)

Then there exist two positive constants $C,\beta>0$ such that

$|\hat{f}(r)|<C\,e^{-\beta|r|}\ .$ (47)

######  Lemma 3.

If $f:\mathcal{T\to T}$ is an analytic function with winding number $w$ then

$\sum_{r\in\mathbb{Z}}r\,|\hat{f}(r)|^{2}=w\ .$ (48)

##  Appendix C Determinant of $O$

######  Lemma 4.

A quasi-free fermionic unitary commutes with the parity operator (15) if and
only if the corresponding orthogonal matrix (14) has unity determinant.

###### Proof.

First, note that the parity operator (15) can be written as

$Q=a_{1}a_{2}\cdots a_{2N}\ .$ (49)

Second, note that

$W^{\dagger}a_{1}\cdots a_{2N}W=\sum_{i_{1},\ldots,i_{2N}}O_{1,i_{1}}\cdots
O_{2N,i_{2N}}\,a_{i_{1}}\cdots a_{i_{2N}}.$

Using the anti-commutation relations we can write

$a_{i_{1}}\cdots a_{i_{2N}}=\varepsilon_{i_{1}\ldots i_{2N}}a_{1}\cdots
a_{2N}\ ,$ (50)

where $\varepsilon_{i_{1}\ldots i_{2N}}$ is the Levi-Civita symbol. Combining
all of the above we obtain

$1=\sum_{i_{1},\ldots,i_{2N}}O_{1,i_{1}}\cdots
O_{2N,i_{2N}}\,\varepsilon_{i_{1}\ldots i_{2N}}=\textrm{det}(O)\ .$ (51)

∎


# An Adaptive Boosting Technique to Mitigate Popularity Bias in Recommender
System

Ajay Gangwar  Indian Institute of Technology, RoparIndia
<EMAIL_ADDRESS>and  Shweta Jain  Indian Institute of Technology,
RoparIndia<EMAIL_ADDRESS>
###### Abstract.

The observed ratings in most recommender systems are subjected to popularity
bias and are thus not randomly missing. Due to this, only a few popular items
are recommended, and a vast number of non-popular items are hardly
recommended. Not suggesting the non-popular items lead to fewer products
dominating the market and thus offering fewer opportunities for creativity and
innovation. In the literature, several fair algorithms have been proposed
which mainly focused on improving the accuracy of the recommendation system.
However, a typical accuracy measure is biased towards popular items, i.e., it
promotes better accuracy for popular items compared to non-popular items. This
paper considers a metric that measures the popularity bias as the difference
in error on popular items and non-popular items. Motivated by the fair
boosting algorithm on classification, we propose an algorithm that reduces the
popularity bias present in the data while maintaining accuracy within
acceptable limits. The main idea of our algorithm is that it lifts the weights
of the non-popular items, which are generally underrepresented in the data.
With the help of comprehensive experiments on real-world datasets, we show
that our proposed algorithm outperforms the existing algorithms on the
proposed popularity bias metric.

Recommender systems, Popularity Bias, Fairness

††copyright: none††journalyear: 2021††doi:
10.1145/1122445.1122456††conference: 4th FAccTRec Workshop: Responsible
Recommendation; 4th FAccTRec Workshop: Responsible Recommendation††booktitle:
;††price: 15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs: Information Filtering
Systems Recommender systems

##  1\. Introduction

Recommendation systems have become an essential part of our lives, from movies
we watch on OTT platforms to shopping on e-commerce websites to the news we
read on the internet. In this information age, where the information is
available in abundance or huge volumes, recommendation systems make our lives
simpler by filtering information for us according to our taste.

Recommender systems recommend relevant items to the users based on the user’s
previous data. The primary goal of recommendation engines is to identify the
items that a particular customer might buy or be interested in based on the
previous ratings. They do so by anticipating the ratings that users would give
to the items, and the more data we have for a user, the more precisely we can
predict their ratings. Netflix (movie recommendation system), Spotify (music
recommender system), and Amazon (product recommendation system) are some
examples of recommendation systems that we see in our daily lives.

Although recommender systems are quite popular, there are concerns regarding
the fairness of these systems amongst the research community. Recommender
systems primarily face fairness issues at two levels, user-level (Ekstrand and
Kluver, 2021; Mansoury et al., 2020) and item level (Saito, 2020; Ovaisi et
al., 2020). This paper considers the most important fairness issue, namely the
popularity bias at the item level. The items rated by most of the users or
have received high ratings are known as popular items, and the items rated by
very few users or not rated at all are known as non-popular items. The
underlying notion behind popularity bias is that people are more inclined to
offer comments on mainstream or popular products than on non-popular items. As
a result, reported user response is skewed towards popular products rather
than genuine user interest. Thus, frequently rated items or popular items
receive a lot of exposure, but less popular or niche items are
underrepresented in a recommendation. Recommending only popular items is
problematic for two reasons: first, not everyone wants to follow the taste of
the mainstream crowd, and second, it makes it harder for new goods to gain
user attention.

One example of how just suggesting popular things affects companies is that a
lawsuit was recently filed in the United States against Google for only
showing advertising for popular items and not giving opportunities to less
popular and newer items on the market (San Ramon, 2020). Because just a few
items are displayed most of the time, this sort of activity generates a
monopoly in the market, which is unsuitable for any firm. It will hamper the
possibilities of creativity and innovation in the products or items.

Exploring and mitigating popularity bias in recommender systems is not a new
problem and have been previously explored in many works (Abdollahpouri et al.,
2019; Huang et al., 2020; Schnabel et al., 2016; Steck, 2011; Abdollahpouri et
al., 2017) . However, the existing works focus either on improving the
accuracy of the overall recommendation system or exploring the non-popular
items through diversification. Since more ratings are available for the
popular items, the accuracy is inherently biased towards the popular items. A
highly accurate recommender system may have very high accuracy on popular
items but a very low accuracy on non-popular items. On the other hand, naive
diversification may lead to poor accuracy of the overall recommender system.
In this work, we define popularity bias as the difference between errors on
non-popular items and popular items. A fair recommender system should perform
equally well on both popular and non-popular items thus ensuring a balanced
accuracy of the recommender system. The main aim is to reduce popularity bias
from the data and recommend relevant items to the user based on prior data
while keeping the error as low as possible.

###  1.1. Contributions

The main contributions of this paper are as follows:

  * •

We provide a metric that quantitatively measures the popularity bias present
in the data. The previous papers have used error as a metric that is
inherently biased towards the popular items and thus may lead to poor accuracy
on non-popular items.

  * •

To reduce the popularity bias, we propose a novel algorithm, namely FairBoost,
that significantly reduces the popularity bias present in the data without
deteriorating on the error.

  * •

We compare our proposed algorithm with the existing algorithms that claim to
remove popularity bias and show that our algorithm outperforms the existing
algorithms on the proposed popularity bias metric.

##  2\. Related Works

Fairness in recommender systems has been the center of discussion for a long
time since these systems suffer from different types of biases present in the
data. Due to these biases, it can not represent the entire population, and due
to this, the results produced by recommender systems are biased, so there is
an utmost need to solve this problem. The most common types of biases are
gender bias (Ekstrand and Kluver, 2021; Mansoury et al., 2020), racial
bias(Obermeyer et al., 2019), selection bias (Saito, 2020; Ovaisi et al.,
2020), exposure bias (Yadav et al., 2019; Singh and Joachims, 2018), position
bias (Collins et al., 2018; Guo et al., 2019), and popularity bias
(Abdollahpouri et al., 2019; Huang et al., 2020; Schnabel et al., 2016; Steck,
2011; Abdollahpouri et al., 2017). The focus of this paper is on popularity
bias.

It is crucial to recommend non-popular items as they are the ones that are
less likely to be discovered. Abdollahpouri et al. (2019) have proposed a
customized re-ranking diversification strategy that aims to boost non-popular
items representation while retaining acceptable recommendation accuracy. The
focus in (Abdollahpouri et al., 2019) was to achieve a trade-off between
coverage of popular and non-popular items. To achieve this coverage, non-
popular items that consumers could appreciate were randomly selected and
presented to the user without being displayed in some order of preference.
Thus the proposed method could result in poor accuracy for non-popular items.

To address the popularity bias, Huang et al. (2020) and Schnabel et al. (2016)
employed the Inverse Propensity Score (IPS) method. Propensity-based
techniques have previously been used in causal inference and observational
studies but are used for the first time in the recommender system by Schnabel
et al. (2016). To apply the IPS approach, we need propensities, and the
authors in (Schnabel et al., 2016) have discussed primarily two propensity
estimation models, via Naive Bayes, and the other is via logistic regression.
Saito (2020) has also discussed a straightforward approach to measure the
propensities. However, the performance of these propensity-based methods is
heavily influenced by the model used to estimate propensity. The algorithms
require actual probabilities to function correctly which are difficult to
estimate (Saito, 2020). As a result, a more effective strategy is required.

Furthermore, all of the above works focus on reducing recommender system
errors without explicitly being fair to non-popular items.Two exceptions to
this are (Boratto et al., 2021; Zhu et al., 2020). (Boratto et al., 2021)
mitigate the issue of popularity bias by providing a metric similar to the
individual fairness metric used in machine learning. The idea is to equalise
the outcomes across all individual items. On the other hand, (Zhu et al.,
2020) focuses on group fairness by grouping the items together and equalizing
the true positives across all groups. Our work focuses on group fairness, and
instead of just equalizing true positives, we sought to equalize the total
error across the groups.

Machine learning algorithms typically discriminate based on biased historical
data. Data manipulation is used in certain preprocessing techniques to address
this discrimination problem. Massaging (Kamiran and Calders, 2009) , re-
weighting (Calders et al., 2009) , uniform or preferential sampling (Kamiran
and Calders, 2012) , and data augmentation (Iosifidis and Ntoutsi, 2018) are
examples of preprocessing methods that might have an impact on data
distribution. By altering instance labels, giving different weights, under or
oversampling instances, and producing pseudo instances, they attempt to
correct imbalances between protected and non-protected groups within the data.
Iosifidis and Ntoutsi (2019) demonstrated Adafair, an Adaboost-based fairness-
aware classifier that adjusts the weights of the instances in each subsequent
round while explicitly addressing class imbalance by optimizing the number of
ensemble models for balanced classification error. Adapting this idea, we
examine the use of boosting algorithm in the recommender system. The idea is
to increase the weights of non-popular items to maintain the balance between
popular and non-popular items. We show that this technique can significantly
reduce the popularity bias present in the system.

##  3\. The Model

Let us consider the set of users ${U}$ and the set of items ${M}$, with the
number of users $|U|=k$ and the number of items $|M|=l$. $A=U\times{M}$
represents the collection of (user, item) pairs with each entry $A_{u,m}$
denoting the true rating provided by user $u$ to item $m$. Integer values
ranging from $1$ to $5$ can be used for representing each entry, with $1$
representing low interest and $5$ indicating strong interest. Because
consumers only rate a small portion of the items, many ratings are typically
unknown. The main aim is to develop an algorithm that generates the best
possible predicted rating matrix $\hat{A}$. In the predicted rating matrix
$\hat{A}$, each entry $\hat{A}_{um}$ denotes the probable rating given by user
$u$ to item $m$. To achieve the goal, ideal loss function to find predicted
rating matrix is defined as:

(1)
$L_{ideal}(\hat{A})=\frac{1}{|A|}\sum_{u=1}^{k}\sum_{m=1}^{l}{\delta_{um}(A_{um},\hat{A}_{um})}$

where $\delta_{um}$ can be taken as mean squared error (MSE) i.e.

$\delta_{um}^{MSE}(A,\hat{A})=(A_{um}-\hat{A}_{um})^{2}$

or minimum absolute error (MAE) or i.e.

$\delta_{um}^{MAE}(A,\hat{A})=|A_{um}-\hat{A}_{um}|$

.

The ideal loss function cannot be calculated since the majority of the
elements in the actual rating matrix are missing. Therefore, to build an
effective recommender system, we need to estimate the ideal loss as accurately
as possible using only the ratings present in the true rating matrix. The
simplest straightforward estimator for estimating ideal loss, also known as a
naïve estimator is used which is defined as follows:

(2)
$L_{naive}(\hat{A})=\frac{\sum_{(u,m):B_{um}=1}\delta_{um}(A_{um},\hat{A}_{um})}{|\\{(u,m):B_{um}=1\\}|}$

where $B$ is an observed rating matrix and each entry $B_{um}=1$ if user $u$
has given rating to item $m$ otherwise $B_{um}=0$.

The loss over observed ratings is calculated using naive loss, which is
defined as the sum of loss overall user-item pairs in the data divided by the
total number of user-item pairs in the true rating matrix. It is easy to see
that when the data are missing fully at random, this naïve estimator is
unbiased i.e.

$E[L_{naive}(\hat{A})]=L_{ideal}(\hat{A})$

However, when the data are not missing completely at random which happens in
the presence of popularity bias where there is more rating available for
popular items as opposed to that of non popular items, it is shown that the
naive loss may not exactly correspond to the ideal loss function (Schnabel et
al., 2016; Steck, 2010) i.e.

(3) $E[L_{naive}(\hat{A})]\neq L_{ideal}(\hat{A})$

One popular approach to mitigate the popularity bias is by using Inverse
Propensity Score (IPS) approach (Huang et al., 2020; Schnabel et al., 2016).
The main idea behind this approach is to build a pseudo missing completely at
random dataset by weighting all the observed ratings by the inverse of their
propensity score. It is easy to theoretically show that IPS loss is an
unbiased estimator and hence can be used to remove popularity bias (Saito,
2020).

The IPS estimator’s unbiasedness is desired; nevertheless, this feature is
dependent on the true propensity scores which need to be approximated using
various approaches. These methods majorly suffers from two problems (Saito et
al., 2017). One, the IPS estimator is no longer an unbiased estimator if a
propensity estimation model is not stated appropriately. Another, the inverse
of the propensities may be substantial, so the IPS estimator has a high
variance problem. As a result, building learning methods that are resilient to
misspecification of propensity and estimator variation is crucial for applying
the approaches to real-world MNAR situations.

Saito (2020) devised an asymmetric tri-training technique based on the
asymmetric tri-training approach used in unsupervised domain adaptation to
tackle the challenges that an IPS approach faced. It employs three rating
predictors, two of which were used to create a pseudo rating dataset and the
third of which was used to train the model on these pseudo ratings. The
problem with this method is that the size of the dataset reduces after the
algorithm is applied, making it impossible to accurately estimate the ratings
of all the items. Because data in MNAR datasets are typically sparse, we must
make appropriate use of it in order to predict ratings. Boosting technique on
the other hand would be ideal because they utilise the majority of their data
in different iteration stages, and one can also use upweighting to boost non-
popular items in the same way that incorrectly categorised points are boosted
in subsequent iterations.

###  3.1. Quantifying popularity bias

So far, the research community has been focusing on mitigating the popularity
bias so as to improve the overall accuracy of the system. As mentioned in the
Introduction section, resolving popularity bias as a standalone problem is
highly required in order to prevent monopoly in the system. A typical accuracy
measure is not a good metric for popularity bias because it is biased towards
popular items, i.e., it promotes better accuracy for popular items compared to
non-popular items. This is because more number of ratings are given to popular
items.

One naive way to resolve popularity bias is to recommend non-popular items
randomly thereby exploring these items. Another possibility to avoid the
monopoly is through diversification (Parapar and Radlinski, 2021; Koutsopoulos
and Halkidi, 2018) which ensures that all group of items must be recommended
some number of times. However, both naive or diversification methods could
potentially lead to dissatisfaction amongst the users of the recommendation
system. Thus, we need a metric that not only ensures the accuracy of the
overall system, but also prevent the issue of popularity bias.

We define the popularity biasedness of any algorithm by the difference in the
error it achieves between the non-popular items and popular items. Let
$\mathcal{PS}$ be the set of popular items and let $\mathcal{NPS}$ be the set
of non-popular items. We employ the threshold $\tau$ to identify which of the
items are popular and which are not. The set of all those items that obtained
more than the threshold number of ratings will be regarded popular, whereas
the set of items that received less than the threshold number of ratings will
be considered non-popular.$\tau$ is a dataset-specific parameter which can be
taken as input. Then the popularity biasedness of the algorithm with predicted
rating matrix $\hat{A}$ is defined as:

(4) $\begin{split}PB(\hat{A},\tau)&=\frac{\sum_{(u,m)\colon
m\in\mathcal{NPS}}\delta_{u,m}(A_{um},\hat{A}_{um})B_{um}}{\sum_{(u,m)\colon
m\in\mathcal{NPS}}B_{um}}-\\\ &\frac{\sum_{(u,m)\colon
m\in\mathcal{PS}}\delta_{u,m}(A_{u,m},\hat{A}_{u,m})B_{um}}{\sum_{(u,m)\colon
m\in\mathcal{PS}}B_{um}}\end{split}$

where $\delta_{u,m}$ can be taken as MSE(mean squared error) or MAE(minimum
absolute error). This popularity bias metric separates the items into two
groups: popular and non-popular. After that, the mean error is computed
separately for non-popular and popular sets. Finally, the difference between
them is used to determine the popularity bias $PB(\hat{A},\tau)$. The goal is
to minimize the $PB(\hat{A},\tau)$ to reduce the effect of popularity bias.

###  3.2. Proposed Algorithm : FairBoost

Boosting is a method of creating a powerful learner by combining several weak
base learners. Adaboost(Schapire, 1999) is one such approach, which
iteratively calls weak base learners after modifying the weights based on
misclassified data points in each iteration. Because it separates the learning
problem into numerous sub-problems and then combines their answers into an
overall model, we believe boosting is a good fit for our problem. In sub-
models, the popularity bias problem is easier to address than in the entire
complicated model. As popular items frequently appear in the recommender
system’s results, the weights of the non-popular items must be increased in a
way such that it shows up in the recommender system’s results. In the
classification setting, the AdaBoost algorithm internally boosts the weights
of erroneously categorised data points. The way we increase the weights of
incorrectly classified data points in Adaboost, we increase the weights of
non-popular items while keeping the accuracy on these items in mind. Finally,
we adjust the reweighting procedure of the AdaBoost algorithm to make it more
fair and obtain the FairBoost Algorithm. Adaboost is not a novel technique for
recommender systems; it was previously used to improve accuracy when using a
decision tree as the base learner (Golbandi et al., 2011). However, we are the
first to use and modify this notion to reduce popularity bias in recommender
systems.

Algorithm 1 depicts the training phase of FairBoost. To establish its
significance in the training dataset, FairBoost also provides weight to each
training example. When the given weights are high, that set of training user-
item pairs is more likely to influence the training set. Similarly, user-item
pairs with low weights will have a little impact on the training dataset. At
first, all of the user-item pairs will be given same weight of $1/n$, where
$n$ is the number of user item pairs. Fairboost additionally employs a
popularity bias-related cost denoted by $cost_{um}$ for each user-item pair
$(u,m)$, which try to maintain similarity and reduce the popularity bias that
exists between popular and non-popular sets for current learners. We
initialize $cost_{um}$ to zero for all $(u,m)$. Then the user item pairs are
sampled using sample weights $w_{um}$, and a weak learner $\hat{A}^{j}$ is
trained on these sampled points. Let $S_{j}$ is the set of sampled user-item
pairs at $j^{th}$ iteration. The error rate $err_{j}$ is computed using:

(5) $err_{j}=\sum_{(u,m)\in
S_{j}}w_{um}\left(1-e^{\frac{-(A_{um}-\hat{A}^{j}_{um})}{\max(A_{um}-\hat{A}^{j}_{um})}}\right)$

where $w_{um}$ is the weight of $(u,m)$ user-item pair in the sampled user-
item pairs set $S_{j}$, $A_{um}$ is the actual rating of user-item pair and
$\hat{A}^{j}_{um}$ is the rating predicted by the current base learner for
user-item pair $(u,m)$. Following that, we use the following formula to
calculate the the weight $\alpha_{j}$ depicting the influence of the base
learner $j$ in predicting ratings.

(6) $\alpha_{j}=\frac{1}{2}\log{\left(\frac{1-err_{j}}{err{j}}\right)}$

After that popularity bias ${PB(\hat{A}^{j},\tau)}$ is computed for the
current base learner using Equation (4).

Next popularity bias related cost, $cost_{um}$ is computed for all the user-
item pairs $(u,m)$ in the sampled user-item pairs. Since we conduct random
sampling every time, popular items become non-popular in some iterations and
non-popular items become popular in other iterations. We seek to achieve
similarity between popular and non-popular sets or eliminate popularity bias.
As a result, in the current iteration, a popularity bias-related cost is
employed to preserve this similarity or reduce popularity bias. $cost_{um}$ is
given by :

(7)
$cost_{um}=\begin{cases}|{PB(\hat{A}^{j},\tau)}|,if\big{[}(A_{um}-\hat{A}^{j}_{um})>\epsilon_{1}AND\\\
\mspace{28.0mu}|{PB(\hat{A}^{j},\tau)}|>\epsilon_{2},m\in\mathcal{PS},{PB(\hat{A}^{j},\tau)}>0\big{]}\\\
|{PB(\hat{A}^{j},\tau)}|,if\big{[}(A_{um}-\hat{A}^{j}_{um})>\epsilon_{1}AND\\\
\mspace{28.0mu}|{PB(\hat{A}^{j},\tau)}|>\epsilon_{2},m\in\mathcal{NPS},{PB(\hat{A}^{j},\tau)}<0\big{]}\\\
0,otherwise\end{cases}$

It should be noted that $\epsilon_{1}$ is a hyperparameter used to set a bound
on the difference between true and predicted ratings, and $\epsilon_{2}$ is a
hyperparameter used to set a bound on popularity bias. Both of these
parameters must be tuned for the algorithm to work properly. For all user-item
pairs, if item belongs to popular set, popularity bias
${PB(\hat{A}^{j},\tau)}$ is greater than zero and if difference between true
and predicted rating is greater than $\epsilon_{1}$ and absolute value of
popularity bias for current base learner is greater than $\epsilon_{2}$, then
$cost_{um}$ assigned to the pair is $|{PB(\hat{A}^{j},\tau)}|$. Here, a
popular item has become unpopular, therefore we are weighing it more. For all
user-item pairs, if item belongs to non-popular set, popularity bias
${PB(\hat{A}^{j},\tau)}$ is less than zero and if difference between true and
predicted rating is greater than $\epsilon_{1}$ and absolute value of
popularity bias for current base learner is greater than $\epsilon_{2}$, then
$cost_{um}$ assigned to the pair is $|{PB(\hat{A}^{j},\tau)}|$. We are
upweighting this item because it is already unpopular. And for all other user-
item pairs $cost_{um}$ of zero is assigned. Then the weights are updated for
the next round after computing costs using :

(8)
$w_{um}\leftarrow\frac{1}{Z_{j}}w_{um}.e^{\alpha_{j}.(A_{um}-\hat{A}^{j}_{um})}.(1+cost_{um})$

where $Z_{j}$ is a factor used for normalizing weights. The user-item pairs
having the larger error are given greater weight, so they can be predicted
accurately in the following iteration. Weights are updated for all sampled
user-item pairs by multiplying $w_{um}$ with exponential of product of current
base learner weight and difference between true and predicted rating so that
the examples that are having more error gets more weight in next round and the
examples having less errors gets less weight in next round. Weights are also
multiplied by (1+$cost_{um}$). It is done for all those examples that were
treated unfairly during current round.

And once the number of rounds is reached, the algorithm converges. The
algorithm will generate a number of base learners equal to the number of
rounds, which will then be merged using the weights $\alpha_{j}$ or the amount
of influence they have to produce the estimator $\hat{A}^{j}$ as an output.
Now, in order to forecast the rating of a new user-item pair, it will go
through all of the base learners. The predicted ratings from all these
learners are weighted with the corresponding weights of the base learners and
then combined to provide the expected rating for the new user-item pair.

Input: A=$(X_{um},Y_{um})^{N},M,\epsilon_{1},\epsilon_{2}$

Output: Estimator $B$

1 Initialize $w_{um}$ = $\frac{1}{N}$ and $cost_{um}=0$,for all user-item pair
$(u,m)\in A$

2 for  _$j\leftarrow 1$ to $M$_ do

3 a) Weak learner $\hat{A}^{j}$ is trained using weights $w_{um}$ on training
data

4 b) Error rate $err_{j}$ is computed

5 c) Weight is computed for the weak learner,
$\alpha_{j}=\frac{1}{2}\log{\left(\frac{1-err_{j}}{err{j}}\right)}$

6 d) Popularity Bias ${PB(\hat{A}^{j},\tau)}$ is computed

7 e) $cost_{um}$ related to popularity bias is computed.

8 f) Distribution is updated as
$w_{um}\leftarrow\frac{1}{Z_{j}}w_{um}.e^{\alpha_{j}.(A_{um}-\hat{A}^{j}_{um})}.(1+cost_{um})$,
where $Z_{j}$ is a factor used for normalizing weights.

$B(x)=\sum_{j=1}^{M}\alpha_{j}\hat{A}^{j}(x)$

Algorithm 1 FairBoost Algorithm to mitigate popularity bias

##  4\. Experimental Analysis

###  4.1. Experimental setup

####  4.1.1. Datasets

To show the effectiveness of our proposed algorithm, we have used the
following realworld datasets. In all the datasets, training set and testing
set have been created by considering 80% and 20% ratings respectively.

  * •

Netflix dataset(29, 2006) The Netflix dataset consists of about 100 million
MNAR five-star movie ratings (missing not at random). There are 480189 users
and 17770 movies involved. Due to computational restrictions, we took 10
million of the most recent ratings by sorting them according to time. There
are 370811 users and 1962 movies in the training set, whereas there are 258603
users and 1962 movies in the test set.

  * •

Yahoo dataset(31, 2008)
This dataset captures the preferences of the Yahoo! Music community and
comprises almost 717 million ratings on 136 thousand songs provided by 1.8
million users. The data was gathered between 2002 and 2006. The ratings are on
a scale of 1 to 5, with 1 being the lowest and 5 being the highest. We have
sampled around 10 Million ratings from this set. Training set consists of
23179 users and 136737 songs while test set consists of 23179 users and 63261
movies.

  * •

Amazon dataset(R. He, 2016)
This is an Amazon Movies and TV dataset which consists of around 4.6 million
ratings. The ratings are on a scale of 1 to 5, with 1 being the lowest and 5
being the highest. There are 2088620 users and 200941 items involved. There
are 1666901 users and 188083 items in the training set, whereas there are
562187 users and 80112 items in the test set.

  * •

Movielens dataset(30, 1998)
This dataset is made up of 100K five-star movie ratings(missing not at random)
gathered from a movie recommendation service . There are 1682 movies and 943
people involved. The data has been sorted by date. The dataset is divided into
a training set and a test set, with the training set consisting of the
previous 80% of user-item pairs and the test set consisting of the rest or
most recent 20% of user-item pairs. There are 751 users and 1616 movies in the
training set, whereas there are 301 users and 1448 movies in the test set.

####  4.1.2. Compared methods

We conducted comprehensive testing on the above mentioned datasets to
demonstrate that our proposed algorithm lowers the system’s popularity bias.
Three algorithms were tested and their results were compared. We compared the
Fairboost algorithm to that of Matrix Factorization with Inverse Propensity
scoring (Schnabel et al., 2016) and Asymmetric tri-training (Saito, 2020). The
previous papers used error as a measure and focused on lowering the error
rather than explicitly addressing the system’s popularity bias. We compared
several methods using the proposed popularity bias metric in Equation ( 4). We
keep the value of $\tau$ to be $100$ which means if an item received rating
from more than $100$ users, then we call that item as popular item.

####  4.1.3. Hyperparameter Tuning

To adjust the parameters $\epsilon_{1}$ and $\epsilon_{2}$, we used a random
search cv hyperparameter tuning procedure. Both parameters were tuned in the
range $[10^{-5},1]$.

###  4.2. Results and Discussion

Datasets | Algorithm | Error | Error on popular items | Error on non-popular items | Popularity Bias
---|---|---|---|---|---
Netflix | Matrix factorization | 1.0896 | 1.0675 | 1.1315 | 0.0639
|  Matrix factorization with IPS (Schnabel et al., 2016) | 1.0798 | 1.236 | 1.297 | 0.0610
|  Asymmetric tri-training (Saito, 2020) | 1.2165 | 1.1076 | 1.1585 | 0.0508
|  Adaboost (Golbandi et al., 2011) | 1.1221 | 1.1087 | 1.1480 | 0.0393
| FairBoost | 1.1481 | 1.1373 | 1.168 | 0.0315
Movielens | Matrix factorization | 1.0268 | 1.0154 | 1.0272 | 0.0117
|  Matrix factorization with IPS (Schnabel et al., 2016) | 1.0199 | 1.1630 | 1.1728 | 0.0098
|  Asymmetric tr-training (Saito, 2020) | 1.1452 | 1.1519 | 1.1450 | 0.0069
|  Adaboost (Golbandi et al., 2011) | 1.10544 | 1.1011 | 1.1055 | 0.0080
| FairBoost | 1.0577 | 1.0478 | 1.0511 | 0.0025
Amazon | Matrix factorization | 1.1468 | 1.2201 | 1.1223 | 0.0977
|  Matrix factorization with IPS (Schnabel et al., 2016) | 1.1432 | 1.2264 | 1.3140 | 0.0876
|  Asymmetric tr-training (Saito, 2020) | 1.1509 | 1.2126 | 1.1304 | 0.0822
|  Adaboost (Golbandi et al., 2011) | 1.1581 | 1.1821 | 1.1502 | 0.0319
| FairBoost | 1.2984 | 1.3286 | 1.3225 | 0.0119
Yahoo | Matrix factorization | 1.4729 | 1.5436 | 1.4700 | 0.0736
|  Matrix factorization with IPS (Schnabel et al., 2016) | 1.3969 | 1.5876 | 1.5203 | 0.0673
|  Asymmetric tr-training (Saito, 2020) | 1.5693 | 1.6111 | 1.5676 | 0.0434
|  Adaboost (Golbandi et al., 2011) | 1.5712 | 1.6107 | 1.5696 | 0.0411
| FairBoost | 1.5021 | 1.5749 | 1.4991 | 0.0327

Table 1. Results

The results from all the datasets are summarised in Table 1 which displays the
results of 10 iterations of all the algorithms. Table 1 can be used to make
the following observations. The popularity bias was reduced when IPS was used
in conjunction with matrix factorization. However, the error on both popular
and non-popular items had increased when compared to the baseline matrix
factorization algorithm on all the datasets. Another finding was that the
propensity estimation model chosen had a significant impact on the performance
of propensity-based unbiased estimation approaches. We have used naive Bayes
propensity estimation method for the comparison as it gave us the least value
of popularity bias. The other methods such as user-item propensity gave us a
good accuracy, however, the popularity bias was quite high. Another
observation is that the popularity bias was also reduced when asymmetric tri-
training was used when compared to matrix factorization based methods, as
shown in table 1.

(a) Netflix dataset

(b) Yahoo dataset

(c) Amazon dataset

(d) Movielens dataset

Figure 1. Number of estimators vs Popularity Bias

(a) Netflix dataset

(b) Yahoo dataset

(c) Amazon dataset

(d) Movielens dataset

Figure 2. Number of estimators vs Error Plot

As popular items appear frequently in the recommender system’s results. So, in
order for non-popular items to show up in the recommender system’s results,
their weightage must be increased in some way. In the classification setting,
the AdaBoost algorithm internally boosts the weightage of wrongly categorised
data points, so we thought we’d give that a try in ours because we also want
to give non-popular items more weightage. Table 1 shows that after applying
Adaboost, the popularity bias has decreased on almost all of the datasets as
compared to previously implemented methods. Inspired by how well Adaboost
performed, we modified it to create the FairBoost algorithm (Algorithm 1).
When running the FairBoost algorithm, we employed matrix factorization as our
underlying base learner. On all the datasets, Table 1 shows that our proposed
algorithm Fairboost significantly reduced the popularity bias when compared to
other algorithms. In the following steps, the FairBoost algorithm tries to
make non-popular items popular. We increase the weights of non-popular items
in our algorithm in the following steps to make them popular, similar to how
Adaboost increases the weights of inaccurately classified data points. Figures
1(a), 1(b) and 1(c) show that the popularity bias decreases as the number of
estimators increases. This is because the FairBoost algorithm inherently gives
more weight to non-popular items, causing the graph to decrease. In the case
of the movielens dataset, we got a zig-zag graph, as shown in Fig 1(d). This
is because the data is nearly fair; as the table 1 shows, the error difference
is extremely small, and ratings are also uniformly distributed. Figures 2(a),
2(b), 2(c) and 2(d) show that the overall inaccuracy increases slightly. This
is due to the fact that in succeeding steps, popular items weights are given
less weight. As a result, FairBoost Algorithm minimises the popularity bias in
each subsequent step while keeping the error increase within acceptable
bounds.

##  5\. Conclusion and Future Work

Adequate coverage of non-popular or long-tail items is critical to any
business’s success. Because almost all users are familiar with popular items,
a recommender system’s ability to recommend non-popular items will determine
how well it introduces users to new experiences and products; however, it is
well known that recommender systems are biased towards popular items.
This paper proposed and compared FairBoost, an adaptive boosting algorithm, to
previously implemented approaches for mitigating popularity bias in
recommender systems. We also proposed a new metric to quantify popularity
bias, because the error metric was insufficient for this purpose. On the four
datasets, we were able to demonstrate that the FairBoost algorithm
significantly reduces the popularity bias compared to other algorithms while
keeping the error as low as possible. One exciting area for future study would
be to test this algorithm on other sorts of biases that might exist in the
system such as selection bias, gender bias, and so on. Another thing we can
try is to test the algorithm with different base learners.

## References

  * (1)
  * 30 (1998) 1998\.  MovieLens 100K Dataset.  https://grouplens.org/datasets/movielens/.
  * 29 (2006) 2006\.  Netflix Prize data.  https://www.kaggle.com/netflix-inc/netflix-prize-data.
  * 31 (2008) 2008\.  R2 - Yahoo! Music User Ratings of Songs with Artist, Album, and Genre Meta Information.  https://webscope.sandbox.yahoo.com/.
  * Abdollahpouri et al. (2017) Himan Abdollahpouri, Robin Burke, and Bamshad Mobasher. 2017\.  Controlling Popularity Bias in Learning-to-Rank Recommendation. In _Proceedings of the Eleventh ACM Conference on Recommender Systems_ (Como, Italy) _(RecSys ’17)_. Association for Computing Machinery, New York, NY, USA, 42–46.  https://doi.org/10.1145/3109859.3109912
  * Abdollahpouri et al. (2019) Himan Abdollahpouri, Robin Burke, and Bamshad Mobasher. 2019\.  Managing popularity bias in recommender systems with personalized re-ranking.  _arXiv preprint arXiv:1901.07555_ (2019).
  * Boratto et al. (2021) Ludovico Boratto, Gianni Fenu, and Mirko Marras. 2021\.  Connecting user and item perspectives in popularity debiasing for collaborative recommendation.  _Information Processing & Management_ 58, 1 (2021), 102387\.
  * Calders et al. (2009) Toon Calders, Faisal Kamiran, and Mykola Pechenizkiy. 2009\.  Building classifiers with independency constraints. In _2009 IEEE International Conference on Data Mining Workshops_. IEEE, 13–18.
  * Collins et al. (2018) Andrew Collins, Dominika Tkaczyk, Akiko Aizawa, and Joeran Beel. 2018.  Position bias in recommender systems for digital libraries. In _International Conference on Information_. Springer, 335–344.
  * Ekstrand and Kluver (2021) Michael D Ekstrand and Daniel Kluver. 2021.  Exploring author gender in book rating and recommendation.  _User Modeling and User-Adapted Interaction_ (2021), 1–44.
  * Golbandi et al. (2011) Nadav Golbandi, Yehuda Koren, and Ronny Lempel. 2011\.  Adaptive Bootstrapping of Recommender Systems Using Decision Trees. In _Proceedings of the Fourth ACM International Conference on Web Search and Data Mining_ (Hong Kong, China) _(WSDM ’11)_. Association for Computing Machinery, New York, NY, USA, 595–604.  https://doi.org/10.1145/1935826.1935910
  * Guo et al. (2019) Huifeng Guo, Jinkai Yu, Qing Liu, Ruiming Tang, and Yuzhou Zhang. 2019.  PAL: a position-bias aware learning framework for CTR prediction in live recommender systems. In _Proceedings of the 13th ACM Conference on Recommender Systems_. 452–456.
  * Huang et al. (2020) Jin Huang, Harrie Oosterhuis, Maarten de Rijke, and Herke van Hoof. 2020.  Keeping Dataset Biases out of the Simulation: A Debiased Simulator for Reinforcement Learning based Recommender Systems. In _Fourteenth ACM Conference on Recommender Systems_. 190–199.
  * Iosifidis and Ntoutsi (2018) Vasileios Iosifidis and Eirini Ntoutsi. 2018.  Dealing with bias via data augmentation in supervised learning scenarios.  _Jo Bates Paul D. Clough Robert Jäschke_ 24 (2018).
  * Iosifidis and Ntoutsi (2019) Vasileios Iosifidis and Eirini Ntoutsi. 2019.  Adafair: Cumulative fairness adaptive boosting. In _Proceedings of the 28th ACM International Conference on Information and Knowledge Management_. 781–790.
  * Kamiran and Calders (2009) Faisal Kamiran and Toon Calders. 2009.  Classifying without discriminating. In _2009 2nd International Conference on Computer, Control and Communication_. IEEE, 1–6.
  * Kamiran and Calders (2012) Faisal Kamiran and Toon Calders. 2012.  Data preprocessing techniques for classification without discrimination.  _Knowledge and Information Systems_ 33, 1 (2012), 1–33.
  * Koutsopoulos and Halkidi (2018) Iordanis Koutsopoulos and Maria Halkidi. 2018.  Efficient and Fair Item Coverage in Recommender Systems. In _2018 IEEE 16th Intl Conf on Dependable, Autonomic and Secure Computing, 16th Intl Conf on Pervasive Intelligence and Computing, 4th Intl Conf on Big Data Intelligence and Computing and Cyber Science and Technology Congress (DASC/PiCom/DataCom/CyberSciTech)_. IEEE, 912–918.
  * Mansoury et al. (2020) Masoud Mansoury, Himan Abdollahpouri, Jessie Smith, Arman Dehpanah, Mykola Pechenizkiy, and Bamshad Mobasher. 2020.  Investigating potential factors associated with gender discrimination in collaborative recommender systems.  _arXiv preprint arXiv:2002.07786_ (2020).
  * Obermeyer et al. (2019) Ziad Obermeyer, Brian Powers, Christine Vogeli, and Sendhil Mullainathan. 2019.  Dissecting racial bias in an algorithm used to manage the health of populations.  _Science_ 366, 6464 (2019), 447–453.
  * Ovaisi et al. (2020) Zohreh Ovaisi, Ragib Ahsan, Yifan Zhang, Kathryn Vasilaky, and Elena Zheleva. 2020.  Correcting for selection bias in learning-to-rank systems. In _Proceedings of The Web Conference 2020_. 1863–1873.
  * Parapar and Radlinski (2021) Javier Parapar and Filip Radlinski. 2021.  Diverse User Preference Elicitation with Multi-Armed Bandits. In _Proceedings of the 14th ACM International Conference on Web Search and Data Mining_. 130–138.
  * R. He (2016) J. McAuley R. He. 2016\.  Ups and downs: Modeling the visual evolution of fashion trends with one-class collaborative filtering.
  * Saito et al. (2017) Kuniaki Saito, Yoshitaka Ushiku, and Tatsuya Harada. 2017\.  Asymmetric tri-training for unsupervised domain adaptation. In _International Conference on Machine Learning_. PMLR, 2988–2997.
  * Saito (2020) Yuta Saito. 2020\.  Asymmetric Tri-training for Debiasing Missing-Not-At-Random Explicit Feedback. In _Proceedings of the 43rd International ACM SIGIR Conference on Research and Development in Information Retrieval_. 309–318.
  * San Ramon (2020) Marcy Gordon San Ramon. 2020\.  Ten states sue Google for ’anti-competitive’ online ad sales.  https://brandequity.economictimes.indiatimes.com/news/digital/ten-states-sue-google-for-anti-competitive-online-ad-sales/79771479.
  * Schapire (1999) Robert E Schapire. 1999\.  A brief introduction to boosting. In _Ijcai_ , Vol. 99. Citeseer, 1401–1406.
  * Schnabel et al. (2016) Tobias Schnabel, Adith Swaminathan, Ashudeep Singh, Navin Chandak, and Thorsten Joachims. 2016\.  Recommendations as treatments: Debiasing learning and evaluation. In _international conference on machine learning_. PMLR, 1670–1679.
  * Singh and Joachims (2018) Ashudeep Singh and Thorsten Joachims. 2018.  Fairness of exposure in rankings. In _Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_. 2219–2228.
  * Steck (2010) Harald Steck. 2010\.  Training and testing of recommender systems on data missing not at random. In _Proceedings of the 16th ACM SIGKDD international conference on Knowledge discovery and data mining_. 713–722.
  * Steck (2011) Harald Steck. 2011\.  Item popularity and recommendation accuracy. In _Proceedings of the fifth ACM conference on Recommender systems_. 125–132.
  * Yadav et al. (2019) Himank Yadav, Zhengxiao Du, and Thorsten Joachims. 2019\.  Fair learning-to-rank from implicit feedback.  _arXiv preprint arXiv:1911.08054_ (2019).
  * Zhu et al. (2020) Ziwei Zhu, Jianling Wang, and James Caverlee. 2020.  Measuring and mitigating item under-recommendation bias in personalized ranking systems. In _Proceedings of the 43rd International ACM SIGIR Conference on Research and Development in Information Retrieval_. 449–458.