# Big Bounce cosmology for Palatini $R^{2}$ gravity with a Nieh-Yan term

Flavio Bombacigno<EMAIL_ADDRESS>Physics Department,
“Sapienza” University of Rome, Piazzale Aldo Moro 5, 00185 (Roma), Italy
Giovanni Montani<EMAIL_ADDRESS>ENEA, FSN-FUSPHY-TSM, R.C.
Frascati, Via E. Fermi 45, 00044 Frascati, Italy.
Physics Department, “Sapienza” University of Rome, P.le Aldo Moro 5, 00185
(Roma), Italy

###### Abstract

We analyze the cosmological implementation of Palatini $f(R)$ theories,
constructed with a Nieh-Yan term and solved with respect to the torsion. We
consider the relevant case of the quadratic correction to the Hilbert-Palatini
action in the Ricci scalar, mimicking the Starobinsky model of the metric
formulation. We point out the emergence of peculiar cosmological scenarios,
depending on the sign of such correction, able to reproduce bouncing settings
and to restore the standard Universe dynamics in the late asymptotic limit.
Furthermore, we outline the settling of Little-Rip dynamics, which calls for a
deeper investigation in order to be regularized via matter creation. Finally,
we also show that in our model the Immirzi field is asymptotically frozen in
time, resembling the morphology of Loop Quantum Gravity standard formulation.

##  I Introduction

General Relativity is a very rigorous and self-consistent construction for the
geometrical representation of the gravitational interaction, from the point of
view of the kinematic theory. In this sense, the tensor language arises as the
mathematical implementation of the General Relativity Principle and the
geodesic motion as the natural implication of the Equivalence Principle CQG .
However, the Einsteinian dynamics, associated to the Einstein-Hilbert action
is physically grounded only from the point of view of being the simplest
choice, leading to equations which contain second derivatives of the metric
tensor field. Indeed, simple generalizations of the gravitational action can
be easily constructed by adding other scalar invariants to the Ricci scalar
Nojiri:2010wj , and of particular impact over the last two decades it has been
the so-called $f(R)$ gravity, where the Ricci scalar is replaced by a generic
function of it Sotiriou:2008rp ; Nojiri:2017ncd . Moreover, the success of
this revised gravitational theory must be also attributed to the possibility
of a Brans-Dicke reformulation ST  of the $f(R)$ theory in the so-called
Jordan frame.
Besides the peculiar form of the action, another non-trivial ambiguity
concerning the gravitational interaction is the possibility to consider a
priori the affine connection as an independent variable with respect to the
metric field (Palatini or first order formulation) Pal . Indeed, if the metric
is related to the local causal structure of the space-time, the connection is
in general a different geometric object, responsible for the transport of
tensorial quantities across the space-time manifold. In general relativity,
taking the affine connection as an independent entity leads to a dynamically
equivalent description as it can completely solved in terms of the metric,
i.e. one simply recovers the Levi-Civita connection. However, when the
Palatini formalism is implemented for a $f(R)$ model, although the connection
could be still considered an auxiliary field devoid of a proper dynamics, its
form is affected by the specific form of the function $f(\cdot)$
Sotiriou:2008rp . Especially, it can be seen that the additional contribution
due to the function $f(\cdot)$ allows to recast Palatini $f(R)$ theories into
metric ones endowed with torsion Olmo:2011uz . This means that the affine
connection is equipped with an antisymmetric component depending of the
function $f(\cdot)$, and we deal with a Riemann-Cartan space-time Cartan1 ;
Cartan2 ; Hehl:2013qga ; Hehl:1976kj ; Shapiro:2001rz .
Now, since in Palatini $f(R)$ models torsion naturally emerges, in
Bombacigno:2018tbo  we proposed the idea that for formulating $f(R)$ gravity
in the connection language we have to include torsional contribution already
into the Lagragian. In particular, also in relation with features of Loop
Quantum Gravity (LQG) formalism Ashtekar1 ; Barbero1 ; Holst:1995pc ;
Immirzi:1996di ; Thiemann:2007zz ; Rovelli , we considered in the Lagrangian
density a Nieh Yan term NY1 ; NY2 ; Mercuri:2007ki  with the Immirzi parameter
promoted to be a field Rovelli:1997na ; Calcagni:2009xz ; TorresGomez:2008fj ;
Cianfrani:2009sz ; Bombacigno:2016siz .
This choice allowed to fully solve torsion in terms of the function $f(\cdot)$
and the Immirzi field, reducing the original model to a scalar-tensor theory,
characterized by an interesting phenomenology for the gravitational waves
polarizations Bombacigno:2018tih .
Here, we explore the cosmological implementation of the theory proposed in
Bombacigno:2018tbo  for a flat Friedman Universe, in order to shed light on
the dynamical and physical implications that our revised Palatini $f(R)$
theory can have on the Universe history.
Specifically, we analyze the $f(R)=R+\alpha R^{2}$ model, outlining very
different evolutionary scenarios according the sign of the parameter $\alpha$,
ruling the correction term to the Palatini-Hilibert part of the action
Antoniadis:2018ywb ; Tenkanen:2019jiq ; Edery:2019txq .
For $\alpha>0$ we obtain a modified Friedman dynamics, marked by an effective
gravitational constant depending on the matter content considered, and
restoring the general relativity framework in the asymptotic limit.
Much more intriguing turns out to be the case $\alpha<0$, where we stress the
settling of bouncing cosmologies as a purely classic effect, due to the non
minimal coupling between the Immirzi field and the extended gravity sector of
the theory. In this case our analysis also points out the existence of closed
Universe solutions, even in the absence of spatial curvature, still affected
by the singularity and where general relativity is never restored. Eventually,
it is worth noting that in our model the Immirzi field can be dynamically
relaxed to a constant by the Universe expansion, recovering the standard LQG
perspective.
The paper is structured as follows. In Sec. II we discuss Palatini $f(R)$
gravity and the role played by torsion in this framework; in Sec. III we
briefly recall the main features of our model in the general case, and in Sec.
IV we specialize to the isotropic and homogeneous background. In Sec. V we
consider the implications of the correction term in the Universe evolution.
Finally, in Sec. VI conclusions are drawn.

##  II The role of torsion in Palatini $f(R)$ theories

The action for generic $f(R)$ models in Palatini formulation is given by111We
set $\chi=8\pi G$ and $c=1$.

$S=\frac{1}{2\chi}\int d^{4}x\;\sqrt{-g}\,f(R)+S_{M}[g_{\mu\nu},\psi],$ (1)

where $S_{M}$ represents the matter action and $\psi$ collects globally the
matter fields, which minimally couples to the metric field only. The function
$f(R)$ depends on the Ricci scalar $R$, which according a first order analysis
reads as

$R=g^{\mu\nu}R_{\mu\nu}(\Gamma,\partial\Gamma)=g^{\mu\nu}R\indices{{}^{\rho}_{\mu\rho\nu}}(\Gamma,\partial\Gamma),$
(2)

the Riemann tensor $R\indices{{}^{\mu}_{\nu\rho\sigma}}$ being function of the
affine connection and its derivative, i.e.

$R\indices{{}^{\mu}_{\nu\rho\sigma}}=\partial_{\rho}\Gamma\indices{{}^{\mu}_{\nu\sigma}}-\partial_{\rho}\Gamma\indices{{}^{\mu}_{\nu\rho}}+\Gamma\indices{{}^{\mu}_{\tau\rho}}\Gamma\indices{{}^{\tau}_{\nu\rho}}-\Gamma\indices{{}^{\mu}_{\tau\sigma}}\Gamma\indices{{}^{\tau}_{\nu\rho}}.$
(3)

It is worth remarking that the form of the connection is not established in
the well-know Levi-Civita solution, as in the metric approach, but is
determined properly by the equation of motion obtained from (1). Indeed, if we
assume the affine connection to be symmetric in its lower indices, which a
priori could be not guaranteed, the variation of (1) with respect to the
metric field leads to

$f^{\prime}(R)R_{(\mu\nu)}-\frac{1}{2}f(R)=\chi T_{\mu\nu},$ (4)

with a prime denoting differentiation with respect to the argument and
brackets symmetrization on the indices. The equation for the connection is
given instead by

$\nabla_{\rho}\left(\sqrt{-g}f^{\prime}(R)g^{\mu\nu}\right)=0,$ (5)

where $\nabla_{\mu}$ is the covariant derivative from
$\Gamma\indices{{}^{\rho}_{\mu\nu}}$ and the stress-energy tensor $T_{\mu\nu}$
is defined as

$T_{\mu\nu}\equiv-\frac{2}{\sqrt{-g}}\frac{\delta S_{M}}{\delta g^{\mu\nu}}.$
(6)

Eventually, condition (5) can be still restated as the Levi-Civita definition
for the connection $\Gamma\indices{{}^{\rho}_{\mu\nu}}$, provided we perform a
conformal transformation of the metric $g_{\mu\nu}$, that is

$\tilde{g}_{\mu\nu}\equiv f^{\prime}(R)g_{\mu\nu}.$ (7)

Then, a solution for (5) is given by

$\begin{split}\Gamma\indices{{}^{\rho}_{\mu\nu}}=&\frac{1}{2}\tilde{g}^{\rho\sigma}\left(\partial_{\nu}\tilde{g}_{\mu\sigma}+\partial_{\mu}\tilde{g}_{\nu\sigma}-\partial_{\sigma}\tilde{g}_{\mu\nu}\right)=\\\
=&\frac{1}{2}g^{\rho\sigma}\left(\partial_{\nu}g_{\mu\sigma}+\partial_{\mu}g_{\nu\sigma}-\partial_{\sigma}g_{\mu\nu}\right)+\\\
+&\frac{1}{2}\left(\delta\indices{{}^{\rho}_{\mu}}\partial_{\nu}\ln
f^{\prime}(R)+\delta\indices{{}^{\rho}_{\nu}}\partial_{\mu}\ln
f^{\prime}(R)-g_{\mu\nu}\partial^{\rho}\ln f^{\prime}(R)\right).\end{split}$
(8)

However, if we do not impose at the very beginning any particular symmetries
properties on the form on the connection, solution (8) is not the most general
form the connection we can have Olmo:2011uz . In fact, before variation of the
action be performed, we expect that connection could be endowed with an
antisymmetric component, i.e. torsion tensor might be present:

$T\indices{{}^{\rho}_{\mu\nu}}\equiv\frac{1}{2}\left(\Gamma\indices{{}^{\rho}_{\mu\nu}}-\Gamma\indices{{}^{\rho}_{\nu\mu}}\right)\neq
0,$ (9)

and explicit calculations show that the solution (8) has to be enlarged to
include a vector-like contribution, namely

$\Gamma\indices{{}^{\rho}_{\mu\nu}}\rightarrow\Gamma\indices{{}^{\rho}_{\mu\nu}}-\frac{2}{3}\delta\indices{{}^{\rho}_{\mu}}V_{\nu}.$
(10)

Now, since the symmetric part of the Ricci tensor $R_{(\mu\nu)}$ is not
affected by (10), neglecting or considering any torsional contribution due to
$V_{\nu}$ seems to be dynamically equivalent and we always recover the metric
field equation (4). However, by the inspection of (8) it is clear that also
the function $f(R)$ is responsible for the vector part

$-\frac{2}{3}\delta\indices{{}^{\rho}_{\mu}}V^{f(R)}_{\nu}=\frac{1}{2}\delta\indices{{}^{\rho}_{\mu}}\partial_{\nu}\ln
f^{\prime}(R),$ (11)

so we can imagine to fix $V_{\nu}=-V^{f(R)}_{\nu}$, in order to deal with a
total vanishing vector component in the connection. This choice allows us to
recast (8) in the more suitable form

$\Gamma\indices{{}^{\rho}_{\mu\nu}}=\bar{\Gamma}\indices{{}^{\rho}_{\mu\nu}}+K\indices{{}^{\rho}_{\mu\nu}},$
(12)

where $\bar{\Gamma}\indices{{}^{\rho}_{\mu\nu}}$ is the ordinary Levi-Civita
connection, defined in terms of the metric field $g_{\mu\nu}$, and
$K\indices{{}^{\rho}_{\mu\nu}}$ the so-called contorsion tensor related to
torsion by:

$K\indices{{}^{\rho}_{\mu\nu}}=\frac{1}{2}\left(T\indices{{}^{\rho}_{\mu\nu}}-T\indices{{}_{\mu}^{\rho}{}_{\nu}}-T\indices{{}_{\nu}^{\rho}{}_{\mu}}\right),$
(13)

that from (8) is recognized in

$K_{\rho\mu\nu}=\frac{1}{2}\left(g_{\rho\nu}\partial_{\mu}\ln
f^{\prime}(R)-g_{\mu\nu}\partial_{\rho}\ln f^{\prime}(R)\right).$ (14)

Then, since the contorsion in general still depends on $f^{\prime}(R)$, the
definition (14) (or (8) itself), is well-posed only if one is able to express
$R$ as a function of quantities which do not depend on the connection. With
this regard, if we trace the equation for the metric field (4), we obtain the
structural equation:

$f^{\prime}(R)R-2f(R)=\chi T,$ (15)

which in principle can be solve algebraically for $R=R(T)$, allowing us to
completely determine the connection in terms of the metric field and the
matter source. Now, taking into account (12), the equation for the
gravitational field can be rearranged as:

$\begin{split}\bar{G}_{\mu\nu}(g)=&\frac{\chi}{f^{\prime}(R)}T_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\left(R-\frac{f(R)}{f^{\prime}(R)}\right)+\\\
-&\frac{3}{2f^{\prime}(R)^{2}}\left(\bar{\nabla}_{\mu}f^{\prime}(R)\bar{\nabla}_{\nu}f^{\prime}(R)-\frac{1}{2}g_{\mu\nu}(\bar{\nabla}f^{\prime}(R))^{2}\right)+\\\
+&\frac{1}{f^{\prime}(R)}g_{\mu\nu}\left(\bar{\nabla}_{\mu}\bar{\nabla}_{\nu}-g_{\mu\nu}\bar{\Box}\right)f^{\prime}(R),\end{split}$
(16)

where $\bar{G}_{\mu\nu}$ is the Einstein tensor constructed via the metric
field $g_{\mu\nu}$ only, the d’Alambert operator defined by $\bar{\Box}\equiv
g^{\mu\nu}\bar{\nabla}_{\mu}\bar{\nabla}_{\nu}$, and
$\left(\bar{\nabla}f^{\prime}(R)\right)^{2}\equiv\bar{\nabla}_{\mu}f^{\prime}(R)\bar{\nabla}^{\mu}f^{\prime}(R)$.
In (16), being ultimately $f^{\prime}(R)$ a function of $T$, derivatives of
the stress-energy tensor appear, changing how the matter source generates
space-time curvature. Moreover, when $T=0$, as in vacuum, relation (15) may
admit a set of constant solutions $R=R_{0}^{(i)}$ and in this case equation
(16) simply reduces to GR equation with an effective cosmological constant

$\Lambda_{0}=\frac{1}{2}\left(R_{0}-\frac{f(R_{0})}{f^{\prime}(R_{0})}\right).$
(17)

Therefore, the scenarios offered by Palatini formulation of $f(R)$ theories
depart significantly from standard predictions just in the presence of matter,
where the connection is not an independent variable, but an auxiliary field
affecting the way metric and matter interact. Furthermore, by virtue of (12),
we see that first order $f(R)$ models can equivalently restated as metric
theories endowed with torsion, primarily given by the specific form of the
function $f$. However, if torsion is present at very fundamental level, it
seems reasonable to include torsional contributions already into the
Lagrangian. In this respect, a simple way for achieving that in a LQG-oriented
analysis is offered by $f(R)$ extensions of the well-known Nieh-Yan action,
which represents the starting point of our work.

##  III The Nieh-Yan $f(R)$ model

Let us consider the following extension of the action (1) (see
Bombacigno:2018tbo ), where the so-called Nieh-Yan term is considered in the
presence of a dynamical Immirzi field $\beta(x)$:

$\begin{split}S_{NY}=&\frac{1}{2\chi}\int d^{4}x\;\sqrt{-g}\,f(R)+\\\
+&\frac{1}{4\chi}\int
d^{4}x\;\sqrt{-g}\,\beta(x)\epsilon^{\mu\nu\rho\sigma}\left(g_{\tau\lambda}T\indices{{}^{\tau}_{\mu\nu}}T\indices{{}^{\lambda}_{\rho\sigma}}-R_{\mu\nu\rho\sigma}\right)+\\\
+&S_{M}[g_{\mu\nu,}\psi].\end{split}$ (18)

When $f(R)=R$, action (18) resembles for $\beta$ constant the Nieh-Yan action
usually adopted in Loop Quantum Gravity, and the Immirzi parameter rules a
total divergence that does not affect classically the equations of motion.
However, as stressed in Calcagni:2009xz ; Bombacigno:2016siz ;
Bombacigno:2018tih , if $\beta$ is space-time dependent, it behaves as a
source of torsion and the theory can be reformulated as General Relativity in
the presence of a minimally coupled massless scalar field. Therefore, if we
keep $f(R)\neq R$ generic, we expect that both types of torsion could generate
a dynamic larger framework with respect to (16), by virtue of the non trivial
coupling between the Immirzi field and the gravitational degrees of freedom
(d.o.f).
Then, following the analysis of Bombacigno:2018tbo  that here we are widening
for including matter as well, with a bit of algebra action (18) can be
rewritten in the more convenient scalar-tensor form:

$\begin{split}S_{NY}=&\frac{1}{2\chi}\int
d^{4}x\;\sqrt{-g}\left(\phi\bar{R}-g^{\mu\nu}\Pi_{\mu\nu}(\phi,\beta)-V(\phi)\right)\\\
&+S_{M}[g_{\mu\nu,}\psi],\end{split}$ (19)

where $\phi\equiv f^{\prime}(R)$ and $\bar{R}$ represents the metric Ricci
scalar, depending only on metric variables. In particular, we introduced the
quantities

$\displaystyle\Pi_{\mu\nu}(\phi,\beta)=\frac{3}{2\phi}\left(\bar{\nabla}_{\mu}\beta\bar{\nabla}_{\nu}\beta-\bar{\nabla}_{\mu}\phi\bar{\nabla}_{\nu}\phi\right)$
(20) $\displaystyle V(\phi)\equiv\phi R(\phi)-f(R(\phi)).$ (21)

By analogy with (14), dealing at the effective level with (19) means having a
contorsion tensor given by

$K_{\rho\mu\nu}=\frac{1}{2\phi}\left(g_{\rho\nu}\bar{\nabla}_{\mu}\phi-
g_{\mu\nu}\bar{\nabla}_{\rho}\phi\right)-\frac{1}{2\phi}\epsilon_{\rho\mu\nu\sigma}\bar{\nabla}^{\sigma}\beta,$
(22)

which for $\beta$ constant boils down to (14).
Then, varying (19) with respect to $g_{\mu\nu}$ carries out:

$\begin{split}\bar{G}_{\mu\nu}=&\frac{\chi}{\phi}T_{\mu\nu}-\frac{1}{2\phi}g_{\mu\nu}V(\phi)+\\\
+&\frac{1}{\phi}\left(\Pi_{\mu\nu}(\phi,\beta)-\frac{1}{2}g_{\mu\nu}\Pi\indices{{}^{\rho}_{\rho}}(\phi,\beta)\right)+\\\
+&\frac{1}{\phi}\left(\bar{\nabla}_{\mu}\bar{\nabla}_{\nu}-g_{\mu\nu}\bar{\Box}\right)\phi,\end{split}$
(23)

while the equations for $\phi$ and $\beta$ are given by, respectively:

$\bar{R}=-\frac{3}{2\phi^{2}}\left(\bar{\nabla}_{\mu}\beta\bar{\nabla}^{\mu}\beta+\bar{\nabla}_{\mu}\phi\bar{\nabla}^{\mu}\phi\right)+\frac{3\bar{\Box}\phi}{\phi}+V^{\prime}(\phi)$
(24)

and

$\bar{\Box}\beta(x)=\frac{\bar{\nabla}_{\mu}\beta(x)\bar{\nabla}^{\mu}\phi}{\phi}.$
(25)

Substituting in (24) the trace of (23), we get the modified structural
equation (see (15) for comparison):

$2\;V(\phi)-\phi\;V^{\prime}(\phi)=\chi
T-\frac{3\bar{\nabla}_{\mu}\beta\bar{\nabla}^{\mu}\beta}{\phi}.$ (26)

that once we chose the specific $f(R)$ model, allows us to solve for
$\phi=\phi((\bar{\nabla}\beta)^{2},T)$. We point out that still in vacuum
relation (26) admits a larger set of solutions for $\phi$, which is not
compelled to relax to a constant value as in original formulation. The reason
for this is the non trivial coupling between the Immirzi field and the
enlarged gravitational sector, which makes the vacuum configuration never
actually devoid of matter content. In particular, the Immirzi term in (26)
only partially resembles the contribution to the structural equation of scalar
fields in ordinary Palatini $f(R)$ gravity, where we would expect a standard
kinetic term deprived of the coupling with $\phi$. Furthermore, we note that
the requirement of recovering to some extent a proper vacuum state, as it is
described by (15) for $T=0$, raises the issue about the relaxation of the
Immirzi field on a constant configuration, able to reproduce standard LQG
theory as well. In this respect, these problems can be properly addressed in
cosmology, where we can ask for the asymptotically freezing of the Immirzi
field during the expansion of the Universe. Of course, such a mechanism does
not accounts for local fluctuations $\delta\beta(x)$ (see Bombacigno:2018tbo ;
Bombacigno:2018tih ), but it can be considered responsible for the evolution
of the background value $\beta_{B}$, where $\beta=\beta_{B}+\delta\beta$,
which we may demand to match with independent LQG estimates Ghosh:2004wq ;
Ghosh:2011fc .

##  IV Modified Friedmann equation

A simple model for analyzing the role played by $\beta_{B}$ is represented by
the homogeneous and isotropic Universe, described by the Friedman-Robertson-
Walker (FRW) line element222We set the lapse function $N(t)=1$.

$ds^{2}=-dt^{2}+a^{2}(t)\left(\frac{dr^{2}}{1-kr^{2}}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2})\right),$ (27)

the scale factor $a(t)$ being the only metric dynamical degree of freedom and
$k$ the curvature of space. Within such a framework, the background value for
the Immirzi field can be considered function of the cosmological time $t$
only, namely $\beta_{B}=\beta_{B}(t)$. Now, be $T_{\mu\nu}$ the stress-energy
tensor for a perfect fluid, i.e.:

$T_{\mu\nu}=\left(\rho+P\right)u_{\mu}u_{\nu}+g_{\mu\nu}P,$ (28)

where $\rho$ and $P$ are the energy density and the pressure, respectively,
and $u_{\mu}=(-1,0,0,0)$. Then, in the presence of the energy density $\rho$
the Friedman equation stemming from (23) can be rearranged as

$H^{2}=\left(\frac{\dot{a}}{a}\right)^{2}=\frac{\chi}{3\phi}\left(\rho_{eff}+\rho\right)-\frac{k}{a^{2}},$
(29)

with dot denoting time derivative, while combining the equation for the $ii$
component with (29) we get the acceleration equation

$\frac{\ddot{a}}{a}=-\frac{\chi}{6\phi}\left(\rho+\rho_{eff}+3(P+P_{eff})\right),$
(30)

where we introduced the effective energy density and pressure given by,
respectively:

$\rho_{eff}\equiv\frac{1}{\chi}\left[\frac{3}{4\phi}\left(\dot{\beta}_{B}^{2}-\dot{\phi}^{2}\right)+\frac{1}{2}V(\phi)-3\frac{\dot{a}}{a}\dot{\phi}\right]$
(31)

and

$P_{eff}\equiv\frac{1}{\chi}\left[\frac{3}{4\phi}\left(\dot{\beta}_{B}^{2}-\dot{\phi}^{2}\right)-\frac{1}{2}V(\phi)+\ddot{\phi}+2\frac{\dot{a}}{a}\dot{\phi}\right].$
(32)

Deriving equation (29) with respect to time and plugging (30) in it, we can
obtain the relation

$\begin{split}&\dot{\rho}+3\left(\frac{\dot{a}}{a}\right)\left(\rho+P\right)+\\\
&+\dot{\rho}_{eff}-\frac{\dot{\phi}}{\phi}\left(\rho_{eff}+\rho\right)+3\left(\frac{\dot{a}}{a}\right)\left(\rho_{eff}+P_{eff}\right)=0.\end{split}$
(33)

Then, in order the standard continuity equation be preserved, that is:

$\dot{\rho}+3\frac{\dot{a}}{a}\left(\rho+P\right)=0,$ (34)

the following condition has to be fulfilled

$\dot{\rho}_{eff}-\frac{\dot{\phi}}{\phi}\left(\rho+\rho_{eff}\right)+3\frac{\dot{a}}{a}\left(\rho_{eff}+P_{eff}\right)=0.$
(35)

By virtue of (31) and (32) this relation can be rewritten as

$\frac{\dot{\phi}}{2}\left(V^{\prime}(\phi)-\bar{R}+\frac{3\dot{\phi}^{2}}{2\phi}-\frac{3\dot{\beta}_{B}^{2}}{2\phi^{2}}+\frac{3\bar{\Box}\phi}{\phi}\right)-\frac{3\dot{\beta}_{B}\bar{\Box}\beta_{B}}{2\phi}=0,$
(36)

where we used the expressions of the Ricci scalar and the d’Alambert operator
for the background (27). Hence, using (25), relation (36) takes the form

$\frac{\dot{\phi}}{2}\left(V^{\prime}(\phi)-\bar{R}+\frac{3\dot{\phi}^{2}}{2\phi}+\frac{3\dot{\beta}_{B}^{2}}{2\phi^{2}}+\frac{3\bar{\Box}\phi}{\phi}\right)=0,$
(37)

which is identically satisfied given (24).
Therefore, if the equation of state $P=w\rho$ holds, where $w$ is the
polytropic index, from (34) the standard solution can be obtained, namely

$\rho(a)=\frac{\mu^{2}}{a^{3(w+1)}},$ (38)

with $\mu$ a constant.
Furthermore, we note that the equation (25) for $\beta_{B}$ can be actually
solved analitically for $\dot{\beta}_{B}$. Indeed, for a FRW background, (25)
simply reads as:

$\ddot{\beta}_{B}(t)+\left(3\frac{\dot{a}}{a}-\frac{\dot{\phi}}{\phi}\right)\dot{\beta}_{B}=0,$
(39)

whose solution is given by:

$\dot{\beta}_{B}(t)=C_{0}\frac{\phi(t)}{a^{3}(t)},$ (40)

where we defined the integration constant
$C_{0}\equiv\frac{\dot{\beta}_{B}(t_{0})a(t_{0})^{3}}{\phi(t_{0})}$ for a
fiducial instant $t_{0}$. Thus, inserting (38) and (40) in (26) yields:

$2V(\phi)-\phi
V^{\prime}(\phi)=\chi\mu^{2}\frac{3w-1}{a^{3(w+1)}}+\frac{3C_{0}^{2}}{a^{6}}\phi,$
(41)

that, once a peculiar $f(R)$ model has been fixed, allows us to express $\phi$
as a function of the scale factor $a(t)$, by virtue of (21). Moreover, given
(40), this implies in turn that $\dot{\beta}_{B}(t)$ depends on time only by
means of the scale factor. Therefore, the requirement that the Immirzi field
relaxes on a constant, can be equivalently restated as

$\lim_{a\rightarrow+\infty}\dot{\beta}_{B}=0.$ (42)

Now, taking into account (40), the Friedman equation can be reformulated as:

$\begin{split}H^{2}=&\frac{\chi}{3\phi}\frac{\mu^{2}}{a^{3(w+1)}}+\frac{C_{0}^{2}}{4a^{6}}+\frac{V(\phi)}{6\,\phi}-\frac{1}{4}\frac{\dot{\phi}^{2}}{\phi^{2}}-\frac{\dot{a}}{a}\frac{\dot{\phi}}{\phi}-\frac{k}{a^{2}}.\end{split}$
(43)

We note that, since now $\phi$ has to be understood by means of (41) as a
function of the scale factor $a$, the terms in the R.H.S. of (43) depending on
the time derivative of $\phi$ always give rise to terms proportional to
$H^{2}$, regardless the $f(R)$ model considered. Therefore, it is possible to
rearrange (43) in the following way:

$H^{2}=\frac{\left(\frac{\chi\rho(a)}{3\phi}+\frac{C_{0}^{2}}{4a^{6}}+\frac{V(\phi)}{6\,\phi}-\frac{k}{a^{2}}\right)_{\phi=\phi(a)}}{F_{1}(a)^{2}},$
(44)

where $F_{1}(a)$ is a function of the scale factor that has the general form

$F_{1}(a)=\left(1+\frac{a}{2}\frac{d\ln\phi(a)}{da}\right),$ (45)

and the term depending on the time derivative of the Immirzi field appears, by
virtue of (40), as a sort of scalar field energy density.
We note that the behaviour of $F_{1}(a)$ and $\phi(a)$ can remarkably affect
the evolution of the scale factor. Indeed, in the presence of any polos and
zeros for $F_{1}(a)$, the function $H$ can vanish or diverge, giving rise to
peculiar cosmological scenarios. Similarly, by virtue of the coupling with the
energy density content and the potential term, also $\phi(a)$ can be in
principle responsible for analogous frameworks.

##  V The model $f(R)=R+\alpha R^{2}$

In the following, we will restrict our attention on a specific Lagrangian,
including a correction term quadratic in the total Ricci scalar $R$, i.e.:

$f(R)=R+\alpha R^{2}.$ (46)

It is worth noting that with respect to the metric approach (the well-
established Starobinsky model Starobinsky ; Hawking:2000bb ; Kehagias:2013mya
; Whitt:1984pd ), in Palatini formulation there are no issues concerning
possible instabilities of the solution Olmo:2011uz ; Sotiriou:2006sf , being
that ultimately due to the non dynamical nature of the field $\phi$. For this
reason, the real parameter $\alpha$ is not compelled a priori to be positive,
and also negative values represent a suitable choice.
Then, when the model (46) is taken into account, the potential $V(\phi)$ can
be easily found, that is

$V(\phi)=\frac{1}{\alpha}\left(\frac{\phi-1}{2}\right)^{2},$ (47)

which inserted in (41) gives us:

$\phi=\frac{a^{6}\,f(a)}{a^{6}+6\alpha C^{2}_{0}},$ (48)

being $f(a)$ a function of the energy density, i.e.

$f(a)=1-2\alpha\chi(3w-1)\rho(a).$ (49)

Eventually, setting $k=0$, by means of (47)–(48) the Friedman equation can be
rearrange as

$H^{2}=\frac{(a^{6}+6\alpha
C_{0}^{2})\left(4\chi\rho+\frac{3C_{0}^{2}f(a)}{a^{6}+6\alpha
C_{0}^{2}}+\frac{2\alpha\left(\chi(3w-1)a^{6}\rho+3C_{0}^{2}\right)^{2}}{(a^{6}+6\alpha
C_{0}^{2})^{2}}\right)}{12a^{6}f(a)\left(\frac{a^{6}+24\alpha
C_{0}^{2}}{a^{6}+6\alpha C_{0}^{2}}+\frac{af^{\prime}(a)}{2f(a)}\right)^{2}}.$
(50)

By first inspection of (50), we see that according the sign of $\alpha$ the
parameter $C_{0}^{2}$, related to the Immirzi energy density, is crucial in
determining the critical points of the Friedman equation. In particular, with
the aim of investigating the possible emergence of bouncing cosmologies ruled
by the Immirzi field Olmo:2008nf ; Barragan:2009sq ; Koivisto:2010jj , it can
be instructive to consider the vacuum case $\rho=0$, where $f(a)=1$ and (50)
takes a very simple form. More complex examples, even if still feasible for
analytic studies, are represented both by the cosmological constant case,
where $\rho$ is constant and the term $f^{\prime}(a)$ in (50) vanishes, and by
the radiation one, where the trace of $T_{\mu\nu}$ is zero and $f(a)=1$ as in
vacuum.

###  V.1 The vacuum case

The vacuum model constitutes a very useful tool for studying the effects, on
the space-time structure, of the Immirzi coupling to gravitational d.o.f. . In
this case relation (48) is simply

$\phi(a)=\frac{a^{6}}{a^{6}+6\alpha C_{0}^{2}}$ (51)

and (50) boils down to

$H^{2}=\frac{C_{0}^{2}}{4a^{6}}\frac{(a^{6}+6\alpha C_{0}^{2})(a^{6}+12\alpha
C_{0}^{2})}{(a^{6}+24\alpha C_{0}^{2})^{2}}.$ (52)

When $\alpha>0$, equation (51) does not exhibit critical points and it can be
recast into the form

$H^{2}=\frac{\chi_{\beta}(a)}{3}\rho_{\beta},$ (53)

which represents the Friedman equation for the scalar field energy density
$\rho_{\beta}\equiv\frac{3C_{0}^{2}}{4\chi a^{6}}$. It is characterized by an
effective gravitational constant

$\chi_{\beta}(a)\equiv\frac{(a^{6}+6\alpha C_{0}^{2})(a^{6}+12\alpha
C_{0}^{2})}{(a^{6}+24\alpha C_{0}^{2})^{2}}\,\chi,$ (54)

and General Relativity is recovered for $a\rightarrow+\infty$, where
$\chi_{\beta}\rightarrow\chi$ and $\phi\rightarrow 1$ in agreement333The value
$\phi=1$ corresponds to $f^{\prime}(R)=1$. with (51).
When instead $\alpha<0$, the presence of $C_{0}^{2}$ affects drastically (51)
. Indeed, in order the condition $H^{2}\geq 0$ be preserved, the scale factor
$a$ cannot assume arbitrary values in $\mathbb{R}^{+}$, but is constrained
into domains

$\begin{split}&\mathcal{D}_{1}:\;a\in[0,(-6\alpha C_{0}^{2})^{1/6}],\\\
&\mathcal{D}_{2}:\;a\in[(-12\alpha C_{0}^{2})^{1/6},+\infty).\end{split}$ (55)

Therefore, we deal with two disconnected branches, denoting two different kind
of Universe. In particular, the region $\mathcal{D}_{1}$ defines a closed
Universe, even for $k=0$, bounded by the turning point $a_{T}=(-6\alpha
C_{0}^{2})^{1/6}$, where the General Relativity limit is never reached
($\phi=0$ for $a=0$ and $\phi\rightarrow-\infty$ for $a\rightarrow a_{T}$) and
it can be disregarded since unphysical.
Instead, the brach $\mathcal{D}_{2}$ is endowed with the critical point
$a_{B}=(-12\alpha C_{0}^{2})^{1/6}$ where $H^{2}=0$ and a bounce occurs,
driven by the Immirzi energy density. That can be further proved by evaluating
(30) at the bounce, where it can be recast into the form

$\frac{\ddot{a}}{a}=\frac{-\frac{\chi(1+3w)}{6\phi}\rho-\frac{C_{0}^{2}}{2a^{6}}+\frac{V(\phi)}{\phi}-\frac{1}{2}H^{2}F_{2}(a)}{F_{1}(a)},$
(56)

with $F_{2}(a)$ given by

$F_{2}(a)\equiv a^{2}\frac{d^{2}}{da^{2}}\ln\phi+a\frac{d}{da}\ln\phi.$ (57)

Now, since $\phi$ is not singular at $a=a_{B}$, at the bounce ($H^{2}=0$)
relation (56) simply gives

$\left.\frac{\ddot{a}}{a}\right|_{a=a_{B}}=-\frac{1}{32\alpha},$ (58)

which is positive for $\alpha<0$.
Remarkably, in this case (53) can be put in the LQC-like form

$H^{2}=\frac{\chi}{3}\rho_{\beta}\left(1-\frac{\rho_{\beta}}{\rho_{crit}^{vac}}\right),$
(59)

where we introduced the critical density

$\rho_{crit}^{vac}\equiv\frac{(a^{6}+24\alpha C_{0}^{2})^{2}}{8\alpha\chi
a^{6}(5a^{6}+84\alpha C_{0}^{2})}.$ (60)

With respect to Olmo:2008nf ; Barragan:2009sq , where analogous results were
discussed, we stress that in our case we are able to reproduce bouncing
cosmology for (46) also in the presence of stiff-like matter ($w=1$) (properly
mimicked by the Immirzi field contribute) when $\alpha<0$. Moreover, requiring
that the bounce occurs for Planckian energy density, allows to set the order
of magnitude of the parameter $\alpha$. Indeed, if at the bounce

$\rho_{\beta}=\rho_{crit}^{vac}(a_{B})=-\frac{1}{16\alpha\chi}\sim\rho_{Planck},$
(61)

where444For the sake of clarity here we show explicitly the speed of light
$c$, that in the rest of the work we set to unity. $\rho_{Planck}=c^{7}/\hbar
G^{2}$, then it follows that $\alpha$ can be estimated by

$|\alpha|\sim\frac{\hbar G}{128\pi c^{3}}.$ (62)

We note that the branch $\mathcal{D}_{2}$ is marked by another peculiar point,
namely $a=a_{R}=(-24\alpha C_{0}^{2})^{1/6}$ where the function $H^{2}$
diverges, and we have a vanishing Hubble radius for a finite scale factor
Odintsov:2018uaw ; Astashenok:2012tv . We expect that this type of singularity
could be regularized taking into account the gravitational particle creation
Montani:2001fp ; Dimopoulos:2018kgl ; Ford:1986sy , related to the presence of
a cosmological horizon, or the non-equilibrium nature of the involved
thermodynamic processes Disconzi:2014oda ; Venanzi:2016pds ; Israel:1976tn ,
like bulk viscosity effects Barrow:1988yc ; Belinskii1 ; Belinskii2 . In
particular, particle creation can be described by means of additional terms in
the Friedman equation, able to stabilize the singular behaviour of the Hubble
parameter Contreras:2018two . Therefore, we hypothesize that the Universe
might evolve smoothly through the critical point $a_{R}$, reaching
asymptotically the General Relativity regime (53), where we also require that
the Immirzi field relaxes on a constant value. That can be easily checked
combining (48) and (40), namely

$\lim_{a\rightarrow+\infty}\dot{\beta}(a)=\lim_{a\rightarrow+\infty}\frac{a^{3}(1-2\alpha\chi(3w-1)\rho)}{a^{6}+6\alpha
C_{0}^{2}}C_{0}^{2}=0,$ (63)

which, providing $w\geq-2$, holds irrespective of the specific energy density
content $\rho$.

###  V.2 The cosmological constant case

For $w=-1$ the energy density does not depend on the scale factor and we can
formally set $\rho=\Lambda/\chi$, where $\Lambda$ is a cosmological constant
term. Then, relation (48) reads as

$\phi(a)=\frac{a^{6}(1+8\alpha\Lambda)}{a^{6}+6\alpha C_{0}^{2}},$ (64)

and (50) takes the form

$H^{2}=\frac{(a^{6}+6\alpha C_{0}^{2})(4\Lambda
a^{12}+3C_{0}^{2}a^{6}+36\alpha C_{0}^{4})}{12a^{6}(a^{6}+24\alpha
C_{0}^{2})^{2}}.$ (65)

By close analogy with (53), for $\alpha>0$ it can be simply recast as

$H^{2}=\frac{\Lambda_{\beta}(a)}{3}+\frac{\chi_{\beta}(a)}{3}\rho_{\beta},$
(66)

with the effective cosmological constant $\Lambda_{\beta}$ given by

$\Lambda_{\beta}(a)\equiv\frac{a^{6}(a^{6}+6\alpha C_{0}^{2})}{(a^{6}+24\alpha
C_{0}^{2})^{2}}\;\Lambda,$ (67)

and for $a\rightarrow+\infty$, the dynamical term $\Lambda_{\beta}$
asymptotically reaches the constant value $\Lambda$ and the standard de Sitter
phase is recovered. On the other hand, near the singularity the $\Lambda$ term
is negligible, i.e. $\Lambda_{\beta}\rightarrow 0$ and the Immirzi field is
the leading contribution to the dynamics.
If instead $\alpha<0$, the requirement of having a positive Hubble parameter
compels once again the scale factor into specific regions of $\mathbb{R}^{+}$.
Specifically, assuming the value of $\Lambda$ fixed, as for instance by
current data Aghanim:2018eyx , it is possible to distinguish two separate
cases, labelled by the size of $\alpha$ with respect to $\Lambda$, i.e.:

$\begin{split}8\Lambda\alpha<-1&\Rightarrow\begin{cases}&\mathcal{D}_{1}^{\Lambda}\;a\in[0,a_{\Lambda}],\\\
&\mathcal{D}_{2}^{\Lambda}\;a\in[(-6\alpha
C_{0}^{2})^{1/6},+\infty);\end{cases}\\\
-1<8\Lambda\alpha<0&\Rightarrow\begin{cases}&\mathcal{D}_{3}^{\Lambda}\;a\in[0,(-6\alpha
C_{0}^{2})^{1/6}],\\\
&\mathcal{D}_{4}^{\Lambda}\;a\in[a_{\Lambda},+\infty);\end{cases}\end{split}$
(68)

where
$a_{\Lambda}=\left(-\frac{3C_{0}^{2}}{8\Lambda}\left(1-\sqrt{1-64\alpha\Lambda}\right)\right)^{1/6}$.
Analogously to the vacuum case, the domains $\mathcal{D}_{1,3}^{\Lambda}$
always designate closed Universes, which do not admit General Relativity as
limit, and they can be overlooked. Instead, branches
$\mathcal{D}_{2,4}^{\Lambda}$ describe bouncing cosmologies, with the Big
Bounce point critically depending on the value of $\alpha$. Especially, when
$-1<8\Lambda\alpha<0$ holds, the bounce takes place for values corresponding
to the turning point $a_{T}$ of the vacuum case (see (55)), while if
$8\Lambda\alpha<-1$ the Big Bounce point is determined by $a_{\Lambda}$ and
also the cosmological constant term is involved in fixing its value. Moreover,
it is easy to see that in both cases the bounce occurs for scale factor values
lower than in vacuum, being $a_{\Lambda}<a_{B}=(-12\alpha C_{0}^{2})^{1/6}$
always satisfied for $\alpha<0$. However, if we assume $\alpha$ fixed by (62)
and $\Lambda$ reproducing the current dark energy phase ($\Lambda\sim
10^{-18}l_{P}^{-2}$), then we see that the condition $8\alpha\Lambda<-1$
cannot be satisfied, and $\mathcal{D}_{4}^{\Lambda}$ is the only valid branch.
Finally, the critical point $a_{R}=(-24\alpha C_{0}^{2})^{1/6}$, where $H$
diverges, is not removed since for negative values of $\alpha$ it is always
contained in the $\mathcal{D}_{2,4}^{\Lambda}$ domains.

###  V.3 The radiation case

When $w=1/3$ the trace of the stress-energy tensor vanishes and relation (51)
is unaltered, whereas the Friedman equation (52) is slightly modified and
reads as:

$H^{2}=\frac{C_{0}^{2}}{4a^{6}}\frac{(a^{6}+6\alpha C_{0}^{2})(a^{6}+12\alpha
C_{0}^{2}+\frac{4\chi\mu_{R}^{2}}{3C_{0}^{2}}\frac{(a^{6}+6\alpha
C_{0}^{2})^{2}}{a^{4}})}{(a^{6}+24\alpha C_{0}^{2})^{2}}.$ (69)

Following (66), when $\alpha>0$ it can be rewritten as

$H^{2}=\frac{\chi_{R}(a)}{3}\rho_{R}+\frac{\chi_{\beta}(a)}{3}\rho_{\beta},$
(70)

where with analogy with (67) we defined the effective gravitational coupling

$\chi_{R}(a)=\frac{(a^{6}+6\alpha C_{0}^{2})^{3}}{a^{6}(a^{6}+24\alpha
C_{0}^{2})^{2}}\;\chi.$ (71)

In particular, by virtue of (54)-(71), we see that near the singularity the
Immirzi energy density is negligible and the Friedman equation behaves like
$H^{2}\sim a^{-10}$, corresponding to an effective superluminal index $w=7/3$.
We observe that such results are quite common in ekpyrotic theories (see
Lehners:2008vx  and references therein for an introduction), where it is in
general requested $w\gg 1$ in order to solve the fine tuning issues of
standard cosmological model.
When $\alpha<0$, it can be demonstrated with bit long calculations that the
effect of the radiation energy density is twofold: It endows the Hubble
function of an additional zero $a_{B_{1}}$ with respect to the vacuum case and
displaces the critical point $a_{B}=(-12\alpha C_{0}^{2})^{1/6}$ in a new root
$a_{B_{2}}$. Even if such two new zeros cannot be analytically evaluated, they
may be still algebraically estimated by

$\begin{split}&a_{B_{1}}\in\left(0;(-6\alpha C_{0}^{2})^{1/6}\right)\\\
&a_{B_{2}}\in\left((-6\alpha C_{0}^{2})^{1/6};(-12\alpha
C_{0}^{2})^{1/6}\right).\end{split}$ (72)

Accordingly, the regions where relation $H^{2}\geq 0$ is valid are changed
into the following new domains:

$\begin{split}&\mathcal{D}_{1}^{R}:\;a\in[a_{B_{1}},(-6\alpha
C_{0}^{2})^{1/6}],\\\
&\mathcal{D}_{2}^{R}:\;a\in[a_{B_{2}},+\infty),\end{split}$ (73)

and we see that the unphysical branch $\mathcal{D}_{1}^{R}$ is now turned in a
cyclic Universe equipped with a proper bounce point. Concerning instead
$\mathcal{D}_{2}^{R}$, we note that the Big Bounce is shifted to lower values,
as for the cosmological constant case, whereas the point of divergence $a_{R}$
is unaffected.

##  VI Concluding remarks

The analysis above provided the cosmological implementation of the Palatini
$f(R)$ model discussed in Bombacigno:2018tbo , where a Nieh-Yan term was
included in the presence of an Immirzi field. The peculiarity of that approach
was the possibility to completely solve torsion in terms of the Immirzi and
gravitational d.o.f., so obtaining a non-minimally coupled scalar-tensor
theory.
In particular, we considered a specific class of $f(R)$ theories, mimicking
the well-established Starobinsky model in metric formalism, by means of a
quadratic correction to the Palatini-Hilbert action. In this regard, we
clearly distinguished two different cosmological scenarios, depending on the
sign of such a correction. Indeed, for $\alpha>0$ the analysis outlined a
slightly modified Friedman dynamics, approaching in the asymptotic limit the
standard description of general relativity and characterized by effective
gravitational couplings, according the type of energy density considered. For
$\alpha<0$ instead, we pointed out the existence of radically different
solutions, consisting in closed and bouncing Universes, respectively.
Especially, the former were obtained even for vanishing spatial curvature, but
in general they turned out to be still singular, and they were ruled out
because of the absence of the general relativity limit. Concerning the latter,
we were able to identify in the non minimal coupling of the Immirzi field with
the gravitational d.o.f. the cause of the classical removal of the initial
singularity. In this respect, it is worth stressing that when a radiation
energy density was taken into account, the combined effect of the Immirzi and
radiation field was to introduce a further bouncing point in the closed
solution, resulting in a cyclic model. We mention that this kind of solutions,
even if a priori disregarded, they could represent a Planckian state of the
Universe, from which the bouncing branch could originate as the result of a
quantum tunneling phenomenon.
Moreover, we pointed out that in general the reliable classical solutions are
always endowed with critical points associated to little rip dynamics, where
the Friedman equation diverges for specific values of the scale factor. Of
course, they must be regularized by reducing to a finite value the expansion
rate, as effect of matter creation, as well as non equilibrium thermodynamics
implications, mainly associated to bulk viscosity effects.
Eventually, we shown that the Immirzi field can asymptotically relaxed to a
constant by the expansion of the Universe, newly recovering the LQG vision of
an Immirzi parameter.

## References

  * (1) F. Cianfrani, O.M. Lecian, M. Lulli, G. Montani, Canonical quantum gravity (World Scientific Pub Co Inc., Singapore, 2014)
  * (2) S. Nojiri and S. D. Odintsov, Phys. Rept. 505, 59 (2011)
  * (3) T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. 82 (2010) 451
  * (4) S. Nojiri, S. D. Odintsov and V. K. Oikonomou, Phys. Rept. 692 (2017) 1
  * (5) Y. Fujii and K. Maeda, The Scalar-tensor theory of Gravitation (Cambridge University Press, Cambridge, 2009)
  * (6) A. Palatini, Rend. Circ. Mat. Palermo 43, 203 (1919) [English translation by R. Hojman and C. Mukku, in Cosmology and Gravitation, edited by P.G. Bergmann and V. De Sabbata (Plenum Press, New York, 1980)]
  * (7) G. J. Olmo, Int. J. Mod. Phys. D 20, 413 (2011)
  * (8) É. Cartan. Sur une généralisation de la notion de courbure de Riemann et les espaces à torsion, Cosmology and Gravitation: Spin, Torsion, Rotation, and Supergravity, Nato Science Series B, edited by P.G. Bergmann and V. De Sabbata (Springer, New York, 1980), p. 489.
  * (9) É. Cartan. Ann. Sci. École Norm. Sup. 40, 325 (1923); Ann. Sci. Éc. Norm. Sup. 41, 1 (1924); Ann. Sci. École Norm. Sup. 42, 17–88 (1925).
  * (10) F. W. Hehl, P. Von Der Heyde, G. D. Kerlick and J. M. Nester, Rev. Mod. Phys. 48, 393, (1976) .
  * (11) I. L. Shapiro, Phys. Rept. 357 (2002) 113
  * (12) F. W. Hehl, Y. N. Obukhov and D. Puetzfeld, Phys. Lett. A 377, 1775 (2013)
  * (13) F. Bombacigno and G. Montani, Phys. Rev. D 97, 124066 (2019)
  * (14) A. Ashtekar, Phys. Rev. Lett 57 2244 (1986)
  * (15) J.F. Barbero G., Phys. Rev. Lett 51 5507 (1995)
  * (16) S. Holst, Phys. Rev. D 53 5966 (1996)
  * (17) G. Immirzi, Class. Quant. Grav. 14, L177 (1997)
  * (18) T. Thiemann, Modern canonical quantum general relativity (Cambridge University Press, Cambridge, 2007)
  * (19) C. Rovelli, Quantum Gravity (Cambridge University Press, Cambridge, 2004)
  * (20) H.T. Nieh and M.L. Yan, J. Math. Phys. 23, 373-374 (1982);
  * (21) H.T. Nieh, Int. J. Mod. Phys. A22, 5237-5244 (2007).
  * (22) S. Mercuri, Phys. Rev. D 77, 024036 (2008)
  * (23) C. Rovelli and T. Thiemann, Phys. Rev. D 57, 1009 (1998)
  * (24) G. Calcagni and S. Mercuri, Phys. Rev. D 79, 084004 (2009)
  * (25) A. Torres-Gomez and K. Krasnov, Phys. Rev. D 79, 104014 (2009)
  * (26) F. Cianfrani and G. Montani, Phys. Rev. D 80, 084040 (2009)
  * (27) F. Bombacigno, F. Cianfrani and G. Montani, Phys. Rev. D 94 (2016) no.6, 064021
  * (28) F. Bombacigno and G. Montani, Phys. Rev. D 99, 064016 (2019)
  * (29) I. Antoniadis, A. Karam, A. Lykkas and K. Tamvakis, JCAP 1811, no. 11, 028 (2018)
  * (30) T. Tenkanen, arXiv:1901.01794 [astro-ph.CO].
  * (31) A. Edery and Y. Nakayama, arXiv:1902.07876 [hep-th].
  * (32) A. Ghosh and P. Mitra, Phys. Lett. B 616, 114 (2005)
  * (33) A. Ghosh and A. Perez, Phys. Rev. Lett. 107, 241301 (2011) Erratum: [Phys. Rev. Lett. 108, 169901 (2012)]
  * (34) A.A. Starobinsky, Phys. Lett. B, 91, 99 (1980)
  * (35) S. W. Hawking, T. Hertog and H. S. Reall, Phys. Rev. D 63, 083504 (2001)
  * (36) A. Kehagias, A. Moradinezhad Dizgah and A. Riotto, Phys. Rev. D 89, no. 4, 043527 (2014)
  * (37) B. Whitt, Phys. Lett. 145B, 176 (1984).
  * (38) T. P. Sotiriou, Phys. Lett. B 645, 389 (2007)
  * (39) G. J. Olmo and P. Singh, JCAP 0901, 030 (2009)
  * (40) C. Barragan, G. J. Olmo and H. Sanchis-Alepuz, Phys. Rev. D 80, 024016 (2009)
  * (41) T. S. Koivisto, Phys. Rev. D 82, 044022 (2010)
  * (42) S. D. Odintsov and V. K. Oikonomou, Phys. Rev. D 98, no. 2, 024013 (2018)
  * (43) A. V. Astashenok, S. Nojiri, S. D. Odintsov and A. V. Yurov, Phys. Lett. B 709, 396 (2012)
  * (44) G. Montani, Class. Quant. Grav. 18, 193 (2001)
  * (45) K. Dimopoulos, Phys. Lett. B 785, 132 (2018)
  * (46) L. H. Ford, Phys. Rev. D 35, 2955 (1987).
  * (47) M. M. Disconzi, T. W. Kephart and R. J. Scherrer, Phys. Rev. D 91, no. 4, 043532 (2015)
  * (48) G. Montani and M. Venanzi, Eur. Phys. J. C 77, no. 7, 486 (2017)
  * (49) W. Israel, Annals Phys. 100, 310 (1976).
  * (50) J. D. Barrow, Nucl. Phys. B 310, 743 (1988).
  * (51) V.A. Belinskii and I.M.‘Khalatnikov, Sov. J. Exp. Theor. Phys. 42, 205 (1975)
  * (52) V.A. Belinskii, E.S. Nikomarov and I.M. Khalatnikov, Sov. J. Exp. Theor. Phys. 50 (2), 21 (1979)
  * (53) F. Contreras, N. Cruz, E. Elizalde, E. González and S. Odintsov, arXiv:1808.06546 [gr-qc].
  * (54) N. Aghanim et al. [Planck Collaboration], arXiv:1807.06209 [astro-ph.CO].
  * (55) J. L. Lehners, Phys. Rept. 465, 223 (2008)


pagerefparensbackrefparens

# Latent Space Representations of Neural Algorithmic Reasoners

Vladimir V. Mirjanić
University of Cambridge
<EMAIL_ADDRESS>Pascanu
Google DeepMind
<EMAIL_ADDRESS>Veličković
Google DeepMind
<EMAIL_ADDRESS>
###### Abstract

Neural Algorithmic Reasoning (NAR) is a research area focused on designing
neural architectures that can reliably capture _classical computation_ ,
usually by learning to execute algorithms. A typical approach is to rely on
Graph Neural Network (GNN) architectures, which encode inputs in high-
dimensional latent spaces that are repeatedly transformed during the execution
of the algorithm. In this work we perform a detailed analysis of the structure
of the latent space induced by the GNN when executing algorithms. We identify
two possible failure modes: (i) loss of resolution, making it hard to
distinguish similar values; (ii) inability to deal with values outside the
range observed during training. We propose to solve the first issue by relying
on a _softmax_ aggregator, and propose to decay the latent space in order to
deal with out-of-range values. We show that these changes lead to improvements
on the majority of algorithms in the standard CLRS-30 benchmark when using the
state-of-the-art Triplet-GMPNN processor.

##  1 Introduction

Algorithms are one of the cornerstones of computer science. Recently, a large
body of work on Neural Algorithmic Reasoning showed how neural networks can be
taught to simulate classical algorithms [18, 14, 17]. This has emerged as a
highly popular framework, and its utility extends beyond replicating steps of
known algorithms. Neural networks trained on max-flow and min-cut have
improved brain vessel classification [13]. In reinforcement learning, similar
approaches are immune to many downsides of methods such as Value Iteration
Networks [5].

Neural Algorithmic Reasoners, often implemented as GNN architectures and
trained to simulate algorithms, use latent spaces to represent complex data,
with each of their message passing steps corresponding to a step in algorithm
execution. The evolution of the latent representation therefore corresponds to
the execution of an algorithm. In this work we exploit our understanding of
classical algorithms in order to analyse this space. Specifically, we use
dimensionality reduction techniques like PCA and explore perturbations of
algorithm inputs, to visualize and provide insights into the structures that
emerge in the representational space. We observe clusters under different
algorithmic symmetries, and we show that trajectories in the latent space tend
to converge to a single attractor, using dynamical systems.

Our analysis also reveals two weaknesses of the typical GNN architectures used
for algorithmic reasoning. The first one is that these systems struggle when
comparing similar values. We hypothesize that due to the _max_ operator used
for aggregation, the GNN tends to make a random choice between the values and
is able to propagate gradients only through that choice. Even if the choice
was suboptimal, it gets reinforced, and the learning process remains blind to
the fact that it could perform better. To address this we propose to use
_softmax_ as an aggregator, which in this scenario will allow gradient to
propagate through all pathways, and the learning process to discover the
optimal choice.

The second weakness is that the model tends to struggle when encountering out-
of-distribution values during algorithm execution. We propose that GNN should
decay magnitude of the representations at each step, allowing slightly out-of-
range values to become within range during the execution of the algorithm. We
show that these changes bring improvements to state-of-the-art models on the
majority of algorithms from the commonly used CLRS-30 benchmark [16, 8].

The success of our improvements indicates that a better understanding of
latent spaces is likely to be crucial to further improving GNN architectures.
In summary, our main contributions are

  * •

The first, comprehensive, study of the latent spaces emerging from training on
algorithmic tasks, which shows that

    * –

learned manifolds of graph embeddings are of much lower dimension compared to
the size of latent spaces

    * –

graphs with similar execution trajectories tend to cluster together

    * –

the embeddings converge to an attractor state corresponding to the end of
algorithm execution

  * •

Showing that GNNs fail to distinguish between branches of relatively similar
values, and when encountering OOD values during execution.

  * •

Softmax aggregation and processor decay as improvements, to ameliorate
discovered failure modes, with an evaluation on CLRS-30 showing improvement on
the majority of algorithms.

##  2 Background and Related Work

Consider the Bellman-Ford algorithm for finding shortest paths [1]. Its input
is a weighted graph with nodes $V$ and edges $E$, as well as a source node
$s$. It outputs $\pi:V\to V$, pointers to predecessor nodes along shortest
paths to source. Values $d:V\to\mathbb{R}$ denote the distances to source. The
main update rule (see Algorithm 1, lines 9-11) can be rewritten as in Eq. 1.

$d_{v}=\min\left(\hat{d}_{v},\
\min_{u\in\mathcal{N}_{v}}\left(\hat{d}_{u}+w(u,v)\right)\right)$ (1)

Here, $\hat{d}$ stands for the shortest distances from the previous step.
Compare it to a general message passing GNN in Eq. 2, with $\psi$ and $\phi$
as trainable MLPs [3].

$\mathbf{z}^{(i+1)}_{u}=\phi\left(\mathbf{z}^{(i)}_{u},\bigoplus_{v\in\mathcal{N}_{u}}\psi\left(\mathbf{z}_{u}^{(i)},\mathbf{z}^{(i)}_{v}\right)\right)$
(2)

This similarity is no coincidence; GNNs are closely aligned with dynamic
programming [17, 6]. Therefore, they are particularly suited for simulating
classical algorithms.

[14] uses this similarity to simulate algorithms by learning to replicate
individual steps of their execution. For example, in Bellman-Ford (Algorithm
1), one step corresponds to one iteration of the loop in lines 6-12, and the
GNN is trained to reproduce intermediate values of $\pi$ and $d$. These
intermediate values are called _hints_. They guide a GNN in learning the
algorithm, and they provide more immediate feedback than only using
algorithm’s outputs. In addition to hints, GNN is also given access to latent
vectors in a high-dimensional space, to simulate working memory.

Algorithm 1 Bellman-Ford

1:

2:function Bellman-Ford($\mathcal{G}(V,E)$, $s$)

3: for $v\in V$ do

4: $\pi_{v}$ $\leftarrow$ None

5: $d_{v}$ $\leftarrow$ $\infty$

6: $d_{s}$ $\leftarrow$ 0

7: repeat

8: $\hat{d}$ $\leftarrow$ $d$

9: for $e(u,v)\in E$ do

10: if $d_{v}>d_{u}+w(u,v)$ then

11: $d_{v}$ $\leftarrow$ $d_{u}+w(u,v)$

12: $\pi_{v}$ $\leftarrow$ $u$

13: until $\hat{d}=d$

14:return $\pi$ hg Step of algorithm execution

Several different approaches to aligning neural networks with algorithms
exist. On a microscopic scale, there are Neural Execution Engines [18], which
learn basic operations, such as summation, $\min$, or $\max$, and combine them
into larger compositional blocks. Macroscopic approaches also exist [17], and
they learn algorithms end-to-end. The model we will focus on in this work,
proposed by [14], that we will refer to as Neural Algorithmic Reasoner (NAR),
operates on a mesoscopic scale, and learns key steps in algorithm execution.

These approaches target different levels of abstraction, and therefore their
latent spaces learn to represent different types of data. When latent spaces
of Neural Execution Engines are visualised with PCA, sequential natural
numbers are found organised in chains in sorting tasks, and “human-
interpretable patterns” are found for other tasks as well [18]. These patterns
extend to values not seen during training, even with large holdout. Thus,
these networks can structure their latent spaces even with low amounts of
training data.

The reason we focus on the NAR is that it allows multi-task learning, while
Neural Execution Engines are hard-wired for a specific algorithm. In fact, NAR
architecture can learn all 30 algorithms from the CLRS-30 benchmark [16] at
once, and this even leads to improvement on some tasks, compared to learning
them separately [8]. This suggests a much richer latent structure, which can
provide insights that hold across multiple algorithms.

##  3 Latent Representations

We focus on a single architecture and a single algorithm. For the algorithm,
we pick Bellman-Ford (as implemented in Algorithm 1). Bellman-Ford trains
quickly, since the number of algorithm steps is bounded by the number of
nodes, and it uses little GPU memory while training. These properties make it
a good choice for rapid prototyping of many different ideas.

For the network, we create a variation of PGN [15]. We simplify the network by
removing all ReLU non-linearities, so that the only nonlinearity remaining is
the $\max$ aggregation inherent in the layer. We name the result LinearPGN.
This new processor has simpler representations than original PGN, due to lack
of non-linearities.

Despite being structurally simple, LinearPGN is expressive enough to simulate
Bellman-Ford. Min aggregation over the edges in Bellman-Ford can be converted
into max aggregation used in message passing by negating the messages both
before and after it. The messages themselves are simply the distance updates
$\hat{d}_{u}+w(u,v)$. This can be seen in Eq. 3, which at the same time
corresponds both to the Bellman-Ford update rule, and to fully linear message
passing.

$\mathbf{z}_{v}^{(i+1)}=-\max_{u}\left(-\mathbf{z}_{u}^{(i)}-w(u,v)\right)$
(3)

###  3.1 Trajectories of Embeddings

To visualise the latent spaces, we record the complete _trajectories_ of input
graphs’ latent embeddings, i.e., we save the hidden vectors
$\mathbf{z}_{v}^{(i)}$ for each node $v$ in the graph, for each time step $i$.
Some graphs may halt their execution earlier than others, and this skews the
distribution of hidden vectors. To avoid this, we filter out a subset of
graphs that all terminate at the same time $T$. Thus, trajectory data has
shape $N\times|V|\times D\times T$, where $N$ is the number of sampled graphs,
$|V|$ is the number of nodes in each graph, $D$ is the dimensionality of the
latent space, and $T$ is the number of steps that the processor executed
before terminating. We first reduce the number of dimensions by eliminating
the (second) $|V|$ axis. We opt to use $\max$-aggregation, so that for each
step $i$ we extract $\max_{v}\mathbf{z}_{v}^{(i)}$ out of the individual node
embeddings. We provide an ablation to validate our choice in Appendix C.2.

After this reduction, the trajectory tensor has shape $N\times D\times T$. We
consider two different types of plot by

  * •

combining the dimension and time step axes; this leads to a matrix of size
$N\times(D\cdot T)$ over which we can apply PCA. We call this the trajectory-
wise plot because entire trajectories for each of $N$ input graphs are
summarised into a single point. We reduce the number of $DT$ down to $d=2$ or
$d=3$.

  * •

grouping time step and sample axes against the latent space, in order to
obtain a matrix of size $(N\cdot T)\times D$. We call this the step-wise plot
because each graph is represented with $T$ points – one for each execution
step. We typically trace red lines through the $T$ points of each graph, in
order to represent trajectories in the latent space, and color the points
based on the time step.

On all PCA visualisations, X axis is the axis of the highest variance, and Y
axis has the second highest variance.

Figure 1: Left: Step-wise PCA visualisation of latent spaces for random
graphs. Middle: Step-wise PCA for a set of graphs derived from scaling
symmetry. Scaling has one degree of freedom. This is captured in the latent
space, as all step embeddings are one-dimensional. Right: Step-wise PCA of 8
classes of random graphs, each containing graphs equivalent under reweighting
symmetry. Clusters corresponding to each class of graphs are clearly visible,
showing that the model learns to group graphs with same execution close
together. Red lines trace graph trajectories. Different colors correspond to
different steps in algorithm execution. Sequential view available in Appendix
F.

###  3.2 Dimensionality

We first study the robustness of latent representations. The latent space in
NAR is $128$-dimensional ($D=128$), but valid trajectories might live in a
considerably lower-dimensional space, where the redundancy might help
robustness to noise.

We measure the amount of variance captured by three most dominant dimensions
for trajectory-wise and step-wise PCA. We sample random Erdős-Rényi graphs,
and measure that in trajectory-wise PCA the three dimensions capture $63.5\%$
of variance. For step-wise PCA (Fig. 1 left), this number is $96.4\%$. If the
embeddings were uniformly spread, then each PCA dimension would account for
only $\frac{1}{128}$ of the total variance, and therefore the first three
dimensions would contribute with only $2.3\%$. The $63.5\%$ that they
contribute with is therefore quite high, and indicates that the
representations are low-dimensional. The even higher value of over $96.4\%$
for step-wise PCA shows that there is a large amount of structure in the step-
wise view of the latent space, and that the movement between steps is
significantly stronger than the variance in graph trajectories.

###  3.3 Representations Under Symmetries

GNNs are designed to be invariant or equivariant to graph permutations.
Therefore, they can operate on all types of graphs, including algorithm
inputs. However, algorithms and their data often have more structure than mere
permutation invariance. After showing the structure of latent spaces for
random graphs (Fig. 1 left), we explore how the model deals with additional
symmetries of the learned algorithm. Recall that inputs to Bellman-Ford are
weighted graphs, with graph weights in the interval $(0,1)$. There are many
non-trivial transforms that we can perform on the inputs that keep algorithm
execution constant, but change the embeddings in the latent space.

We start by formalizing the concept of keeping algorithm execution constant.
We take it to mean that if we apply a transformation to the inputs, then hints
at each time step should change accordingly. Bellman-Ford has two kinds of
hints, $\pi:V\to V$, pointers to other nodes, and $d:V\to\mathbb{R}$,
distances to source. The difference between them impacts how they change with
changes in inputs. By analogy with the permutation equivariance of GNNs, we
require that both $\pi$ and $d$ be equivariant to permutations of input nodes.
Meanwhile, $\pi$ should be invariant and $d$ equivariant to transformations
that change weights only. To understand how symmetries are encoded in the
latent space, we study collections of inputs that have equivalent algorithm
executions. We skip over permutation symmetries because they are general to
all GNNs, and look at those symmetries specific to Bellman-Ford.

####  3.3.1 Scaling Symmetry

First we consider the scaling symmetry of graph weights $w_{uv}\mapsto\lambda
w_{uv}$, where $\lambda>0$, which trivially preserves algorithm execution. As
discussed above, this symmetry preserves pointers $\pi$, while transforming
distances $d$ as $d_{u}\mapsto\lambda d_{u}$. To visualise the symmetry (Fig.
1 middle), we generate one random graph, and then scale it with
$\lambda\in(\frac{1}{2},1)$. We avoid $\lambda>1$ as that would create edge
weights larger than those seen in training. We also avoid small $\lambda$
because $w=0$ is used to encode lack of edges, and not edges of weight $0$.
This means that representations around $w=0$ can behave unpredictably.

Note that scaling by positive $\lambda$ is necessary because introducing
negative values creates negative cycles that dramatically change algorithm
execution. Also, observe that adding a constant $c$ to each edge is not an
algorithm symmetry. A single degree of freedom of this symmetry, $\lambda$,
gives rise to one-dimensional embeddings (Fig. 1 middle). The embeddings are
not in a straight line, but rather on a non-trivial curve. This is not as
visible in the step-wise plot, as movement between steps overshadows variance
between graphs.

####  3.3.2 Reweighting Symmetry

Next we turn to a more complex symmetry of Bellman-Ford. We name it
reweighting symmetry after the “reweighting” trick from the Johnson’s
algorithm [9]. We outline the trick in Lemma 3.1.

######  Lemma 3.1.

(Reweighting trick) Let $h:V\to\mathbb{R}$ assign arbitrary values to nodes.
Suppose the input graph $\mathcal{G}$ with source node $s$ is reweighted as
follows

$\hat{w}(u,v)\triangleq w(u,v)+h(u)-h(v)$

Then the new re-weighted graph $\mathcal{\hat{G}}$ has identical execution
trace, meaning that at each execution step $\hat{\pi}=\pi$ and
$\hat{d}(u)=d(u)+h(s)-h(u)$ hold, i.e. while values change, branches taken by
the algorithm do not.

Figure 2: Left: Trajectory-wise PCA of eight clusters of reweighted graphs
showing that they all contain a single dominant direction. Different clusters
have different colors. Middle: Many embedding clusters with dominant
directions overlaid in red. Right: Step-wise PCA of random graphs with the
dominant cluster directions overlaid in red.

In Johnson’s algorithm, this trick removes negative edges from the graph, so
that Dijkstra’s algorithm can be used to find shortest paths. We use the trick
to generate many graphs with the same execution.

The reweighting trick operates on graphs with real-valued weights, but our
Bellman-Ford NAR is trained on graphs with weights in $(0,1)$. Therefore, we
must ensure that weights stay in this distribution when applying Lemma 3.1 to
generate perturbations. That is why we sample weights in the reduced range
$(c,1-c)$, and then sample $h$ randomly from $(0,c)$. Thus, we preserve high
variance between weights from the original graph and can safely apply
reweighting, as we are guaranteed to stay in the $(0,1)$ range afterwards. We
experiment with $c=\frac{1}{2}$ and $c=\frac{1}{4}$, and observe no
difference.

We sample several such random graphs and then generate multiple variations
through reweighting (Fig. 1 right). These reweighting symmetries have $|V|$
degrees of freedom because each node $v$ gets assigned an independent value
$h(v)$. In line with that, variance of $94\%$ for the three most dominant
dimensions of reweighted graphs is similar to that of random graphs. However,
it is highly significant that these clusters are clearly separated.

###  3.4 Cluster Directions in Latent Spaces

Another interesting observation to be made regarding reweighted graphs is that
their embeddings seem to have one main direction of variance. This is more
easily observed in the trajectory-wise plot (Fig. 2 left) than in the step-
wise plot. Note that the variance along these directions is larger than the
difference between cluster embeddings. We are interested in learning whether
these directions are parallel or perpendicular to the manifold describing the
execution of the algorithm. If they are perpendicular, moving along them
should not affect the execution of the algorithm and hence the behavior of the
model.

To proceed, we build a database of hundreds of reweighting clusters in order
to cover as much of the embeddings manifold as possible. Then, we apply PCA to
each cluster at each time step, and extract the first several principal
components. We use this fine-grained approach because the manifold itself is
curved, and therefore these directions vary along it (Fig. 2 right).
Qualitatively, these principal directions seem to be parallel to the manifold,
which would mean that moving along them changes the meaning of the embeddings.
Nevertheless, we perform a quantitative analysis using our database of
directions.

We modify the embeddings by adding noise at test-time along these principal
directions versus along random directions. We also consider removing the
variance along the directions of interest, and also preserving only the
variance along these directions. To pick the correct direction representing
the symmetry, we pick the principal direction of the closest cluster from the
database. In a control experiment, we pick a mean of all directions instead.

Table 1: Accuracy with perturbations along principal directions of variance.
Perturbations include: the noise-free baseline (Noise-free), Gaussian noise
along the direction (Directional), Gaussian noise along random direction
(Random), projecting out the direction (Out), and projecting the embedding
onto the direction (Onto).

| Perturbation Type
---|---
Direction | Noise-free | Directional | Random | Out | Onto
L2 Closest |  |  | Mean |  |

Adding noise along these directions is more harmful than adding noise along a
random direction. Also, projecting embeddings onto the directions is better
than projecting the directions out of them (Tab. 3.4). Both of these results
contradict the hypothesis that these directions might be perpendicular to the
manifold. One interpretation could be that the manifold is highly curved, or
that we can view it as divided into small patches. Any movement within the
path preserves execution, while larger movements reach different patches with
different algorithm trajectories.

##  4 Mispredict Analysis

We now investigate how close the model learns to be aligned with the ground-
truth Bellman-Ford algorithm. To do so, we consider three faulty
implementations of the Bellman-Ford algorithm, all of which represent
simplifications that the model might have learned, and mostly change how
distances are being compared (Alg. 1 line 9). Firstly, we discuss the Greedy-
Ford algorithm, which compares only weights and not full distances, and is
used to see whether the model learns a simpler algorithm than what is
required. Secondly, the Decay-Ford algorithm scales distances by a constant,
and is used to check whether the model learns the correct algorithm while
falsely learning a few weights. Finally, the Noisy-Ford algorithm, which adds
small amount of noise to distances, is there to check whether errors of the
model result from noise in the embeddings.

Figure 3: Left: Similarity between final predictions for our model, true
Bellman-Ford, and multiple faulty implementations. Our model is closer to
ground truth than to any others. Right: Predictions of the model split on the
basis of whether they match ground truth at the end (y-axis), and during
execution (x-axis). Almost all mispredicts were correctly predicted at an
earlier point.

Our model is close in its predictions to all variants except Greedy-Ford,
which means that it is learning something quite similar to the true algorithm
(Fig. 3 left). In fact, out of all variants we have considered, our model is
the closest to ground truth. After this, we separate predictions based on
whether they are correct at the end of algorithm execution, and whether they
are correct at any intermediate step (Fig. 3 right). Intermediate predictions
also include final predictions, which means that the upper-right field is
exactly $0\%$ (if our model’s prediction for a node is never correct, then it
is also not correct in the end). However, the bottom-left field is
surprisingly high.

This shows that nearly all mispredicts were predicted correctly at some point
during the simulation, and were then changed to a wrong value afterwards. Only
one percent of nodes were never predicted correctly. The ground-truth Bellman-
Ford algorithm, meanwhile, never changes the correct output.

We observe that these wrong predictions typically occur when the network has
to decide between two different paths with similar distances, and give an
example in Appendix E. To address this problem, we propose as a solution to
use softmax aggregation instead of hard $\max$. With $\max$ aggregation,
derivatives are backpropagated only through the largest value, even if others
are very close to it, and information can be lost. With softmax, all values
impact the output and all values are backpropagated through. The effect would
be particularly strong in cases where there are two large values and many
small ones. We regard this as one main outcome of our analysis and validate it
at scale in Section 5.

###  4.1 Mispredicts and Value Generalisation

Consider random weighted graphs where each pair of nodes is connected with
probability $p$, such as those used in our evaluations. Let weights be sampled
uniformly from $(0,1)$, and $p=0.5$. We measure the average distance between
nodes to be $0.15$. However, if $p=0.25$, the average distance rises to $1.1$.

This shift in graph statistics means that a model trained on random graphs
where $p_{avg}=0.5$ does not see long distances as often as short ones and
struggles when handling them. Indeed, as we change connectivity of graphs from
$p=0.5$ to $p=0.25$ during testing, we notice a drop in accuracy from $93\%$
to $87\%$. On the other hand, if we also scale weights down before evaluating,
to mitigate the distance increase, the accuracy is restored to previous
levels. This shows that the drop in performance is due to the model failing to
handle paths longer than those it has seen during training.

Figure 4: Left: Misprediction rate as a function of the shortest path
distances. Its theoretical range is from $0$ to $1$, and lower is better.
Right: Misprediction rate as a function of the _predicted_ distances.

In fact, wrong predictions of the model on Bellman-Ford algorithm with
LinearPGN as the processor are directly correlated with the nodes’ distances
to source (Fig. 4). We observe that the model is failing to correctly
reproduce distances larger than $1$. Furthermore, we show that misprediction
rate steadily increases for predicted distances between $0.75$ and $1.5$, and
for paths with predicted distance above $1.5$ the model wrongly predicts every
single one (Fig. 4 right).

As the model is struggling with large out-of-distribution values, we propose
using a decay-like regularisation where, at every message passing step, we
scale the embeddings by a constant $c<1$. We show in the following section
that this provides improvements not only on the Bellman-Ford algorithm, but on
other algorithms in the CLRS-30 benchmark [16] as well.

##  5 Evaluation on CLRS-30

Our investigation produced two proposed changes that have the potential to
improve the ability of GNN to model algorithms.

The first is softmax aggregation. In GNNs, once all messages to node $v_{i}$
are computed, they are aggregated in order to produce the new latent
embeddings. Common choices for aggregation include summing over the messages,
or finding the mean, but $\max$ aggregation is most commonly used in
algorithmic reasoning tasks. That is, latents are calculated via
$\mathbf{z}_{i}=\max_{j}\mathbf{m}_{ji}$. The way messages $\mathbf{m}$ are
computed depends on the layer (PGN [15], Triplet-GMPNN [8], etc.). We argue
that $\max$ aggregation cuts off messages that are close but not equal to the
maximum value, and that this explains why the model struggles when comparing
paths of similar lengths in the Bellman-Ford algorithm.

Softmax aggregation [10] weighs each message based on its softmax coefficient
and then sums them up:

$z_{i}=\sum_{j}\sigma(\frac{\mathbf{m}_{-i}}{T},j)\cdot\mathbf{m}_{ji}$ (4)

We write softmax of $x_{1}$, $x_{2}$, …, as
$\sigma(x_{-},i)=\frac{\exp{x_{i}}}{\sum_{j}\exp{x_{j}}}$. This is further
parameterised with temperature $T$. As $T\to 0$, softmax approaches the values
of (hard) $\max$, while for $T=1$ it behaves as standard unparametrised
softmax, and as $T\to\infty$ it behaves as a mean aggregator.

The second improvement we proposed is the use of processor decay. This means
that after each reduction, we scale the values down with a constant factor
$c$.

These modifications can be applied individually, or together. For evaluating
them, we use the CLRS-30 benchmark [16], which contains a curated list of 30
classical algorithms [4]. These algorithms cover a wide range of concepts –
searching, sorting, dynamical programming, divide and conquer, greedy, as well
as many graph-specific topics, such as shortest paths, spanning trees, and
more. All algorithms are tested on graphs that are four times larger than
those seen during training. We estimate standard deviation with five runs. We
use Triplet-GMPNN [8] as the NAR processor for evaluation, as it is currently
state-of-the-art on CLRS-30.

Table 2: Performance of Triplet-GMPNN with pipeline modifications on CLRS-30 algorithms. Pipeline modifications are softmax aggregation (sft) and processor decay (dec). Summary | Baseline | +sft | +dec | +sft+dec
---|---|---|---|---
Best on #/30 | 11 | 8 | 2 | 9
Above baseline | N/A | 15 | 10 | 11

We observe that our changes achieve top accuracy on the majority of tasks
(Tab. 2). Furthermore, softmax is most advantageous when used alone. Decay is
more situational, and more work is needed in order to understand how to
separate its benefits from the downsides. We use no hyperparameter tuning for
training. Softmax temperature is set to $0.01$, and processor is decayed with
a factor of $0.9$.

Finally, we perform multiple ablations. We show that small decay and small
temperature perform better than large ones. We also show that our analysis of
latent spaces extends to Triplet-GMPNN. Full results can be found in the
Appendix.

##  6 Discussion

Conclusions We provide an extensive analysis of the latent representations
learned by GNN trained on algorithmic reasoning tasks. We show that these
latent spaces have rich structure, in particular in how they reflect different
symmetries in executing algorithms. We show that the execution trajectories in
the latent space tend to converge to an attractor-like state, and graphs of
similar executions tend to cluster together. Furthermore, from our analysis we
identify two issues with current architectures: (i) dealing with values of
similar magnitude, and (ii) dealing with values that are out of range compared
to those seen during training.

We hypothesise that the first issue is due to the use of the max aggregator
function, which backpropagates gradients only along the largest of the similar
values, making it harder for the learning process to identify whether it made
a suboptimal choice. We propose to use _softmax_ instead of _max_ as
aggregator, allowing gradients to flow through all paths proportional to their
magnitude. We empirically validate this choice at scale on the CLRS-30
benchmark.

The second issue for the Bellman-Ford algorithm happens when accumulating
distances between nodes. The issue is that depending on the graph connectivity
the distribution of distances and the embeddings in latent space can change
drastically. We propose a simple fix – decaying the magnitude of the embedding
by a fixed rate at every execution step. This allows the embeddings to
consistently stay in a similar range. Decay, coupled with the choice of
softmax aggregator, provides improvement across many algorithms in CLRS-30
benchmark.

We regard our analysis as the first step in trying to understand how neural
architectures model the algorithms they are trained to mimic. Our work shows
that such analysis can reveal weaknesses of current parametrisations, and it
points towards potential ways of addressing them.

Future Work We demonstrated that the latent space contains a single attractor
to which the trajectories converge. This convergence could be safely exploited
as a termination criterion. Currently, all experiments using CLRS-30 rely on
knowing the exact number of steps needed for each input, even at test time.
This is unrealistic, as such information is not available in the real world.
An alternative, explored in [14], is to use a dedicated neural network to
learn when to terminate. The authors proposed using embeddings after each step
as an input to the network. Our research, meanwhile, motivates the use of
differences between successive encodings. Another direction for future work is
investigating the impact of encoder and decoder on our analysis. In this
project, we mainly focused on the processor, assuming that it learns the bulk
of the algorithm. However, the importance of encoders and decoders could also
be investigated (e.g. by freezing them during training). Finally, this study
can be extended to multi-task learning.

## Acknowledgements

We would like to thank Prof. Pietro Liò (University of Cambridge) for
supervising the project while it was being done in partial fulfillment of the
requirements for the Computer Science Tripos, Part III (MEng). We would like
to thank Dr. Jovana Mitrović and Prof. Karl Tuyls (Google DeepMind) for
conducting the tech and strategy reviews, and providing useful comments.

## References

  * [1] Richard Bellman  “On a Routing Problem” Publisher: Brown University  In _Quarterly of Applied Mathematics_ 16.1, 1958, pp. 87–90  URL: https://www.jstor.org/stable/43634538
  * [2] Christopher M. Bishop  “Pattern Recognition and Machine Learning” Google-Books-ID: kTNoQgAACAAJ  Springer, 2006
  * [3] Michael M. Bronstein, Joan Bruna, Taco Cohen and Petar Veličković  “Geometric Deep Learning: Grids, Groups, Graphs, Geodesics, and Gauges” arXiv:2104.13478 [cs, stat]  arXiv, 2021  URL: http://arxiv.org/abs/2104.13478
  * [4] Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest and Clifford Stein  “Introduction to Algorithms, Third Edition”  The MIT Press, 2009
  * [5] Andreea Deac et al.  “XLVIN: eXecuted Latent Value Iteration Nets” arXiv:2010.13146 [cs, stat]  arXiv, 2020  DOI: 10.48550/arXiv.2010.13146
  * [6] Andrew Dudzik and Petar Veličković  “Graph Neural Networks are Dynamic Programmers” arXiv:2203.15544 [cs, math, stat]  arXiv, 2022  URL: http://arxiv.org/abs/2203.15544
  * [7] Ian Goodfellow, Yoshua Bengio and Aaron Courville  “Deep Learning” Google-Books-ID: omivDQAAQBAJ  MIT Press, 2016
  * [8] Borja Ibarz et al.  “A Generalist Neural Algorithmic Learner” arXiv:2209.11142 [cs, stat]  arXiv, 2022  URL: http://arxiv.org/abs/2209.11142
  * [9] Donald B. Johnson  “Efficient Algorithms for Shortest Paths in Sparse Networks”  In _Journal of the ACM_ 24.1, 1977, pp. 1–13  DOI: 10.1145/321992.321993
  * [10] Guohao Li, Chenxin Xiong, Ali Thabet and Bernard Ghanem  “DeeperGCN: All You Need to Train Deeper GCNs” arXiv:2006.07739 [cs, stat]  arXiv, 2020  DOI: 10.48550/arXiv.2006.07739
  * [11] Wolfgang Maass, Prashant Joshi and Eduardo Sontag  “Computational aspects of feedback in neural circuits”  In _PLoS Computational Biology_ 3.1 Public Library of Science, 2007, pp. 1–20
  * [12] Laurens van der Maaten and Geoffrey Hinton  “Visualizing Data using t-SNE”  In _Journal of Machine Learning Research_ 9.86, 2008, pp. 2579–2605  URL: http://jmlr.org/papers/v9/vandermaaten08a.html
  * [13] Danilo Numeroso, Davide Bacciu and Petar Veličković  “Dual Algorithmic Reasoning” arXiv:2302.04496 [cs]  arXiv, 2023  URL: http://arxiv.org/abs/2302.04496
  * [14] Petar Veličković et al.  “Neural Execution of Graph Algorithms”, 2020  URL: https://iclr.cc/virtual_2020/poster_SkgKO0EtvS.html
  * [15] Petar Veličković et al.  “Pointer Graph Networks”  In _Advances in Neural Information Processing Systems_ 33 Curran Associates, Inc., 2020, pp. 2232–2244
  * [16] Petar Veličković et al.  “The CLRS Algorithmic Reasoning Benchmark” arXiv:2205.15659 [cs, stat]  arXiv, 2022  DOI: 10.48550/arXiv.2205.15659
  * [17] Keyulu Xu et al.  “What Can Neural Networks Reason About?”, 2020  URL: https://iclr.cc/virtual_2020/poster_rJxbJeHFPS.html
  * [18] Yujun Yan et al.  “Neural Execution Engines: Learning to Execute Subroutines”  In _Advances in Neural Information Processing Systems_ 33 Curran Associates, Inc., 2020, pp. 17298–17308

##  Appendix A Detailed Results on CLRS-30

Here we evaluate softmax aggregation and processor decay on Triplet-GMPNN and
MPNN architectures. Best results are marked in bold, while (non-best) results
above baseline are blue.

###  A.1 Results on Triplet-GMPNN

Table 3: Performance of Triplet-GMPNN with pipeline modifications on CLRS-30 algorithms. Pipeline modifications are softmax aggregation (sft) and processor decay (dec). Standard deviation estimated with five runs. We observe minor differences between the baseline accuracy as measured by us and that of [8]. We expect that the cause lies in [8] performing 10 runs as opposed to our five runs. | Triplet-GMPNN
---|---
Algorithm | baseline | +sft | +dec | +sft+dec
Activity Selector |  ${\mathbf{96.46\%}}$$\pm 1.64$ | Articulation Points |  ${\mathbf{73.10\%}}$$\pm 21.77$ |  $\mathcolor{blue}{68.28\%}$$\pm 7.89$ | Bellman-Ford |  $\mathcolor{blue}{98.80\%}$$\pm 0.14$ |  ${\mathbf{98.96\%}}$$\pm 0.35$
BFS |  ${\mathbf{99.63\%}}$$\pm 0.10$ |  $\mathcolor{blue}{99.54\%}$$\pm 0.16$ |  $\mathcolor{blue}{99.50\%}$$\pm 0.09$ |
Binary Search |  $\mathcolor{blue}{64.64\%}$$\pm 8.75$ |  ${\mathbf{66.49\%}}$$\pm 6.29$ |  |
Bridges |  $\mathcolor{blue}{77.83\%}$$\pm 14.64$ |  $\mathcolor{blue}{74.50\%}$$\pm 30.95$ |  ${\mathbf{84.10\%}}$$\pm 18.15$ |
Bubble Sort |  ${\mathbf{66.34\%}}$$\pm 6.29$ | DAG Shortest Paths |  ${\mathbf{89.17\%}}$$\pm 6.20$ | DFS |  ${\mathbf{19.29\%}}$$\pm 7.94$ | Dijkstra |  ${\mathbf{97.99\%}}$$\pm 0.24$ | Find Maximum Subarray |  ${\mathbf{62.84\%}}$$\pm 6.45$ | Floyd-Warshall |  $\mathcolor{blue}{25.77\%}$$\pm 16.82$ |  $\mathcolor{blue}{18.90\%}$$\pm 9.45$ |  ${\mathbf{32.52\%}}$$\pm 13.77$
Graham Scan |  ${\mathbf{93.39\%}}$$\pm 3.04$ | Heapsort |  ${\mathbf{42.80\%}}$$\pm 18.85$ | Insertion Sort |  ${\mathbf{69.42\%}}$$\pm 7.87$ | Jarvis March |  ${\mathbf{90.86\%}}$$\pm 1.73$ | KMP Matcher |  $\mathcolor{blue}{11.20\%}$$\pm 9.77$ |  ${\mathbf{20.59\%}}$$\pm 11.94$
LCS Length |  $\mathcolor{blue}{87.41\%}$$\pm 0.92$ |  ${\mathbf{87.75\%}}$$\pm 0.97$ |  |
Matrix Chain Order |  ${\mathbf{92.90\%}}$$\pm 1.49$ |  $\mathcolor{blue}{90.45\%}$$\pm 2.77$ |  |
Minimum |  $\mathcolor{blue}{95.77\%}$$\pm 1.06$ |  $\mathcolor{blue}{95.46\%}$$\pm 5.78$ |  ${\mathbf{96.08\%}}$$\pm 2.07$ |
MST Kruskal |  ${\mathbf{88.83\%}}$$\pm 2.29$ | MST Prim |  ${\mathbf{87.42\%}}$$\pm 2.02$ | Naïve String Matcher |  ${\mathbf{17.50\%}}$$\pm 12.11$ | Optimal BST |  ${\mathbf{82.93\%}}$$\pm 1.84$ | Quickselect |  $\mathcolor{blue}{1.37\%}$$\pm 2.14$ |  ${\mathbf{1.62\%}}$$\pm 1.21$
Quicksort |  ${\mathbf{49.22\%}}$$\pm 13.53$ |  $\mathcolor{blue}{46.46\%}$$\pm 24.85$ |  |
Segments Intersect |  $\mathcolor{blue}{98.74\%}$$\pm 0.22$ |  $\mathcolor{blue}{98.68\%}$$\pm 0.29$ |  ${\mathbf{98.81\%}}$$\pm 0.15$ |
Strongly Connected |  ${\mathbf{43.63\%}}$$\pm 7.94$ | Task Scheduling |  ${\mathbf{87.58\%}}$$\pm 0.66$ | Topological Sort |  ${\mathbf{71.29\%}}$$\pm 3.25$

###  A.2 Results on MPNN

On Table 4 we display our results for MPNN instead of Triplet-GMPNN. Generally
speaking, improvements are more significant compared to Triplet-GMPNN.

Table 4: Performance of MPNN with pipeline modifications on CLRS-30 algorithms. Pipeline modifications are softmax aggregation (sft) and processor decay (dec). Standard deviation estimated with five runs. | MPNN
---|---
Algorithm | baseline | +sft | +dec | +sft+dec
Activity Selector |  $\mathcolor{blue}{90.17\%}$$\pm 4.49$ |  ${\mathbf{92.10\%}}$$\pm 1.59$ |  $\mathcolor{blue}{91.56\%}$$\pm 2.04$ |
Articulation Points |  $\mathcolor{blue}{46.94\%}$$\pm 7.27$ |  ${\mathbf{47.37\%}}$$\pm 13.24$ |  |
Bellman-Ford |  ${\mathbf{93.37\%}}$$\pm 1.33$ |  $\mathcolor{blue}{92.76\%}$$\pm 0.72$ |  $\mathcolor{blue}{92.91\%}$$\pm 1.04$ |
BFS |  $\mathcolor{blue}{99.71\%}$$\pm 0.10$ |  $\mathcolor{blue}{99.69\%}$$\pm 0.16$ |  ${\mathbf{99.74\%}}$$\pm 0.14$ |
Binary Search |  $\mathcolor{blue}{48.85\%}$$\pm 7.60$ |  $\mathcolor{blue}{54.30\%}$$\pm 8.30$ |  ${\mathbf{56.89\%}}$$\pm 9.37$ |
Bridges |  $\mathcolor{blue}{23.54\%}$$\pm 13.19$ |  ${\mathbf{24.02\%}}$$\pm 13.43$ |  $\mathcolor{blue}{23.83\%}$$\pm 13.34$ |
Bubble Sort |  $\mathcolor{blue}{65.55\%}$$\pm 12.71$ |  ${\mathbf{66.82\%}}$$\pm 6.76$ | DAG Shortest Paths |  ${\mathbf{92.01\%}}$$\pm 0.95$ |  $\mathcolor{blue}{85.42\%}$$\pm 7.51$ |  $\mathcolor{blue}{84.05\%}$$\pm 3.39$
DFS |  $\mathcolor{blue}{14.66\%}$$\pm 7.05$ |  $\mathcolor{blue}{16.05\%}$$\pm 2.53$ |  ${\mathbf{16.64\%}}$$\pm 6.04$ |
Dijkstra |  ${\mathbf{93.93\%}}$$\pm 1.41$ | Find Maximum Subarray |  ${\mathbf{51.66\%}}$$\pm 7.19$ | Floyd-Warshall |  $\mathcolor{blue}{16.74\%}$$\pm 11.31$ |  $\mathcolor{blue}{18.74\%}$$\pm 7.18$ |  ${\mathbf{19.37\%}}$$\pm 4.00$
Graham Scan |  ${\mathbf{93.87\%}}$$\pm 1.22$ |  $\mathcolor{blue}{93.13\%}$$\pm 2.28$ |  |
Heapsort |  $\mathcolor{blue}{66.51\%}$$\pm 4.48$ |  ${\mathbf{68.87\%}}$$\pm 6.86$ | Insertion Sort |  $\mathcolor{blue}{60.50\%}$$\pm 7.06$ |  ${\mathbf{66.64\%}}$$\pm 5.33$
Jarvis March |  ${\mathbf{95.29\%}}$$\pm 1.37$ | KMP Matcher |  ${\mathbf{8.83\%}}$$\pm 3.11$ | LCS Length |  ${\mathbf{84.76\%}}$$\pm 2.45$ | Matrix Chain Order |  ${\mathbf{84.80\%}}$$\pm 1.37$ |  $\mathcolor{blue}{84.48\%}$$\pm 0.95$
Minimum |  ${\mathbf{95.62\%}}$$\pm 0.15$ | MST Kruskal |  $\mathcolor{blue}{66.65\%}$$\pm 11.25$ |  ${\mathbf{68.41\%}}$$\pm 8.95$ | MST Prim |  ${\mathbf{72.39\%}}$$\pm 11.88$ | Naïve String Matcher |  $\mathcolor{blue}{6.01\%}$$\pm 5.92$ |  ${\mathbf{7.03\%}}$$\pm 6.74$
Optimal BST |  ${\mathbf{74.02\%}}$$\pm 2.80$ | Quickselect |  ${\mathbf{4.07\%}}$$\pm 5.34$ | Quicksort |  ${\mathbf{61.45\%}}$$\pm 17.25$ | Segments Intersect |  ${\mathbf{94.09\%}}$$\pm 0.16$ | Strongly Connected |  ${\mathbf{32.01\%}}$$\pm 8.96$ |  $\mathcolor{blue}{30.99\%}$$\pm 7.06$ |  $\mathcolor{blue}{30.39\%}$$\pm 6.12$
Task Scheduling |  ${\mathbf{81.71\%}}$$\pm 1.54$ | Topological Sort |  $\mathcolor{blue}{63.11\%}$$\pm 10.34$ |  $\mathcolor{blue}{67.81\%}$$\pm 12.78$ |  ${\mathbf{74.81\%}}$$\pm 2.15$

##  Appendix B Attractors in Latent Spaces

We now view latent representations through the lens of dynamical systems. In
dynamical systems, a function governs how points move in ambient space over
time. The connection with NARs is clear – points are embeddings, the ambient
space is the latent space, and steps of NAR are a discretisation of the
function. We have already relied on this analogy implicitly, and have even
been borrowing terminology such as _trajectories_ of points/embeddings.
However, it is useful to make the connection explicit as this can reveal
additional information about the latent spaces.

One important concept in the area of dynamical systems is that of attractors.
Attractors are subsets of the latent space that are invariant under the update
function and to which nearby embeddings are drawn to. Attractors can be
points, cycles, or, as in the famous example of the Lorentz attractor,
neither. Empirically (Fig. 1), we notice that the first few steps of execution
create large movements in the latent space, and all steps from fifth onward
seem to blend into each other. This suggests some form of convergence towards
a single attractor. This can also be seen in terms of distances travelled at
each time step, which quickly converge to zero. Note that different equivalent
instances of the problem do not converge to the same point; we can view this
as some form of partial (or multidimensional) attractor [11].

Figure 5: Left: Model accuracy versus the step count of Bellman-Ford
execution. Graphs with lower step count perform better. Right: Model accuracy
as their NAR execution is forcefully increased or decreased. Both reduce
performance. Higher accuracy is better.

To test the properties of the attractor, we measure model accuracy as
embeddings reach convergence. We separate the samples _post hoc_ based on the
number of steps performed until completion. We discover that the accuracy is
inversely proportional to the number of steps of execution of the algorithm
(Fig. 5 left). Then, we forcefully modify the number of steps the model
executes (Fig. 5 right). We observe that the performance suffers greatly when
we stop earlier. This suggests that the model learns to be closely aligned
with the ground truth Bellman-Ford algorithm. Specifically, it learns to
always reason about new nodes on the same step that Bellman-Ford does.
Stopping early, therefore, prevents it from completing its computation. Even
more interestingly, performance also suffers when we run for longer than
needed, though to a lesser degree. This is consistent with our previous
observations, and shows that the attractor state is somewhat unstable under
the model update.

##  Appendix C Ablations

###  C.1 Latent Spaces of Triplet-GMPNN

Our analysis of latent spaces was primarily focused on the LinearPGN. In
Figure 6, we show the latent spaces for random ER graphs of Triplet-GMPNN. It
is evident that the structure of latent spaces is similar to that of
LinearPGN. Quantitatively, Triplet-GMPNN achieves slightly larger explained
variance ratio, but qualitatively there is no substantial difference.

Figure 6: Latent space visualisation of Triplet-GMPNN.

In fact, the only aspect of our experiments that holds for LinearPGN but not
for Triplet-GMPNN is the analysis of mispredicts in terms of value
generalisation. Triplet-GMPNN achieves over $97\%$ accuracy on Bellman-Ford,
and manages to generalise well to the out-of-distribution graphs. It learns to
correctly handle large distances even without any changes to the pipeline, and
we believe that this is precisely the reason for its improved performance.

Therefore, our work can be used to analyse other types of GNN processors, and
is not limited to just one architecture.

###  C.2 Choice of Node Aggregation

In Section 3.1, we introduced several approaches to using PCA in order to
represent latent spaces. Common to all of the approaches was the need to
reduce trajectory tensor’s dimension from $4$ to $2$. In order to do so, we
aggregated the trajectories for each node by selecting the largest values. As
a refresher, visualisations of trajectory tensors reduced with $\max$
aggregation are shown on Figure 7.

Figure 7: Visualisation of latent spaces with max aggregation.

We could have done this differently, by using mean or $\min$ instead of
$\max$. However, that would not have impacted our analysis. As an example, in
Figure 8 we recreate Figure 7 with $\min$ aggregation over nodes instead of
with $\max$. There is no difference between the two approaches. Qualitatively,
we observe the same visual structure; quantitatively, expected variance ratios
for both trajectory-wise and step-wise PCA are nearly identical.

Figure 8: Visualisation of latent spaces with $\min$ aggregation.

We visualise latent spaces with mean aggregation as well (Figure 9). Again,
there is no qualitative difference, and quantitatively, taking the mean
increases the explained variance ratio of the trajectory-wise PCA from $63\%$
to $94\%$.

Figure 9: Visualisation of latent spaces with mean aggregation.

###  C.3 Decay Strength

The processor decay has a hyperparameter that determines its strength. We
previously argued that only small amounts of decay are needed in order to
sufficiently change the pipeline.

Now, we look at different choices of the scaling factor and compare their
performance. We also juxtapose decaying to zero, where we merely scale with
the factor, against scaling differences to the mean.

From Table C.3 we observe that even decay with a factor of $0.9$, which
corresponds to reducing embeddings by $10\%$ per step, is large enough to
improve performance. Furthermore, both decaying towards zero, and decaying
towards a mean that is calculated on a per-node basis, perform similarly. We
conclude that mere presence of decay is enough to reshape the processor.

| Decay factor
---|---
Decay type | $1$ | $0.9$ | $0.5$
To Zero |  |  ${\mathbf{98.76\%}}$$\pm 0.10$ | To Mean |

Table 5: Ablation of decay factor on Triplet-GMPNN on Bellman-Ford. Factor of
$1$ means no decay, and smaller values increase decay strength. Standard
deviation estimated with three runs.

###  C.4 Softmax Temperature

Our softmax aggregation also has a hyperparameter called temperature.
Temperature controls how similar softmax is to hard $\max$, and how big an
effect small values have on the result. Temperature is always non-negative.
With temperature zero, softmax aggregation behaves exactly as $\max$. On the
other hand, as temperature approaches $+\infty$, softmax aggregation behaves
as mean. Thus, in order to preserve $\max$-like behaviour, temperature should
be very small. But, in order to be differentiable through all elements,
temperature should be larger than zero.

Temperature | 0 | 0.01 | 0.1 | 1
---|---|---|---|---
Accuracy |  ${\mathbf{98.80\%}}$$\pm 0.14$ |  |  |
Table 6: Ablation of softmax temperature on Triplet-GMPNN on Bellman-Ford.
Larger values of temperature increase impact of small arguments to softmax.
Temperature of zero is equivalent to $\max$. Variance estimated with three
runs.

In Table 6, we display the accuracy of Triplet-GMPNN when trained with
different temperatures. We observe that softmax temperature of $0.01$ performs
the best. Coincidentally, this is the temperature we used in our experiments
above. We also note that large temperature also increases variance between
results.

##  Appendix D Distance Distribution Changes

We now turn to studying the effect of softmax aggregation and processor decay
on the distribution of mispredicts.

Figure 10: Distribution of predicted distances when NAR is trained with decay.
Note that all mispredicts are clustered in a region to the right of $d=1$.
Figure 11: Distribution of predicted distances when NAR is trained with
softmax aggregation. Note the scale on the right plot.

Figures 10 and 11 show the distributions when NAR is trained with decay and
softmax aggregation, respectively. In both cases, we observe that the
distributions are drastically altered.

In the case of decay, we observe a clear separation between the correctly
predicted and the mispredicted nodes. Almost all nodes with predicted distance
under $1$ are predicted correctly, while almost all nodes with predicted
distance above $1.2$ are mispredicted. With this approach, we can know with
high confidence how trustworthy NAR predictions are.

Softmax aggregation does not delineate between true and false predictions like
decay does, but it greatly improves the overall performance. With softmax, NAR
achieves impressive $96.5\%$ on the OOD $p=0.25$ dataset. In Figure 12 we plot
the correlations between the true shortest distances and the ones predicted by
NAR. This way, we can quantify how softmax improves the ability of the model
to handle large unseen values.

Figure 12: Left: Correlation between predicted and true distances for vanilla
NAR. Right: Correlation for NAR with softmax aggregation. Correctly predicted
nodes in blue, mispredicts in orange.

##  Appendix E Bellman-Ford Prediction Visualisation

Here we visualise step-by-step predictions of vanilla Triplet-GMPNN NAR
trained on Bellman-Ford algorithm, and how mispredicts can occur when
comparing similar values.

Node 0 is source and is colored red, while others are gray. For each node, its
id ($0$ to $9$) and predicted distance to source are overlaid on top.

To visualise predictions of node pointers $\pi$, we observe that they define a
spanning tree with minimum distances to source. Therefore, we color each edge
$uv$

  * •

Blue if it belongs to a ground truth spanning tree ($\pi_{\text{true}}(v)=u$
or $\pi_{\text{true}}(u)=v$)

  * •

Green if it belongs to a ground truth spanning tree and is correctly predicted
by NAR ($\pi_{\text{true}}(v)=u$ or $\pi_{\text{true}}(u)=v$ and
$\pi_{\text{pred}}(v)=u$ or $\pi_{\text{pred}}(u)=v$)

  * •

Black if it does not belong to a ground truth spanning tree
($\pi_{\text{true}}(v)\neq u$ and $\pi_{\text{true}}(u)\neq v$)

  * •

Red if it does not belong to a ground truth spanning tree but is mispredicted
as one ($\pi_{\text{true}}(v)\neq u$ and $\pi_{\text{true}}(u)\neq v$ but
$\pi_{\text{pred}}(v)=u$ or $\pi_{\text{pred}}(u)=v$)

The graph is undirected.

After step 1, node 5 is correctly directed to source. None of the other nodes
are reachable yet.

After step 2, node 6 is also correctly directed to source, and after step 3,
so are nodes 4, 7 and 8.

With step 4 complete, all discovered nodes are predicted correctly. One more
step needs to be made to reach node 3.

In this final step, neural network suddenly makes three mispredictions. It
sets $\pi(4)=2$ from previous 6, $\pi(9)=2$ from previous 8, and $\pi(1)=9$
from previous 7. Crucially, in all of those cases, predicted distances to
source between these pairs of paths are similar ($1.69$ vs $1.6$, $1.41$ vs
$1.45$, and $1.79$ vs $1.71$ respectively).

This failure mode is present among all Bellman-Ford mispredicts we observed,
and is the argument for using softmax aggregation outlined in Section 6.

##  Appendix F Latent Space Trajectories - Sequential View

Figure 13: Sequential view of embedding trajectories in Figure 1 left. Figure
14: Sequential view of embedding trajectories in Figure 1 middle.

On these figures we visualise the embeddings by ordering the execution steps
from left to right and tracing the trajectories in red. We see that the first
few steps greatly change the step embeddings, but that trajectories switch to
being near-parallel as they approach the attractor.

##  Appendix G 3D Visualisations

Here, we display high-resolution visualisations of latent spaces, both
trajectory-wise (left) and step-wise (right).

Figure 15: 3D visualisation of random ER graphs. Figure 16: 3D visualisation
of permutation symmetric graphs. Figure 17: 3D visualisation of scale
symmetric graphs. Figure 18: 3D visualisation of reweighting symmetric graphs.
Figure 19: 3D visualisation of reweighting symmetric graphs, with eight
different clusters.