DPUR/TH/70

March, 2021

Quadratic Gravity and Restricted Weyl Symmetry

Arata Kamimura111 E-mail address<EMAIL_ADDRESS>and Ichiro
Oda222 E-mail address<EMAIL_ADDRESS>
Department of Physics, Faculty of Science, University of the Ryukyus,

Nishihara, Okinawa 903-0213, Japan

###### Abstract

We investigate the relationship between quadratic gravity and a restricted
Weyl symmetry where a gauge parameter $\Omega(x)$ of Weyl transformation
satisfies a constraint $\Box\Omega=0$ in a curved space-time. First, we
briefly review a model with a restricted gauge symmetry on the basis of QED
where a $U(1)$ gauge parameter $\theta(x)$ obeys a similar constraint
$\Box\theta=0$ in a flat Minkowski space-time, and explain that the restricted
gauge symmetry removes one on-shell mode of gauge field, which together with
the Feynman gauge leaves only two transverse polarizations as physical states.
Next, it is shown that the restricted Weyl symmetry also eliminates one
component of a dipole field in quadratic gravity around a flat Minkowski
background, leaving only a single scalar state. Finally, we show that the
restricted Weyl symmetry cannot remove any dynamical degrees of freedom in
static background metrics by using the zero-energy theorem of quadratic
gravity. This fact also holds for the Euclidean background metrics without
imposing the static condition.

##  1 Introduction

We have recently watched a renewed interest in scale invariant gravitational
theories since scale symmetries are expected to play a dominant role at high
energies [1]-[9]. For instance, let us consider the era of Big Bang which is
in very high energies. At the Big Bang, it was extremely hot, so the kinetic
energy of particles would be so enormous as to have completely overwhelmed the
particles’ rest mass. In such a situation, the rest mass would be effectively
irrelevant, that is, as good as zero, so that the universe at early times,
could be described in terms of effectively massless particles, for which scale
symmetries naturally emerge. If scale symmetries are exact ones which are only
spontaneously broken as in the gauge symmetries in the standard model (SM) of
particle physics, they might shed some light on various important unsolved
problems such as cosmological constant problem, the origin of particles’ mass
and the gauge hierarchy problem etc.

In this article, we further explore an idea that there might be an
intermediate scale symmetry between local (or Weyl) and global scale
symmetries, which is dubbed a restricted Weyl symmetry in a curved space-time
[10]-[13]. In the restricted Weyl symmetry, a gauge transformation parameter,
which we will call $\it{scale\,factor}$ $\Omega(x)$, is constrained by a
harmonic condition $\Box\Omega=0$ whereas in a conventional or full Weyl
transformation the scale factor is an unconstrained and free function. In the
restricted Weyl symmetry, we are allowed to work with actions with
dimensionless coupling constants as in globally scale invariant gravitational
theories, which should be contrasted to the situation in the full Weyl
symmetry where only the conformal (or Weyl) tensor squared is an invariant
action if we neglect the other fields except for the metric $g_{\mu\nu}$.

Another advantage of exploring theories with restricted symmetries is that we
can understand spontaneous symmetry breakdown of global symmetries [13]. For
instance, as explained in the next section, a QED with a restricted gauge
symmetry has been used to understand why the photon is exactly massless by
studying the zero-mode solutions to the harmonic condition on the gauge
parameter [14]. In fact, it has been already shown that a gravitational theory
with the restricted Weyl symmetry naturally gives rise to spontaneous symmetry
breakdown of a global scale symmetry [13].

The structure of this article is the following: In Section 2, we review a
restricted gauge symmetry in QED which appears after we fix the gauge
invariance by the Lorenz gauge. The restricted gauge symmetry resembles the
restricted Weyl symmetry in the sense that both the symmetries are constrained
by an equation $\Box\Theta=0$ where $\Box$ denotes the d’Alembertian operator
and $\Theta$ stands for a generic gauge transformation parameter although
there is a big difference in that in the former case the operator $\Box$ is
the d’Alembertian operator in a flat Minkowski space-time whereas in the
latter case it is defined in a curved Riemannian space-time.333To distinguish
the difference of the d’Alembertian operators, we will henceforth use
$\Box_{\eta}$ for the Minkowski space-time and $\Box_{g}$ for the Riemannian
space-time, respectively. In Section 3, we introduce quadratic gravity with
the restricted Weyl symmetry and discuss a relation between the quadratic
gravity and conformal gravity on the basis of the BRST formalism. In Section
4, using conformal gravity, we explain how the restricted Weyl symmetry can be
understood in the framework of the BRST quartet mechanism. In Section 5, we
present a proof that the restricted Weyl symmetry in fact kills one dynamical
degree of freedom in the quadratic gravity. In Section 6, the physical meaning
of the constraint in the restricted Weyl symmetry is investigated, for which
the zero-energy theorem of the quadratic gravity [15] gives us information on
the falloff behavior at infinity. The final section is devoted to conclusion.

##  2 Review of restricted gauge symmetry

In order to clarify the problems which we attempt to understand in this
article, it is worthwhile to start with QED which is gauge-fixed by the Lorenz
gauge whose Lagrangian takes the form

$\displaystyle{\cal{L}}=-\frac{1}{4}F_{\mu\nu}^{2}+\bar{\psi}(i\gamma^{\mu}\partial_{\mu}-m)\psi+eA_{\mu}\bar{\psi}\gamma^{\mu}\psi+B\partial_{\mu}A^{\mu}+\frac{\alpha}{2}B^{2},$
(1)

where $A_{\mu},\psi$ and $B$ are respectively the electromagnetic field,
spinor field and Nakanishi-Lautrup field, the field strength is defined as
$F_{\mu\nu}\equiv\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$, and $\alpha$ is
a real number called “gauge parameter” which can be chosen at our disposal.
Even if the $U(1)$ gauge symmetry has been already gauge-fixed by the Lorenz
gauge, this Lagrangian still has a residual symmetry called a restricted
$U(1)$ gauge invariance444This residual gauge transformation was also called
the “c-number” gauge transformation long ago in comparison to the “q-number”
gauge transformation where the gauge parameter is also changed [16, 17]. which
is explicitly given by

$\displaystyle\delta
A_{\mu}=\partial_{\mu}\theta,\quad\delta\psi=ie\theta\psi,\quad\delta\bar{\psi}=-ie\theta\bar{\psi},\quad\delta
B=0,$ (2)

where $\Box_{\eta}\theta=0$ where
$\Box_{\eta}=\eta^{\mu\nu}\partial_{\mu}\partial_{\nu}$ is the d’Alembertian
operator in a flat Minkowski space-time with the flat metric
$\eta_{\mu\nu}=\textrm{diag}(-1,1,1,1)$.

The restricted $U(1)$ gauge symmetry has two physically interesting aspects.
The one is that using this symmetry one can understand why the photon is
exactly massless. Actually, we can construct the generator corresponding to
the restricted $U(1)$ gauge transformation (2)

$\displaystyle Q=\int d^{3}x(\theta\partial_{0}B-\partial_{0}\theta\cdot B).$
(3)

Based on canonical quantization procedure, we can calculate a four-dimensional
commutation relation

$\displaystyle[A_{\mu}(x),B(y)]=i\partial_{\mu}D(x-y),$ (4)

where $D(x)$ is called “invariant D function” which is defined as

$\displaystyle D(x)\equiv-\frac{i}{(2\pi)^{3}}\int
d^{4}k\,\varepsilon(k^{0})\delta(k^{2})e^{ik\cdot x},$ (5)

and obeys the relations

$\displaystyle\Box_{\eta}D(x)=0,\qquad
D(0,\vec{x})=0,\qquad\partial_{0}D(0,\vec{x})=-\delta(\vec{x}).$ (6)

Using Eqs. (3), (4) and (6), we can derive the following vacuum expectation
value:

$\displaystyle\langle 0|[iQ,A_{\mu}(x)]|0\rangle=\partial_{\mu}\theta(x),$ (7)

which implies that $Q$ is necessarily broken spontaneously if $\theta(x)$ is
not a constant. Then, the Goldstone theorem tells us that $A_{\mu}$ contains a
massless state, which can be identified with the photon after a detailed
analysis [14]. In other words, the restricted gauge symmetry provides us with
useful information on spontaneous symmetry breakdown of global
symmetries.555In the detailed analysis, we need to consider zero-mode
solutions to the constraint $\Box_{\eta}\theta=0$, which are explicitly given
by $\theta(x)=a_{\mu}x^{\mu}+b$ where $a_{\mu}$ and $b$ are constants. It is
global symmetries corresponding to global parameters $a_{\mu}$ and $b$ that
spontaneous symmetry breakdown could occur. Precisely speaking, the Ferrari-
Piccaso’s proof is not completely correct, but the claim is true that the
required massless particle must be a vector when the U(1) symmetry
corresponding to the $\theta(x)=b$ transformation remains unbroken.

The other interesting aspect is that we can gauge away one on-shell mode of
$A_{\mu}$ by means of the symmetry (2). With regards to this point, let us
note that since $\Box_{\eta}\theta=0$ is nothing but the Klein-Gordon equation
for a massless real scalar field, a general solution is given by

$\displaystyle\theta(x)=\int\frac{d^{3}k}{\sqrt{(2\pi)^{3}2k_{0}}}[a(k)e^{ikx}+a^{\dagger}(k)e^{-ikx}],$
(8)

where $k_{0}=|\vec{k}|$.

The field equations from the Lagrangian (1) read

$\displaystyle\partial_{\nu}F^{\nu\mu}-\partial^{\mu}B+ej^{\mu}=0,$
$\displaystyle\partial^{\mu}A_{\mu}+\alpha B=0,$
$\displaystyle[i\gamma^{\mu}(\partial_{\mu}-ieA_{\mu})-m]\psi=0,$ (9)

where the electric current $j^{\mu}$ is defined by
$j^{\mu}=\bar{\psi}\gamma^{\mu}\psi$. Using these field equations, we can
derive the following equation

$\displaystyle\Box_{\eta}A^{\mu}-(1-\alpha)\partial^{\mu}B+ej^{\mu}=0.$ (10)

In the absence of the electric current $j^{\mu}=0$ and with the Feynman gauge
$\alpha=1$, the gauge field obeys the same field equation
$\Box_{\eta}A_{\mu}=0$ as in $\Box_{\eta}\theta=0$. This fact implies that in
terms of the restricted gauge symmetry we can gauge away one on-shell mode of
the gauge field, for instance, a temporal component of the gauge field,
$A_{0}=0$, which together with the Feynman gauge condition leaves two physical
massless states for $A_{\mu}$ identified with two transverse polarizations of
the photon.666In the Lorenz gauge, this statement obviously holds, but in the
Feynman gauge it is a bit sutble. Namely, in the Feynman gauge, the NL field
is determined as $B=-\partial^{\mu}A_{\mu}$ and it appears in the physical
Hilbert space, but it is not observed owing to a zero-norm state, leading to
the desired result.

Our purposes in this article is two-fold, one of which is to investigate
whether a similar mechanism for eliminating dynamical degrees of freedom works
even for the restricted Weyl symmetry or not. We will see that the restricted
Weyl symmetry can indeed gauge away one component of the metric tensor around
a flat Minkowski background as in the restricted gauge symmetry. The other
purpose is to understand the physical meaning of the constraint of the
restricted Weyl symmetry. We find that in case of static background metrics
with the timelike Killing vector, the restricted Weyl symmetry reduces to a
global scale symmetry by using the falloff condition at infinity which is
required by the zero-energy theorem of quadratic gravity [15]. It also turns
out that this holds for any Euclidean metrics without the static condition.

##  3 Quadratic gravity with restricted Weyl symmetry

In this section, we consider the most general gravitational theory without
matter fields which is invariant under a restricted Weyl transformation (as
well as a global scale transformation) and explain its relation to conformal
gravity. This gravitational theory, which we will call
${\it{quadratic\,gravity}}$, is described by the Lagrangian777We follow the
conventions and notation of the MTW textbook [18].:

$\displaystyle{\cal{L}}=\sqrt{-g}\,(\alpha\,R^{2}+\beta\,C_{\mu\nu\rho\sigma}^{2}),$
(11)

where $\alpha,\beta$ are dimensionless constants and $C_{\mu\nu\rho\sigma}$
the conformal tensor (or Weyl tensor), and $R$ the scalar curvature. Note that
this Lagrangian is invariant under both a global scale transformation and a
restricted Weyl transformation, but is not so under a local scale (or Weyl)
transformation because of the presence of an $R^{2}$ term. The Gauss-Bonnet
topological invariant enables us to rewrite (11) into the form

$\displaystyle{\cal{L}}=\sqrt{-g}\,(\alpha\,R^{2}+\gamma\,R_{\mu\nu}^{2}),$
(12)

where $\gamma$ is also a dimensionless constant. In this case, both $R^{2}$
and $R_{\mu\nu}^{2}$ terms break the local scale symmetry even if they are
invariant under the restricted Weyl and global scale transformations.

Here the restricted Weyl transformation is defined by the Weyl transformation

$\displaystyle g_{\mu\nu}\rightarrow
g^{\prime}_{\mu\nu}=\Omega^{2}(x)g_{\mu\nu},$ (13)

where the gauge transformation parameter, i.e., scale factor $\Omega(x)$,
obeys a harmonic condition $\Box_{g}\Omega=0$ in a curved space-time
[10]-[13]. Note that in proving the invariance of the Lagrangian (11) under
the restricted Weyl transformation, we need to use the transformation law of
the scalar curvature under (13):

$\displaystyle R\rightarrow
R^{\prime}=\Omega^{-2}(R-6\Omega^{-1}\Box_{g}\Omega).$ (14)

The field equations coming from (11) read

$\displaystyle\alpha\Biggl{(}R_{\mu\nu}R-\frac{1}{4}g_{\mu\nu}R^{2}-\nabla_{\mu}\nabla_{\nu}R+g_{\mu\nu}\Box_{g}R\Biggr{)}-\beta\Biggl{(}\frac{2}{3}R_{\mu\nu}R-\frac{1}{6}g_{\mu\nu}R^{2}$
$\displaystyle+\frac{1}{3}\nabla_{\mu}\nabla_{\nu}R+\frac{1}{6}g_{\mu\nu}\Box_{g}R-\Box_{g}R_{\mu\nu}-2R_{\mu}\,^{\alpha}R_{\alpha\nu}+\frac{1}{2}g_{\mu\nu}R_{\alpha\beta}^{2}\Biggr{)}=0.$
(15)

Taking the trace of Eq. (15) produces a harmonic equation for the scalar
curvature:

$\displaystyle\Box_{g}R=0.$ (16)

Now we wish to clarify a relation between the quadratic gravity (11) and
conformal gravity:

$\displaystyle{\cal{L}}_{c}=\sqrt{-g}\,\beta\,C_{\mu\nu\rho\sigma}^{2},$ (17)

which is the unique locally scale invariant gravitational theory in the
absence of matter fields. In fact, we will see that the Lagrangian (11) can be
obtained from conformal gravity by fixing the local scale symmetry in terms of
a gauge condition $R=0$ up to the FP-ghost term [13].

In order to understand this fact, it is useful to make use of the BRST
formalism. The BRST transformation for the local scale transformation reads

$\displaystyle\delta_{B}g_{\mu\nu}$ $\displaystyle=$ $\displaystyle
2cg_{\mu\nu},\qquad\delta_{B}\sqrt{-g}=4c\sqrt{-g},\qquad\delta_{B}R=-2cR-6\Box_{g}c,$
$\displaystyle\delta_{B}\bar{c}$ $\displaystyle=$ $\displaystyle
iB,\qquad\delta_{B}c=\delta_{B}B=0.$ (18)

Then, the Lagrangian for the gauge condition and the FP ghosts can be obtained
by a standard method [19]:

$\displaystyle{\cal{L}}_{GF+FP}$ $\displaystyle=$
$\displaystyle-i\delta_{B}\left[\sqrt{-g}\bar{c}\left(R+\frac{\alpha}{2}B\right)\right]$
(19) $\displaystyle=$
$\displaystyle\sqrt{-g}\left(\hat{B}R+\frac{\alpha}{2}\hat{B}^{2}-6i\bar{c}\Box_{g}c\right)$
$\displaystyle=$
$\displaystyle\sqrt{-g}\left(-\frac{1}{2\alpha}R^{2}+6ig^{\mu\nu}\partial_{\mu}\bar{c}\partial_{\nu}c\right),$

where we have defined $\hat{B}\equiv B+2i\bar{c}c$ and in the last step we
performed the path integral over the auxiliary field $\hat{B}$ and integration
by parts [13].

Finally, adding (19) to the classical Lagrangian (17), we arrive at a gauge-
fixed and BRST-invariant Lagrangian:

$\displaystyle{\cal{L}}_{q}$ $\displaystyle=$
$\displaystyle\sqrt{-g}\Biggl{(}\beta\,C_{\mu\nu\rho\sigma}^{2}-\frac{1}{2\alpha}R^{2}+6ig^{\mu\nu}\partial_{\mu}\bar{c}\partial_{\nu}c\Biggr{)}.$
(20)

Dropping the last FP-ghost term and rewriting $-\frac{1}{2\alpha}$ as $\alpha$
leads to the Lagrangian (11), which is invariant under the restricted Weyl
transformation (as well as a global scale transformation). In this way, we
have succeeded in deriving the quadratic gravity (11) with the restricted Weyl
symmetry and the global scale symmetry by beginning with the conformal gravity
(17) with the local scale symmetry by taking the gauge condition $R=0$.

Incidentally, it is of interest to point out that the FP-ghost Lagrangian is
also invariant under not the full Weyl transformation but the restricted Weyl
transformation. For this aim, let us first assume that both FP-ghost and FP-
antighost have the Weyl weight $-1$, that is, under the Weyl transformation
they transform as

$\displaystyle c\rightarrow
c^{\prime}=\Omega^{-1}(x)c,\qquad\bar{c}\rightarrow\bar{c}^{\prime}=\Omega^{-1}(x)\bar{c}.$
(21)

Then, we find that under the Weyl transformation the ghost term transforms as

$\displaystyle\sqrt{-g}g^{\mu\nu}\partial_{\mu}\bar{c}\partial_{\nu}c$ (22)
$\displaystyle\rightarrow$
$\displaystyle\sqrt{-g}g^{\mu\nu}\left[\partial_{\mu}\bar{c}\partial_{\nu}c+\Omega^{-1}\nabla_{\mu}\nabla_{\nu}\Omega\cdot\bar{c}c-\nabla_{\mu}(\Omega^{-1}\nabla_{\nu}\Omega\cdot\bar{c}c)\right].$

Hence, the ghost kinetic term is invariant under the restricted Weyl
transformation up to a surface term.

##  4 The BRST quartet mechanism

The derivation of the BRST-invariant Lagrangian (20) from conformal gravity
(17) via the gauge fixing procedure demonstrates that the physical content
between the two theories is exactly the same since Eq. (19) is BRST-exact. In
order to explain how the restricted Weyl symmetry can be understood in the
framework of the BRST formalism, it is useful to check the quantum equivalence
between (17) and (20) through the BRST quartet mechanism [20] in an explicit
manner.

To this aim, let us expand the metric around the flat Minkowski background
metric $\eta_{\mu\nu}$:

$\displaystyle g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu},$ (23)

where $h_{\mu\nu}$ indicates a small fluctuation. Moreover, we use the York
decomposition for the metric fluctuation field $h_{\mu\nu}$ [21]:

$\displaystyle
h_{\mu\nu}=h^{TT}_{\mu\nu}+\partial_{\mu}\xi_{\nu}+\partial_{\nu}\xi_{\mu}+\partial_{\mu}\partial_{\nu}\sigma-\frac{1}{4}\eta_{\mu\nu}\Box_{\eta}\sigma+\frac{1}{4}\eta_{\mu\nu}h,$
(24)

where $h^{TT}_{\mu\nu}$ is both transverse and traceless, and $\xi_{\mu}$ is
transverse:

$\displaystyle\partial^{\mu}h^{TT}_{\mu\nu}=\eta^{\mu\nu}h^{TT}_{\mu\nu}=\partial^{\mu}\xi_{\mu}=0.$
(25)

From the BRST transformation (18), we find the following BRST transformation
for asymptotic fields:

$\displaystyle\delta_{B}s=8c,\qquad\delta_{B}c=0,\qquad\delta_{B}\bar{c}=i\hat{B},\qquad\delta_{B}\hat{B}=0,$
(26)

where $s\equiv h-\Box_{\eta}\sigma$ is invariant under the general coordinate
transformation as shown in the next section.

Next, let us consider the field equations for the above asymptotic fields. The
scalar curvature can be written in terms of $s(x)$ at the linear order as

$\displaystyle R=-\frac{3}{4}\Box_{\eta}s.$ (27)

The field equation (16) then produces the field equations for $s(x)$888Here we
can neglect contributions from the FP-ghost term since they start from the
second order in fields.

$\displaystyle\Box_{\eta}^{2}s=0,$ (28)

which means that $s(x)$ is a dipole field. From Eq. (19), we have the field
equation for the NL field $\hat{B}$:

$\displaystyle\hat{B}=-\frac{1}{\alpha}R.$ (29)

With the help of Eq. (16), Eq. (29) yields

$\displaystyle\Box_{\eta}\hat{B}=0,$ (30)

which implies that the NL field $\hat{B}$ is a simple pole field. As for the
field equations for the FP-ghosts, from Eq. (19) we can easily find that

$\displaystyle\Box_{\eta}c=\Box_{\eta}\bar{c}=0.$ (31)

Thus, the FP-ghosts are also simple pole fields.

The dipole field $s(x)$ can be decomposed into two massless scalar fields.
Actually, the dipole equation is expressed in terms of two coupled equations
by introducing an additional scalar field $\phi(x)$ as

$\displaystyle\Box_{\eta}s=\phi,\qquad\Box_{\eta}\phi=0.$ (32)

Then, by introducing a nonlocal operator ${\cal{D}}$ [20], which works as an
inverse of the d’Alembertian operator $\Box_{\eta}$, defined by

$\displaystyle{\cal{D}}\equiv-\frac{1}{2\Delta}\left(x_{0}\partial_{0}-\frac{1}{2}\right),$
(33)

we can extract a simple pole field $\tilde{s}$ as

$\displaystyle\tilde{s}=s-{\cal{D}}\phi.$ (34)

Using $\tilde{s}$ instead of $s$, the BRST transformation (26) becomes

$\displaystyle\delta_{B}\tilde{s}=8c,\qquad\delta_{B}c=0,\qquad\delta_{B}\bar{c}=i\hat{B},\qquad\delta_{B}\hat{B}=0,$
(35)

and all the asymptotic fields satisfy the simple pole field equation:

$\displaystyle\Box_{\eta}\tilde{s}=\Box_{\eta}c=\Box_{\eta}\bar{c}=\Box_{\eta}\hat{B}=0,$
(36)

which clearly indicates that the set of the asymptotic fields
$\\{\tilde{s},c,\bar{c},\hat{B}\\}$ belongs to a BRST quartet so that they
drop from physical Hilbert state by the Kugo-Ojima’s subsidiary condition
[20].999The physical Hilbert space ${\cal{H}}_{phys}$ is defined as the
quotient space of ${\cal{V}}_{phys}$ defined by the subsidiary condition
$Q_{B}{\cal{V}}_{phys}=0$ ($Q_{B}$ is the BRST charge corresponding to the
Weyl symmetry) with respect to its zero-norm subspace ${\cal{V}}_{0}$, i.e.,
${\cal{H}}_{phys}=\frac{{\cal{V}}_{phys}}{{\cal{V}}_{0}}$. The BRST quartet
belongs to the zero-norm subspace ${\cal{V}}_{0}$. Here it is worth noticing
that using Eqs. (27), (29) and (32), the NL field $\hat{B}$ is proportional to
the scalar field $\phi$:

$\displaystyle\hat{B}=-\frac{1}{\alpha}R=\frac{3}{4\alpha}\Box_{\eta}s=\frac{3}{4\alpha}\phi.$
(37)

Since both $\tilde{s}(x)$ and $\phi(x)$ belong to the BRST quartet, the dipole
field $s(x)$ completely drops from the physical Hilbert space by the Kugo-
Ojima’s subsidiary condition as expected.

An important lesson obtained by means of the BRST quartet mechanism is that it
is essential to decompose the dipole field into two scalar fields as in Eq.
(32), both of which can be eliminated from the physical Hilbert space by the
Kugo-Ojima’s subsidiary condition [20]. In the quadratic gravity with the
restricted Weyl symmetry discussed in the next section, it turns out that it
is also important to make use of this decomposition even if there is no BRST
symmetry.

##  5 Restricted Weyl symmetry

In this section, we will prove that the restricted Weyl symmetry can remove
one dynamical degree of freedom included in the dipole field $s(x)$.

To do so, let us begin by the Lagrangian (11) and rewrite it in terms of the
York variables up to quadratic terms in the fields:

$\displaystyle{\cal{L}}=\frac{9}{16}\alpha\,s\Box_{\eta}^{2}s+\frac{1}{2}\beta\,h^{TT\mu\nu}\Box_{\eta}^{2}h^{TT}_{\mu\nu}.$
(38)

The fields equations are therefore given by two dipole equations:

$\displaystyle\Box_{\eta}^{2}s=0,\qquad\Box_{\eta}^{2}h^{TT}_{\mu\nu}=0.$ (39)

Under the infinitesimal general coordinate transformation, $h_{\mu\nu}$
transforms as

$\displaystyle h_{\mu\nu}\rightarrow
h_{\mu\nu}^{\prime}=h_{\mu\nu}+\partial_{\mu}\varepsilon_{\nu}+\partial_{\nu}\varepsilon_{\mu},$
(40)

where $\varepsilon_{\mu}$ is the infinitesimal gauge parameter. Here let us
divide $\varepsilon_{\mu}$ into transverse and longitudinal modes:

$\displaystyle\varepsilon_{\mu}=\varepsilon_{\mu}^{T}+\partial_{\mu}\varepsilon,$
(41)

where $\partial^{\mu}\varepsilon_{\mu}^{T}=0$. Then, each component in the
York decomposition (24) transforms like

$\displaystyle h_{\mu\nu}^{\prime\,TT}$ $\displaystyle=$ $\displaystyle
h_{\mu\nu}^{TT},\qquad\xi_{\mu}^{\prime}=\xi_{\mu}+\varepsilon_{\mu}^{T},$
$\displaystyle\sigma^{\prime}$ $\displaystyle=$
$\displaystyle\sigma+2\varepsilon,\qquad
h^{\prime}=h+2\Box_{\eta}\varepsilon.$ (42)

Hence, it is obvious that the Lagrangian (38) is manifestly invariant under
the general coordinate transformation owing to $s^{\prime}=s$.

Under the infinitesimal restricted Weyl transformation (13), we have

$\displaystyle h_{\mu\nu}^{\prime\,TT}=h_{\mu\nu}^{TT},\qquad
s^{\prime}=s+8\Lambda,$ (43)

where we have put the gauge parameter to be its infinitesimal form,
$\Omega(x)=e^{\Lambda(x)}$ ($|\Lambda|\ll 1$), and $\Lambda$ also obeys a
harmonic condition $\Box_{\eta}\Lambda=0$.

Now let us focus our attention to the dipole equation for $s(x)$ in (39). As
mentioned in Eq. (32), the dipole field $s(x)$ can be expressed in terms of
two coupled equations by introducing a new scalar field $\phi(x)$. By using
the Green’s function $G(x,y)$ satisfying $\Box_{\eta}G(x,y)=\delta^{4}(x-y)$,
a general solution for the former equation in (39) is provided by101010This
equation is in essence equivalent to Eq. (34).

$\displaystyle s(x)=\int d^{4}y\,G(x,y)\phi(y)+s_{0}(x),$ (44)

where $s_{0}(x)$ satisfies the d’Alembert’s equation $\Box_{\eta}s_{0}=0$. The
infinitesimal restricted Weyl transformation for the two scalar fields $\phi$
and $s_{0}$ becomes

$\displaystyle\phi^{\prime}=\phi,\qquad s^{\prime}_{0}=s_{0}+8\Lambda.$ (45)

It is then possible to take a gauge condition for the restricted Weyl symmetry

$\displaystyle s_{0}(x)=0,$ (46)

since $s_{0}(x)$ is a gauge-variant field and both $s_{0}(x)$ and $\Lambda(x)$
are harmonic functions. As a result, the $R^{2}$ term in the quadratic gravity
propagates just a single spin-0 state associated with the massless scalar
field $\phi(x)$ around a flat Minkowski background, and the field $s(x)$ is
described by $\phi(x)$ as

$\displaystyle s(x)=\int d^{4}y\,G(x,y)\phi(y).$ (47)

Note that unlike the case of conformal gravity in the previous section, in the
quadratic gravity we have no BRST transformation for the Weyl transformation,
so the massless scalar field $\phi(x)$ in $s(x)$ remains as a propagating
state while the scalar field $s_{0}(x)$ is eliminated by the restricted Weyl
symmetry.

##  6 Restricted Weyl symmetry and zero-energy theorem of quadratic gravity

The d’Alembertian operator in a curved space-time is a geometrical operator in
the sense that its definition needs information on geometry through the metric
tensor. It is therefore natural to conjecture that the restricted Weyl
symmetry might not be able to remove dynamical degrees of freedom for certain
background geometries. Indeed, it is well known that there are only static,
asymptotically constant solutions to the massless Klein-Gordon equation for a
scalar field if the first derivative of the scalar field approaches zero
sufficiently fast at infinity [22]-[25]. In this section, we wish to show
explicitly that our conjecture is valid if the space-time is static in the
sense that there is the timelike Killing vector field. It is of interest to
note that the falloff conditions necessary for our derivation stems from a
proof of the zero-energy theorem of the quadratic gravity [15]. Moreover, we
will show that this fact holds for the Euclidean metrics as well although the
space-time is not static.

We would like to consider a geometry where the space-time is static in the
sense that ${\cal{L}}_{t}g_{\mu\nu}=0$ for which $t_{\mu}$ is the timelike
Killing vector field. For clarity, let us consider the following static line
element:

$\displaystyle ds^{2}\equiv
g_{\mu\nu}dx^{\mu}dx^{\nu}=-f^{2}(x^{k})dt^{2}+g_{ij}(x^{k})dx^{i}dx^{j},$
(48)

where $f$ and $g_{ij}$ are functions of only spatial coordinates $x^{k}$ and
$i,j,k,\cdots=1,2,3$.

With these assumptions, we multiply the constraint $\Box_{g}\Omega=0$ by
$f\Omega$ and integrate over spatial coordinates:

$\displaystyle\int d^{3}x\sqrt{{}^{(3)}g}\,f\Omega\Box_{g}\Omega=0,$ (49)

where ${}^{(3)}g\equiv\det g_{ij}$. It turns out that $\Box_{g}\Omega=0$ can
be rewritten into the form111111Here we have considered that the scale factor
is also a function of spatial coordinates, i.e., $\Omega=\Omega(x^{k})$ since
a locally scaled metric $g^{\prime}_{\mu\nu}=\Omega^{2}(x)g_{\mu\nu}$ must be
static as well.

$\displaystyle\Box_{g}\Omega=\frac{1}{f}g^{ij}\partial_{i}f\partial_{j}\Omega+g^{ij}D_{i}D_{j}\Omega,$
(50)

where $D_{i}$ is the covariant derivative with respect to the spatial metric
$g_{ij}$. Substituting Eq. (50) into Eq. (49) leads to

$\displaystyle\int d^{3}x\sqrt{{}^{(3)}g}\left[D^{i}(f\Omega
D_{i}\Omega)-f(D_{i}\Omega)^{2}\right]=0.$ (51)

If $D_{i}\Omega\rightarrow 0$ sufficiently fast at infinity, then the surface
term makes no contribution, and the non-positivity of the remaining term
produces $D_{i}\Omega=0$, which means that the scale factor $\Omega$ is a
constant, thereby the restricted Weyl symmetry reducing to the global scale
symmetry.

Next, let us show that $D_{i}\Omega\rightarrow 0$ sufficiently fast at
infinity on the basis of the zero-energy theorem of quadratic gravity [15]. In
the quadratic gravity, the definition of asymptotic flatness permitting a
proof of the zero-energy theorem is given by [15]

$\displaystyle
g_{ij}\rightarrow\delta_{ij}+{\cal{O}}\left(r^{-1-\varepsilon}\right),$ (52)

where $\varepsilon$ is a small positive quantity. Since the zero-energy
theorem should hold for both a metric $g_{\mu\nu}$ and its locally scaled
metric
$g^{\prime}_{\mu\nu}=\Omega^{2}(x)g_{\mu\nu}=e^{2\Lambda(x)}g_{\mu\nu}$, the
metric $g^{\prime}_{\mu\nu}$ must obey the same falloff behavior as
$g_{\mu\nu}$. Using an infinitesimal gauge parameter $\Lambda(x)$, we obtain

$\displaystyle g^{\prime}_{ij}$ $\displaystyle\simeq$
$\displaystyle(1+2\Lambda)\left[\delta_{ij}+{\cal{O}}\left(r^{-1-\varepsilon}\right)\right]$
(53) $\displaystyle=$
$\displaystyle\delta_{ij}+2\Lambda\delta_{ij}+{\cal{O}}\left(r^{-1-\varepsilon}\right).$

Thus, $\Lambda(x)$ must have at least the asymptotic behavior

$\displaystyle\Lambda={\cal{O}}\left(r^{-1-\varepsilon}\right).$ (54)

Then, we have

$\displaystyle D_{i}\Omega={\cal{O}}\left(r^{-2-\varepsilon}\right).$ (55)

Hence, at infinity $D_{i}\Omega$ goes to zero so rapidly that the surface term
in (51) vanishes.

Of course, the argument done thus far cannot be applied to the case where
there is no timelike Killing vector. However, in case of the positive semi-
definite Euclidean metrics, we can show a similar result without imposing the
static condition as follows: Start with the constraint $\Box_{g}\Omega=0$ of
the restricted Weyl transformation, multiply it by $\Omega$, and then
integrate over space-time coordinates with a four dimensional Euclidean
measure $\sqrt{g}$ as

$\displaystyle 0=\int d^{4}x\sqrt{g}\,\Omega\Box_{g}\Omega=\int
d^{4}x\sqrt{g}\left[\nabla^{\mu}(\Omega\nabla_{\mu}\Omega)-(\nabla_{\mu}\Omega)^{2}\right].$
(56)

If $\nabla_{\mu}\Omega\rightarrow 0$ sufficiently fast at infinity, then the
surface term can be dropped, and consequently we obtain

$\displaystyle\nabla_{\mu}\Omega=0.$ (57)

This equation implies that $\Omega$ is a constant so that the restricted Weyl
symmetry reduces to the global scale symmetry. This fact might be understood
from the simplest example, the flat Euclidean metric
$\eta_{\mu\nu}=\delta_{\mu\nu}\equiv\textrm{diag}(+1,+1,+1,+1)$. In this
specific case, the momenta $k_{\mu}$ identically vanish since $k_{\mu}^{2}=0$
simply means $k_{\mu}=0$ so that we cannot obtain an infinite number of
oscillator solutions as in Eq. (8) but we have only the constant solution
$\Omega=b$ with $b$ being a constant.121212Note that $\Omega=a_{\mu}x^{\mu}$
($a_{\mu}$ is a constant vector) is not a solution since it does not satisfy
the condition $\nabla_{\mu}\Omega\rightarrow 0$ at infinity. However, in the
flat Euclidean space-time $g_{\mu\nu}=\delta_{\mu\nu}$, this falloff condition
might be too strong since $\Omega=a_{\mu}x^{\mu}$ corresponds to the special
conformal group.

From the viewpoint of the conventional QFT, it might appear to be strange at
first sight that the possibility of eliminating dynamical degrees of freedom
by using the restricted Weyl symmetry depends on the property of the
background metrics at infinity. In order to account for this issue, let us
suppose that the infinitesimal gauge parameter $\Lambda$ has the falloff
behavior like

$\displaystyle\Lambda={\cal{O}}\left(r^{+1-\varepsilon}\right),$ (58)

instead of Eq. (54). Note that with this gauge parameter $\Lambda$, the
restricted Weyl transformation becomes ill-defined at infinity owing to its
divergence, so we cannot require the invariance under such a singular gauge
transformation to define the Hilbert space in QFT. To put it differently,
using such the singular gauge transformation, we cannot eliminate any
dynamical degrees of freedom existing at infinity. Since gauge symmetries
should in general remove dynamical degrees of freedom over the whole space-
time, the singular gauge transfomations are not allowed as the genuine gauge
symmetries.131313Of course, it is possible to consider asymptotic gauge
symmetries holding only at infinity.

##  7 Conclusion

In this article, we have studied the quadratic gravity, which is invariant
under the restricted Weyl transformation as well as the global scale
transformation. It has been explicitly shown that the dipole field $s(x)$,
which is involved in the $R^{2}$ term, has only a single scalar state from the
restricted Weyl symmetry. This mechanism of the elimination of a dynamical
degree of freedom in the quadratic gravity is very similar to that of the QED
with the restricted gauge symmetry. Since both the theories also share a
similar mechansim of spontaneous symmetry breakdown of global symmetries, the
quadratic gravity might be regarded as a gravitational analog of the QED which
is gauged-fixed by the Lorenz gauge.

Moreover, we have investigated a relation between the quadratic gravity and
conformal gravity, which is the unique gravitational theory with the local
scale symmetry, by fixing the Weyl symmetry by the gauge condition $R=0$. In
case of the absence of the FP-ghost term, the two theories have the same
Lagrangian.

We also have investigated the physical meaning of the constraint
$\Box_{g}\Omega=0$ of the restricted Weyl transformation. What we have shown
is that when the background metrics are static, the restricted Weyl
transformation reduces to the global scale transformation, thereby making it
for the restricted Weyl symmetry impossible to remove a dynamical degree of
freedom. This observation has been also generalized to the case where the
space-time metric has the Euclidean signature. Of course, the scale factor is
in general a function of all the space-time coordinates $x^{\mu}$ and the
space-time metric does not have the Euclidean signature but specific cases
such as the instantons, so the restricted Weyl symmetry can indeed eliminate
one degree of freedom except for static space-times. Nevertheless, it is of
interest to notice that the constraint $\Box_{g}\Omega=0$ provides us with a
different physical meaning depending on the property of the space-time.

Acknowledgements

We would like to thank K. Taniguchi and K. Uryu for useful discussions, and in
particular, T. Kugo for valuable discussions and reading of this manuscript.

## References

  * [1] M. Shaposhnikov and D. Zenhausern, “Quantum Scale Invariance, Cosmological Constant and Hierarchy Problem”, Phys. Lett. B 671 (2009) 162.
  * [2] I. Oda, “Higgs Mechanism in Scale-Invariant Gravity”, Adv. Stud. Theor. Phys. 8 (2014) 215\.
  * [3] L. Alvarez-Gaume, A. Kehagias, C. Kounnas and and D. Lüst, “Aspects of Quadratic Gravity”, Fortschr. Phys. 64 (2016) 176.
  * [4] I. Oda, “Classical Weyl Transverse Gravity”, Eur. Phys. J. C 77 (2017) 284.
  * [5] A. Salvio, “Quadratic Gravity”, Front. in Phys. 6 (2018) 77.
  * [6] D. M. Ghilencea, “Spontaneous Breaking of Weyl Quadratic Gravity to Einstein Action and Higgs Potential”, JHEP 1903 (2019) 049.
  * [7] I. Oda, “Planck Scale from Broken Local Conformal Invariance in Weyl Geometry”, PoS CORFU2019 (2020) 070.
  * [8] T. Kugo, “Necessity and Insufficiency of Scale Invariance for Solving Cosmological Constant Problem”, PoS CORFU2019 (2020) 071.
  * [9] I. Oda, “Higgs Potential from Weyl Conformal Gravity”, Mod. Phys. Lett. A 35 (2020) 2050304.
  * [10] A. Edery and Y. Nakayama, “Restricted Weyl Invariance in Four-dimensional Curved Spacetime”, Phys. Rev. D 90 (2014) 043007.
  * [11] A. Edery and Y. Nakayama, “Generating Einstein Gravity, Cosmological Constant and Higgs Mass from Restricted Weyl Invariance”, Mod. Phys. Lett. A 30 (2015) 1550152.
  * [12] A. Edery and Y. Nakayama, “Critical Gravity from Four Dimensional Scale Invariant Gravity”, JHEP 1911 (2019) 169.
  * [13] I. Oda, “Restricted Weyl Symmetry”, Phys. Rev. D 102 (2020) 045008.
  * [14] R. Ferrari and L. E. Picasso, “Spontaneous Breakdown in Quantum Electrodynamics”, Nucl. Phys. B 31 (1971) 316.
  * [15] D. G. Boulware, G. T. Horowitz and A. Strominger, “Zero Energy Theorem for Scale Invariant Gravity”, Phys. Rev. Lett. 50 (1983) 1726.
  * [16] N. Nakanishi and I. Ojima, “Covariant Operator Formalism of Gauge Theories and Quantum Gravity”, World Scientific, Singapore, 1990.
  * [17] K. Yokoyama, “Canonical Quantum Electrodynamics with Invariant One-parameter Gauge Families”, Prog. Theor. Phys. 51 (1974) 1956.
  * [18] C. W. Misner, K. S. Thorne and J. A. Wheeler, “Gravitation”, W H Freeman and Co (Sd), 1973.
  * [19] T. Kugo and S. Uehara, “General Procedure of Gauge Fixing Based on BRS Invariance Principle”, Nucl. Phys. B 197 (1982) 378.
  * [20] T. Kugo and I. Ojima, “Local Covariant Operator Formalism of Non-Abelian Gauge Theories and Quark Confinement Problem”, Prog. Theor. Phys. Suppl. 66 (1979) 1.
  * [21] J. W. York, “Conformally Invariant Orthogonal Decomposition of Symmetric Tensors on Riemannian Manifolds and the Initial Value Problem of General Relativity”, J. Math. Phys. 14 (1973) 456.
  * [22] H.-J. Schmidt, “On Static Asymptotically Flat Solutions of Fourth Order Gravitational Field Equations”, Annalen Phys. 499 (1987) 361.
  * [23] W. Nelson, “Static Solutions for Fourth Order Gravity”, Phys. Rev. D 82 (2010) 104026.
  * [24] H. Lü, A. Perkins, C. N. Pope and K. S. Stelle, “Black Holes in Higher Derivative Gravity”, Phys. Rev. Lett. 114 (2015) 171601.
  * [25] A. Kehagias, C. Kounnas, D. D. Lüst and A. Riotto, “Black Hole Solutions in $R^{2}$ Gravity”, JHEP 05 (2015) 143.


# Fixed Point Networks:
Implicit Depth Models with Jacobian-Free Backprop

Samy Wu Fung, Howard Heaton11footnotemark: 1, Qiuwei Li, Daniel McKenzie,
Stanley Osher, Wotao Yin
Department of Mathematics
University of California Los Angeles
Equal contribution

###### Abstract

A growing trend in deep learning replaces fixed depth models by approximations
of the limit as network depth approaches infinity. This approach uses a
portion of network weights to prescribe behavior by defining a limit
condition. This makes network depth implicit, varying based on the provided
data and an error tolerance. Moreover, existing implicit models can be
implemented and trained with fixed memory costs in exchange for additional
computational costs. In particular, backpropagation through implicit depth
models requires solving a Jacobian-based equation arising from the implicit
function theorem. We propose fixed point networks (FPNs), a simple setup for
implicit depth learning that guarantees convergence of forward propagation to
a unique limit defined by network weights and input data. Our key contribution
is to provide a new Jacobian-free backpropagation (JFB) scheme that
circumvents the need to solve Jacobian-based equations while maintaining fixed
memory costs. This makes FPNs much cheaper to train and easy to implement. Our
numerical examples yield state of the art classification results for implicit
depth models and outperform corresponding explicit models.111Codes are
available on Github: github.com/howardheaton/fixed_point_networks

##  1 Introduction

Figure 1: Explicit networks act by computing $S_{\Theta}\circ Q_{\Theta}$.
Implicit networks add an equilibrium condition using $R_{\Theta}$, which is
applied repeatedly to update a latent variable $u^{k}$ that estimates
$u^{\star}=R_{\Theta}(u^{\star};Q_{\Theta}(d))$.

A growing paradigm shift has emerged from explicit to implicit neural networks
[winston2020monotone, bai2019deep, bai2020multiscale, chen2018neural,
ghaoui2019implicit, dupont2019ANODEs, jeon2021differentiable,
zhang2020implicitly, lawrence2020almost, revay2020contracting,
look2020differentiable, gould2019deep]. In the entirely explicit setting, a
neural network prescribes a series of computations that map input data $d$ to
an inference $y$. Networks can also explicitly leverage the assumption that
high dimensional signals typically admit low dimensional representations in
some latent space [van2008visualizing, osher2017low, peyre2009manifold,
elad2010role, udell2019big]. This may be done by designing the network to
first map data to a latent space via a mapping $Q_{\Theta}$ and then apply a
second mapping $S_{\Theta}$ to map the latent variable to the inference. Thus,
an explicit network ${\cal E}_{\Theta}$ may take the compositional form

${\cal E}_{\Theta}(d)=S_{\Theta}(Q_{\Theta}(d)),$ (1)

which is illustrated by the red arrows in Figure 1. The mapping $S_{\Theta}$
may perform a number of computations to manipulate the latent variable $u$
before mapping it to the inference. Inspired by this, implicit networks arose
from the subsequent inquiry:

How can we model the limit as the latent space portion of network depth
approaches infinity?

The limit of latent space computations may be equivalently viewed as imposing
some form of constraint on the latent space portion of the network. The
repeated computations iteratively refine the current latent space variable
until this condition, dependent on the input data $d$ and the weights
$\Theta$, is satisfied. To formulate this limit condition, we introduce a
third mapping $R_{\Theta}$. Given data $d$, $R_{\Theta}(\cdot\ ;Q(d))$ forms
an operator on the latent space. The latent space condition is expressed using
a fixed point equation, i.e. , we define an implicit fixed point network (FPN)
${\cal M}_{\Theta}$ by

${\cal M}_{\Theta}(d)\triangleq S_{\Theta}(u^{\star})\ \ \ \mbox{where}\ \ \
u^{\star}=R_{\Theta}(u^{\star},Q_{\Theta}(d)).$ (2)

Although the definition in (2) may be well-posed (with appropriate
assumptions), it does not provide any indication about how to evaluate the
network ${\cal M}_{\Theta}$. This matter is discussed in Section 2 below.

#### Why Implicit?

Implicit networks, and in particular FPNs, generalize the transition from
residual neural networks (ResNets) [he2016deep] to continuous depth models
[haber2017stable, chen2018neural, bai2019deep, winston2020monotone]. Implicit
networks do not need to store intermediate quantities of the forward pass of a
deep neural network for backpropagation. Instead, implicit models solve a
Jacobian-based equation [winston2020monotone, bai2019deep,
jeon2021differentiable] or adjoint equation [chen2018neural] arising from the
implicit function theorem. Consequently, implicit models enable model training
at constant memory cost – relieving a major bottleneck of training deep
models. Moreover, infinite/implicit depth models are at least as expressible
as explicit models [jeon2021differentiable, ghaoui2019implicit,
frecon2021bilevel]. This can also be observed by setting $R_{\Theta}$ in (2)
to simply return $Q_{\Theta}$; in this case, the implicit network ${\cal
M}_{\Theta}$ reduces to the explicit network ${\cal E}_{\Theta}$ in (1).
Similar arguments are presented in [bai2019deep, Theorem 3] to argue for the
expressiveness of implicit depth networks. Universal approximation properties
of implicit models then follow immediately from such properties of
conventional deep neural networks (e.g. see [csaji2001approximation,
lu2017expressive, kidger2020universal]).

#### Contribution

We first present a new class of neural networks, FPNs, that may be applied to
several types of machine learning problems. FPNs are based on a fixed point
problem formulation that ensures its forward propagation converges to a unique
limit. Our primary contribution is to propose a new and simple Jacobian-free
backpropagation (JFB) technique for training implicit depth networks that
avoids solving the Jacobian-based equation arising from the implicit function
theorem [krantz2012implicit]. Instead, our scheme backpropagates by omitting
the Jacobian term, resulting in a form of preconditioned gradient. This allows
us to maintain the low-memory advantage of implicit depth models _and_ avoid
the extra computational costs that come with solving the Jacobian-based
equation.

##  2 Fixed Point Network Formulation

Herein we formally cover our implicit depth formulation. All terms presented
in this section are provided in a general context, which is made concrete in
subsequent sections for each application. We include a subscript $\Theta$ on
various terms to emphasize the indicated mapping will ultimately be
parameterized in terms of tunable weights222We use the same subscript for all
terms, noting each operator typically depends on a portion of the weights.
$\Theta$. At the highest level, we are interested in constructing a neural
network ${\cal M}_{\Theta}:{\cal D}\rightarrow{\cal Y}$ that maps from a data
space333Each space in this work is assumed to be a real-valued finite
dimensional Hilbert space (e.g. $\mathbb{R}^{n}$) endowed with a product
$\left<\cdot,\cdot\right>$ and norm $\|\cdot\|$. It will be clear from context
which space is being used. ${\cal D}$ and into an inference space ${\cal Y}$.
The implicit portion of the computations occur in a latent space ${\cal U}$.
We define a mapping $T_{\Theta}:{\cal U}\times{\cal D}\rightarrow{\cal U}$ and
refer to this as the network operator, noting that fixing the input data
yields an operator $T_{\Theta}(\cdot\ ;d)$ on the latent space ${\cal U}$.
Provided input data $d$, our aim is to find the unique fixed point
$u^{\star}_{d}$ of $T_{\Theta}(\cdot\ ;d)$ and then map $u_{d}^{\star}$ to the
inference space ${\cal Y}$ via a mapping $S_{\Theta}:{\cal U}\rightarrow{\cal
Y}$. That is, we define

${\cal M}_{\Theta}(d)\triangleq S_{\Theta}(u_{d}^{\star})\ \ \ \mbox{where}\ \
\ u_{d}^{\star}=T_{\Theta}(u_{d}^{\star};d).$ (3)

Algorithm 1 Fixed Point Network (Abstract Form)

1: ${\cal M}_{\Theta}(d):$ |  $\vartriangleleft$ Input data is $d$
---|---

2: $u^{1}\leftarrow\tilde{u}$ |  $\vartriangleleft$ Initialize inference
---|---

3: while $\|u-T_{\Theta}(u;d)\|>\varepsilon$ |  $\vartriangleleft$ Loop to fixed point
---|---

4: $u^{k+1}\leftarrow T_{\Theta}(u^{k};d)$ |  $\vartriangleleft$ Refine latent variable
---|---

5: $k\leftarrow k+1$ |  $\vartriangleleft$ Increment counter
---|---

6: return $S_{\Theta}(u^{k})$ |  $\vartriangleleft$ Output inference
---|---

For simplicity and computational efficiency, we assume $T_{\Theta}$ is defined
as a composition so that it relates to the mapping $R_{\Theta}:{\cal
U}\times{\cal U}\rightarrow{\cal U}$ in (2) by

$T_{\Theta}(u;d)\triangleq R_{\Theta}(u,Q_{\Theta}(d)).$ (4)

Implementation considerations for (4) are discussed below. Below we also
introduce assumptions on $T_{\Theta}$ that yield sufficient conditions to use
the simple procedure in Algorithm 1 to approximate ${\cal M}_{\Theta}(d)$. In
this algorithm, the latent variable initialization $\tilde{u}$ can be anything
(e.g. the zero vector) since convergence is established independently of this
choice. The inequality in Step 3 gives a fixed point residual condition that
measures convergence. Step 4 implements a fixed point update. The estimate of
the inference ${\cal M}_{\Theta}(d)$ is computed by applying $S_{\Theta}$ to
the latent variable $u$ in Step 5. In summary, Algorithm 1 describes the
implicit network (blue path) in Figure 1.

Figure 2: Diagram of a possible architecture for network operator $T_{\Theta}$
(in large rectangle). Data $d$ and latent $u$ variables are processed in two
streams by nonlinearities (denoted by $\sigma$) and affine mappings (denoted
by rectangles). These streams merge into a final stream that may also contain
transformations. Light gray and blue affine maps are $\ell$-Lipschitz and
1-Lipschitz, respectively. The mapping $Q_{\Theta}$ from data space to latent
space is enclosed by the red rectangle.

#### Interpretation

Multiple views of FPNs yield useful intuitions. Figure 1 provides the core
diagram to keep in mind. The network repeatedly performs an operation,
iteratively refining its “understanding” before outputting the final
inference. For example, in image classification image data (encoded as $d$)
are transformed into a latent space variable (encoded as $u$) in a “thought
space” or “feature space.” The repeated application of $T_{\Theta}$
effectively acts as looking over parts of the image again and again until
everything “sinks in.” The inference space is the translation of
thoughts/features into classifications. Similarly, when analyzing sentiment of
a sentence in natural language processing, FPNs read a sentence repeatedly
until its understanding of the sentiment “makes sense.” Upon the final reading
through the sentence, the network’s impression of the sentiment (in thought
space) is consistent. Many instances require a single reading while more
complex sentences may require rereads (e.g. see [bai2019deep] for natural
language processing experiments related to this).

A more classic viewpoint is to think of each fixed point as a solution to an
optimization problem that models the learning task (as is typical in inverse
problems). Network weights $\Theta$ are tuned to make the optimizer closely
resemble the solution to the learning task. In this case, the network operator
$T_{\Theta}$ behaves as an algorithmic update for iteratively solving the
optimization problem (see Algorithm 1).

#### Convergence

The guarantee of finitely many loops in Steps 3 and 4 of Algorithm 1 is a
direct consequence of a classic functional analysis result by Banach
[banach1922operations]. This same approach is used by several implicit models
[ghaoui2019implicit, winston2020monotone, jeon2021differentiable]. A variation
of this theorem for our setting is presented below.

######  Assumption 2.1.

The mapping $T_{\Theta}$ is $L$-Lipschitz with respect to its inputs $(u,d)$,
i.e. ,

$\|T_{\Theta}(u_{1};\ d_{1})-T_{\Theta}(u_{2};\ d_{2})\|\leq
L\|(u_{1},d_{1})-(u_{2},d_{2})\|,\ \ \ \mbox{for all
$(u_{1},d_{1}),(u_{2},d_{2})\in{\cal U}\times{\cal D}$}.$ (5)

Holding $d$ fixed, the operator $T_{\Theta}(\cdot;d)$ is a contraction, i.e.
there exists $\gamma\in[0,1)$ such that

$\|T_{\Theta}(u_{1};\ d)-T_{\Theta}(u_{2};\ d)\|\leq\gamma\|u_{1}-u_{2}\|,\ \
\ \mbox{for all $u_{1},u_{2}\in{\cal U}$.}$ (6)

######  Remark 2.1.

The $L$-Lipschitz condition on $T_{\Theta}$ is used since recent works show
Lipschitz continuity with respect to inputs improves generalization
[sokolic2017robust, Gouk2021Regularisation, finlay2018lipschitz] and
adversarial robustness [cisse2017parseval, anil2019sorting].

######  Theorem 2.1.

(Banach) For any $u^{1}\in{\cal U}$, if the sequence $\\{u^{k}\\}$ is
generated via the update relation

$u^{k+1}=T_{\Theta}(u^{k};\ d),\ \ \ \mbox{for all $k\in{\mathbb{N}}$},$ (7)

and if Assumption 2.1 holds, then $\\{u^{k}\\}$ converges to the unique fixed
point $u_{d}^{\star}$ with $R$-linear rate.

#### Explicit $\to$ Implicit

As stated in the introduction, an implicit network is at least as expressive
as an explicit network. The most straightforward way to transform an explicit
network into an implicit network is to convert the hidden portion of the
network into a fixed point operator. Given any explicit network (e.g. ResNet56
[he2016deep]), we may transform it into an implicit network in three steps.
First, extract the hidden portion of the network and denote it by
$T_{\Theta}$. Next, make $T_{\Theta}$ a $\gamma$-contraction in $u$ so that
Assumption 2.1 holds (e.g. by replacing residual updates with convex
combinations and by bounding singular values of each linear mapping). Lastly,
denote the input layer which maps to the latent space by $Q_{\Theta}$, and the
final layer (e.g. the classification layer) by $S_{\Theta}$. Completing these
steps, our explicit network has transformed into the implicit FPN ${\cal
M}_{\Theta}$ in Equation (2).

More intricate ways to go from explicit to implicit involve making the hidden
portion of the given explicit network into _multiple_ fixed point conditions.
For instance, if a network transforms the features across different
resolutions, we may set a fixed point condition at every resolution
[bai2020multiscale].

##  3 Backpropagation

We present a new and simple way to backpropagate with FPNs, called Jacobian-
free backprop (JFB). Traditional backpropagation will not work effectively for
FPNs since forward propagation during training could entail hundreds or
thousands of iterations, requiring ever growing memory to store computational
graphs. On the other hand, implicit methods maintain fixed memory costs by
backpropagating “through the fixed point” and solving a Jacobian-based
equation (at potentially substantial added computational costs). The key step
to circumvent this Jacobian-based equation in our proposed method is to tune
weights by using a preconditioned gradient. Let $\ell:{\cal Y}\times{\cal
Y}\rightarrow{\mathbb{R}}$ be a smooth loss function, denoted by $\ell(x,y)$,
and consider the problem

$\min_{\Theta}{\mathbb{E}}_{d\sim{\cal D}}\big{[}\ell\left(y_{d},{\cal
M}_{\Theta}(d)\right)\big{]},$ (8)

where we abusively write ${\cal D}$ to also mean a distribution. For clarity
of presentation, in the remainder of this section we notationally suppress the
dependencies on weights $\Theta$ by letting $u^{\star}_{d}$ denote the fixed
point in (3). Unless noted otherwise, mapping arguments are implicit in this
section; in each implicit case, this will correspond to the entries in (3). We
begin with two standard assumptions that allow us to differentiate the network
${\cal M}_{\Theta}$.

######  Assumption 3.1.

The mappings $S_{\Theta}$ and $T_{\Theta}$ are continuously differentiable
with respect to $u$ and $\Theta$.

######  Assumption 3.2.

The weights $\Theta$ may be written as a tuple
$\Theta=(\theta_{S},\theta_{T})$ such that weight paramaterization of
$S_{\Theta}$ and $T_{\Theta}$ depend only on $\theta_{S}$ and $\theta_{T}$,
respectively.444This assumption is easy to ensure in practice. For notational
brevity, we use the subscript $\Theta$ throughout.

Let ${\cal J}_{\Theta}$ be defined as the identity operator, denoted by
$\mathrm{I}$, minus the Jacobian555Under Assumption 2.1, the Jacobian ${\cal
J}_{\Theta}$ exists almost everywhere. However, the presentation is much
cleaner by assuming smoothness. of $T_{\Theta}$ at $(u,d)$, i.e.

${\cal
J}_{\Theta}(u;d)\triangleq\mathrm{I}-\frac{\mathrm{d}T_{\Theta}}{\mathrm{d}u}(u;d).$
(9)

Following backpropagation approaches for existing equilibrium models
[winston2020monotone, bai2019deep], we can differentiate both sides of the
fixed point relation in (3) to obtain, by the implicit function theorem,

$\dfrac{\mathrm{d}u_{d}^{\star}}{\mathrm{d}\Theta}=\dfrac{\partial
T_{\Theta}}{\partial
u}\dfrac{\mathrm{d}u_{d}^{\star}}{\mathrm{d}\Theta}+\dfrac{\partial
T_{\Theta}}{\partial\Theta}\ \ \ \Longrightarrow\ \ \
\dfrac{\mathrm{d}u_{d}^{\star}}{\mathrm{d}\Theta}={\cal
J}_{\Theta}^{-1}\cdot\frac{\partial T_{\Theta}}{\partial\Theta},$ (10)

where ${\cal J}_{\Theta}^{-1}$ exists whenever ${\cal J}_{\Theta}$ exists (see
Lemma A.1). Using the chain rule gives the loss gradient

$\frac{\mathrm{d}}{\mathrm{d}\Theta}\left[\ell(y_{d},{\cal
M}_{\Theta}(d))\right]=\frac{\partial\ell}{\partial
y}\left[\frac{\mathrm{d}S_{\Theta}}{\mathrm{d}u}\frac{\mathrm{d}u^{\star}}{\mathrm{d}\Theta}+\frac{\partial
S_{\Theta}}{\partial\Theta}\right]=\frac{\partial\ell}{\partial
y}\left[\frac{\mathrm{d}S_{\Theta}}{\mathrm{d}u}{\cal
J}_{\Theta}^{-1}\frac{\partial T_{\Theta}}{\partial\Theta}+\frac{\partial
S_{\Theta}}{\partial\Theta}\right].$ (11)

The matrix ${\cal J}_{\Theta}$ satisfies the coercivity inequality (see Lemma
A.1)

$\left<u,{\cal
J}_{\Theta}^{-1}u\right>\geq\dfrac{1-\gamma}{(1+\gamma)^{2}}\|u\|^{2},\ \ \
\mbox{for all $u\in{\cal U}$.}$ (12)

Intuitively, this property makes it seem possible to remove ${\cal
J}_{\Theta}^{-1}$ from (11) and backpropagate using

$p_{\Theta}\triangleq-\frac{\partial\ell}{\partial
y}\left(\frac{\mathrm{d}S_{\Theta}}{\mathrm{d}u}\frac{\partial
T_{\Theta}}{\partial\Theta}+\frac{\partial
S_{\Theta}}{\partial\Theta}\right).$ (13)

The omission of ${\cal J}_{\Theta}^{-1}$ admits to two straightforward
interpretations. First note $p_{\Theta}$ is precisely the gradient of the
expression $\ell(y_{d},S_{\Theta}(T_{\Theta}(u_{d}^{\star};d))),$ treating
$u_{d}^{\star}$ as a constant independent of $\Theta$. Furthermore, the
argument $S_{\Theta}(T_{\Theta}(u_{d}^{\star};d))$ yields the same output as
${\cal M}_{\Theta}(d)$. The distinction is that using
$S_{\Theta}(T_{\Theta}(u_{d}^{\star};d))$ assumes, perhaps by chance, the user
chose the first iterate $u^{1}$ in their fixed point iteration (see Algorithm
1) to be precisely the fixed point $u_{d}^{\star}$. This makes the iteration
trivial, “converging” in one iteration. We can simulate this behavior by using
the fixed point iteration to find $u_{d}^{\star}$ and only backpropagating
through the final step of the fixed point iteration, as shown in Figure 3.

The second interpretation of our approach is to view our backpropagation as a
form of preconditioned gradient descent. Since the weights $\Theta$ typically
lay in a space of much higher dimension than the latent space ${\cal U}$, the
Jacobians $\partial S_{\Theta}/\partial\Theta$ and $\partial
T_{\Theta}/\partial\Theta$ will effectively always have full column rank. We
leverage this fact via the following assumption.

######  Assumption 3.3.

Under Assumption 3.2, given any weights $\Theta=(\theta_{S},\theta_{T})$ and
data $d$, the matrix

$M\triangleq\ \left[\begin{array}[]{cc}\frac{\partial
S_{\Theta}}{\partial\theta_{S}}&0\\\ 0&\frac{\partial
T_{\Theta}}{\partial\theta_{T}}\end{array}\right]$ (14)

has full column rank and is sufficiently well conditioned to satisfy the
inequality666The term $\gamma$ here refers to the contraction factor in (6).

$\kappa(M^{\top}M)=\dfrac{\lambda_{\mathrm{max}}(M^{\top}M)}{\lambda_{\mathrm{min}}(M^{\top}M)}\leq\dfrac{1}{\gamma}.$
(15)

######  Remark 3.1.

The conditioning portion of the above assumption is useful for bounding the
worst-case behavior in our analysis. However, we found it unnecessary to
enforce this in our experiments for effective training (e.g. see Figure 5),
which we hypothesize is justified because worst case behavior rarely occurs in
practice and we train using averages of $p_{\Theta}$ for samples drawn from a
large data set.

The full column rank of $M$ enables us to rewrite $p_{\Theta}$ as a
preconditioned gradient, i.e.

$p_{\Theta}=\underbrace{\left(M\left[\begin{array}[]{cc}\mathrm{I}&0\\\
0&{\cal J}_{\Theta}\end{array}\right]M^{+}\right)}_{\mbox{preconditioning
term}}\frac{\mathrm{d}\ell}{\mathrm{d}\Theta},$ (16)

where $M^{+}$ is the Moore-Penrose pseudo inverse [moore1920reciprocal,
penrose1955generalized]. These insights lead to our main result.

Figure 3: Diagram of backpropagation schemes for recurrent implicit depth and
fixed depth models. Forward propagation is tracked via solid arrows point to
the right. Backpropagation is shown via dashed arrows pointing to the left.
For sake of illustration, in this diagram we assume the depth of feed forward
propagation happens to agree between the traditional fixed depth feed forward
approach and the implicit depth approaches; however, this is not true in
general.

######  Theorem 3.1.

If Assumptions 2.1, 3.1, 3.2, and 3.3 hold for given weights $\Theta$ and data
$d$, then

$p_{\Theta}\triangleq-\left[\frac{\partial\ell}{\partial
y}\left(\frac{\mathrm{d}S_{\Theta}}{\mathrm{d}u}\frac{\partial
T_{\Theta}}{\partial\Theta}+\frac{\partial
S_{\Theta}}{\partial\Theta}\right)\right]_{(u,d)=(u_{d}^{\star},d)}$ (17)

forms a descent direction for $\ell(y_{d},{\cal M}_{\Theta}(d))$ with respect
to $\Theta$.

Theorem 3.1 has two practical implications. The first is we can leverage the
memory efficiency of implicit depth models where memory cost is constant with
respect to the “depth”, i.e. number of applications of $T_{\Theta}$. The
second takeaway is we can avoid difficult computations associated with ${\cal
J}_{\Theta}^{-1}$ in (11) (i.e. solving an associated linear system/adjoint
equation) in implicit depth literature [chen2018neural, dupont2019ANODEs,
bai2019deep, winston2020monotone]. Thus, our scheme more naturally applies to
general multilayered $T_{\Theta}$ and is substantially simpler to code. Our
scheme is juxtaposed in Figure 3 with classic and Jacobian-based schemes.

Two additional considerations must be made when determining the efficacy of
training a network using (17) rather than Jacobian-based gradients (11):

  * $\blacktriangleright$

Does use of (17) degrade training/testing performance relative to (11)?

  * $\blacktriangleright$

Is the gradient (17) resilient to errors in estimates of the fixed point
$u_{d}^{\star}$?

The first answer is that our training scheme simply takes a different path to
minimizers than using gradients with the implicit model. Thus, for nonconvex
problems, one should not expect the results to be the same. Our experiments
(see Section 5) using (17) for training yield state-of-the-art results.
Furthermore, (17) outperforms (11) in our experiments (when applied to
identical models). This leads us to believe, on experimental grounds, that
(17) is at least as effective as (11). The second inquiry is answered in part
by the following corollary, which states our proposed approach yields descent
even when an approximation of the fixed point $u_{d}^{\star}$ is used.

######  Corollary 3.1.

Given weights $\Theta$ and data $d$, there exists $\varepsilon>0$ such that if
$u_{d}^{\varepsilon}\in{\cal U}$ satisfies
$\|u_{d}^{\varepsilon}-u_{d}^{\star}\|\leq\varepsilon$ and the assumptions of
Theorem 3.1 hold, then

$p_{\Theta}^{\varepsilon}\triangleq-\left[\frac{\partial\ell}{\partial
y}\left(\frac{\mathrm{d}S_{\Theta}}{\mathrm{d}u}\frac{\partial
T_{\Theta}}{\partial\Theta}+\frac{\partial
S_{\Theta}}{\partial\Theta}\right)\right]_{(u,d)=(u_{d}^{\varepsilon},d)}$
(18)

is a descent direction for the loss function $\ell(y_{d},{\cal
M}_{\Theta}(u_{d}^{\star},d))$ with respect to $\Theta$.

We are not aware of any analogous results for error tolerances in the implicit
depth literature.

#### Coding Backpropagation

Implicit Forward + Proposed Backprop

    u_fxd_pt = find_fixed_point(in_feature)
    y = explicit_model(u_fxd_pt, in_feature)
    loss = criterion(y, labels)
    loss.backward()
    optimizer.step()

Figure 4: Sample PyTorch code for backpropagation888Please see our implemented
codes for complete details.

The backpropagation of our scheme is very similar to that of a standard
backpropagation. We illustrate this in the sample of PyTorch
[paszke2017automatic] code in Figure 4. Here `explicit_model` represents
$S_{\Theta}(T_{\Theta}(u;d))$. The fixed point `u_fxd_pt` is computed by
successively applying $T_{\Theta}$ (see Algorithm 1) in `eval` mode. With this
fixed point, `explicit_model` evaluates and returns
$S_{\Theta}(T_{\Theta}(u_{d}^{\star},d))$ to `y` in `train` mode (to create
the computational graph). Thus, our scheme coincides with standard
backpropagation through an explicit network with _one_ latent space layer. On
the other hand, standard implicit models backpropagate by solving a linear
system to apply ${\cal J}_{\Theta}^{-1}$ as in (11). That approach also
requires users to manually update the parameters, use more computational
resources, and make considerations (e.g. conditioning of ${\cal
J}_{\Theta}^{-1}$) for the particular network architecture used.

##  4 Related Works

#### Continuous Network Models

A well-known example of models based on a continuous view of deep networks is
neural ODEs [chen2018neural, dupont2019ANODEs, ruthotto2021introduction].
Neural ODEs leverage known connections between deep residual networks and
discretizations of differential equations [haber2017stable,
weinan2017proposal, ruthotto2019deep, chang2018reversible, finlay2020train,
lu2018beyond], and replace these discretizations by black-box ODE solvers in
forward and backward passes. The implicit property of these models arise from
their method for computing gradients. Rather than backpropagate through each
layer, backpropagation is instead done by solving the adjoint equation
[jameson1988aerodynamic] using a blackbox ODE solver as well. This allows the
user to alleviate the memory costs of backpropagation through deep neural
networks by solving the adjoint equation at additional computational costs. A
drawback is that the adjoint equation must be solved to high-accuracy;
otherwise, a descent direction is not necessarily guaranteed
[gholami2019anode, zhang2019anodev2, onken2020discretize, onken2020ot].
Another related work [gomez2017reversible] uses reversible residual networks
(ResNets) to achieve fixed memory costs; this is similar in spirit to the
present work, but it is specialized to ResNets and, as the authors noted, is
prone to accumulating numerical errors since it backpropagates through all
layers (unlike FPNs). We further note the adjoint equation for these models is
analogous to the linear system that is solved in fixed point models (discussed
next).

Figure 5: CIFAR10 performance using identical models/configurations, but with
two backpropagation schemes: our proposed method (blue) and the standard
Jacobian-based backpropagation in (11) (green), with fixed point tolerance
$\epsilon=10^{-4}$. Our method is faster and yields better test accuracy.
MNIST
---
Method | Model size | Acc.
Explicit (Fully Connected) | 84K | 98.2%
MON (Fully Connected)‡ | 84K | 98.2%
FPN (Fully Connected) | 84K | 98.6%
[3.5pt/5pt] |  |
Explicit (Convolution) | 84K | 99.4%
Neural ODE† (Convolution) | 84K | 96.4%
Aug. Neural ODE† (Convolution) | 84K | 98.2%
MON (Convolution)‡ | 84K | 99.2%
FPN (Convolution) | 84K | 99.5%
SVHN
Method | Model size | Acc.
Explicit | 154K | 92.4%
Neural ODE† | 172K | 81.0%
Aug. Neural ODE† | 172K | 83.5%
MON (Single conv)‡ | 172K | 88.8%
MON (Multi-tier lg)‡ | 170K | 92.3%
FPN (ours) | 154K | 92.8%
CIFAR-10
Method | Model size | Acc.
Explicit | 172K | 80.1%
Neural ODE† | 172K | 53.7%
Aug. Neural ODE† | 172K | 60.6%
MON (Single conv)‡ | 172K | 74.1%
FPN (ours) | 172K | 80.7%
[3.5pt/5pt] |  |
MON (Multi-tier lg)‡∗ | 1.01M | 89.7%
FPN (ours)∗ | 860K | 92.1%
Table 1: Test accuracy of FPNs (our method) compared to Neural ODEs, Augmented
Neural ODEs, and MONs; †as reported in [dupont2019ANODEs]; ‡as reported in
[winston2020monotone]; *with data augmentation

#### Fixed Point Models

Another class of implicit models involve formulating machine learning tasks as
fixed point problems. Deep Equilibrium (DEQ) Networks [bai2019deep] find a
fixed point of a nonlinear dynamical system corresponding to an effectively
infinite-depth “weight-tied, input-injected” network. In this setting, forward
propagation consists of applying a root-finding method to find a fixed point,
and backpropagation involves solving a Jacobian-based equation arising from
the implicit function theorem [krantz2012implicit]. Multiscale DEQs (MDEQs)
[bai2020multiscale] expands upon DEQs by learning multiple feature resolutions
simultaneously. As is well known in fixed point theory, naïve forward
propagation of a fixed point scheme is not necessarily stable. To this end,
Monotone Operator Equilibrium Networks (MON) [winston2020monotone] leverage
monotone operator theory to create stable forward and backward propagation by
using special network architectures. These networks must be 1-layer networks
with activations based on proximal operators (e.g. ReLU). Similarly to MON,
FPNs leverage fixed point theory. However, our approach generalizes MONs by
allowing a wider class of network architectures. And, while our approach falls
under the category of implicit depth models, our backpropagation scheme
_circumvents_ the need to solve the Jacobian-based equation presented in (9).
In [ghaoui2019implicit] the need to solve a Jacobian-based equation is also
avoided by eschewing backprop entirely and using the so-called “Fenchel
divergence formulation” of the training loss. However, that approach appears
to be computationally taxing, and the authors do not implement it on any
large-scale benchmark data sets.

##  5 Experiments

We demonstrate the effectiveness our backpropagation scheme on FPNs.
Experiments are done using PyTorch [paszke2017automatic]. All networks were
designed such that Assumptions 2.1 and 3.2 hold.999A weaker version of
Assumption 3.1 also holds in practice, i.e. differentiability almost
everywhere. In particular, we bound the largest singular values of latent
space network mappings and choose standard 1-Lipschitz activations functions
(e.g. see [gao2017properties, combettes2020lipschitz, combettes2020deep]). All
experiments are run on a single NVIDIA TITAN X GPU with 12GB RAM. More details
about experimental settings can be found in Appendix B.

###  5.1 Classification

We train FPNs on three benchmark image classification datasets: SVHN
[netzer2011reading], MNIST [lecun2010mnist], and CIFAR-10
[krizhevsky2009learning]. Table 1 compares our results with state-of-the-art
results for implicit models, including Neural ODEs [chen2018neural], Augmented
Neural ODEs [dupont2019ANODEs], and MONs [winston2020monotone]. We also
compare with an explicit version of our FPNs given in (1). The explicit models
are trained with the exact same setup as those used for FPN (with the
exception of using one latent space layer). Table 1 shows FPNs, trained with
our scheme, substantially outperform all the ODE-based models as well as MONs
using similar or fewer parameters. Appendix C contains further experimental
details.

###  5.2 Comparing Backpropagation

Table 2 compares the performance between using the standard Jacobian-based
backpropagation and our proposed backpropagation. The experiments are
performed on all the datasets described in Section 5.1. To apply the Jacobian-
based in (10), we use the conjugate gradient (CG) method on an associated set
of normal equations in similar manner to [liao2018reviving]. We set a maximum
of 500 CG iterations. The remainder of the experimental settings are kept the
same as those from our proposed approach (and are therefore not tuned to the
best of our ability). The main purpose of the Jacobian-based results in Figure
5 and Table 2, are the significant speedups in training time (in each case
obtaining roughly 2X or more speedup) while maintaining a competitive accuracy
that outperforms previous state-of-the-art implicit models. Further plots and
comparisons are provided in Appendix B.

| Dataset | Avg time per epoch (s) | # of $\mathcal{J}$ mat-vec products | Accuracy %
---|---|---|---|---
With Jacobian (Standard) | MNIST | 218.5 | $5.9\times 10^{8}$ | 99.0
SVHN |  186.3 |  $2.3\times 10^{8}$ |  90.6
CIFAR10 | 293.8 | $4.4\times 10^{8}$ | 90.5
Jacobian-free (Proposed) | MNIST | 50.2 | 0 | 99.5
SVHN |  56.7 |  0 |  92.7
CIFAR10 | 105.4 | 0 | 92.1
Table 2:  Backpropagation comparison for FPN using the Jacobian-based equation
(first three rows) and our proposed Jacobian-free backpropagation. “Mat-vecs”
denotes matrix-vector products.

##  6 Conclusion

We present Fixed Point Networks (FPNs), a class of models that generalizes
existing implicit models to a wider choice of architectures (while maintaining
convergence guarantees). A new and simple Jacobian-free backpropagation (JFB)
scheme is proposed that avoids solving a Jacobian-based equation. This enables
our approach to enjoy fixed memory costs (regardless of depth), be easy to
code, and yield faster backpropagation than previously possible. We show our
scheme is theoretically justified (even with errors present). Our experiments
also show our scheme is efficient and effective. Extensions to FPNs will
enable satisfaction of additional constraints for imaging and phase retrieval
[klibanov1986determination, fienup1982phase, heaton2020projecting,
fung2020multigrid, kan2020pnkh], geophysics [haber2014computational,
fung2019multiscale, fung2019uncertainty], and game theory [von1959theory,
myerson2013game, lin2020apac, li2019cubic, ruthotto2020machine].

## Broader Impact

The rapid rise of deep learning has been largely due to the growth in
computational capabilities. Despite this, a key shortcoming associated with
creating deep learning models is that making them deeper and, thus, more
expressible requires ever increasing memory. Thus, for practical devices,
computer memory limits model depth. Existing implicit depth networks present a
way of circumventing this memory limit; however, this comes with the trade off
of much more computational costs to train networks. This work presents a third
option, implicit depth without increasing the computational costs of training.
This approach increases the efficiency of computational resources both with
respect to memory and computation. Moreover, by using implicit depth models,
the same model can be deployed with settings attuned to both high power (e.g.
servers) and low power (e.g. mobile) devices.

## Acknowledgments and Disclosure of Funding

Samy Wu Fung and Stanley Osher are supported by AFOSR MURI FA9550-18-1-0502,
AFOSR Grant No. FA9550-18-1-0167, and ONR Grants N00014-18-1-2527 snf
N00014-17-1-21. Daniel McKenzie and Qiuwei Li are supported by AFOSR MURI
FA9550-18-1-0502. Howard Heaton is supported by the National Science
Foundation (NSF) Graduate Research Fellowship under Grant No. DGE-1650604. Any
opinion, findings, and conclusions or recommendations expressed in this
material are those of the authors and do not necessarily reflect the views of
the NSF.

##  Appendix A Proofs

This section provides proofs for the results of Section 3. For the reader’s
convenience, we restate all results before proving them.

######  Lemma A.1.

If Assumption 2.1 and 3.1 hold, then ${\cal J}_{\Theta}$ in (9) exists and

$\left<u,{\cal J}_{\Theta}u\right>\geq(1-\gamma)\|u\|^{2},\ \ \ \mbox{for all
$u\in{\cal U}$.}$ (19)

Additionally, ${\cal J}_{\Theta}$ is invertible, and its inverse ${\cal
J}_{\Theta}^{-1}$ satisfies the coercivity inequality

$\left<u,{\cal
J}_{\Theta}^{-1}u\right>\geq\dfrac{1-\gamma}{(1+\gamma)^{2}}\|u\|^{2},\ \ \
\mbox{for all $u\in{\cal U}$.}$ (20)

###### Proof.

We proceed in the following manner. First we establish the coercivity
inequality (19) (Step 1). This is used to show ${\cal J}_{\Theta}$ is
invertible (Step 2). The previous two results are then combined to establish
the inequality (20) (Step 3). All unproven results that are quoted below about
operators are standard and may be found in any standard functional analysis
text (e.g. [kreyszig1978introductory]).

Step 1. To obtain our coercivity inequality, we identify a bound on the
operator norm for $\partial T_{\Theta}/\partial u$. Fix any unit vector
$v\in{\cal U}$. Then, by the definition of differentiation,

$\frac{\mathrm{d}T_{\Theta}}{\mathrm{d}u}v=\lim_{\varepsilon\rightarrow
0^{+}}\dfrac{T_{\Theta}(u^{\star}+\varepsilon
v;d)-T_{\Theta}(u^{\star};d)}{\|(u^{\star}+\varepsilon
v)-u^{\star}\|}=\lim_{\varepsilon\rightarrow
0^{+}}\dfrac{T_{\Theta}(u^{\star}+\varepsilon
v;d)-T_{\Theta}(u^{\star};d)}{\varepsilon}.$ (21)

Thus,

$\left\|\frac{\mathrm{d}T_{\Theta}}{\mathrm{d}u}v\right\|=\left\|\lim_{\varepsilon\rightarrow
0^{+}}\dfrac{T_{\Theta}(u^{\star}+\varepsilon
v;d)-T_{\Theta}(u^{\star};d)}{\varepsilon}\right\|=\lim_{\varepsilon\rightarrow
0^{+}}\dfrac{\|T_{\Theta}(u^{\star}+\varepsilon
v;d)-T_{\Theta}(u^{\star};d)\|}{\varepsilon},$ (22)

where the first equality follows from (21) and the second holds by the
continuity of norms. Combining (23) with the Lipschitz assumption (6) gives
the upper bound

$\left\|\frac{\mathrm{d}T_{\Theta}}{\mathrm{d}u}v\right\|\leq\lim_{\varepsilon\rightarrow
0^{+}}\dfrac{\gamma\|(u^{\star}+\varepsilon
v)-u^{\star}\|}{\varepsilon}=\gamma.$ (23)

Because the upper bound relation in (23) holds for an arbitrary unit vector
$v\in{\cal U}$, we deduce

$\left\|\frac{\mathrm{d}T_{\Theta}}{\mathrm{d}u}\right\|\triangleq\sup\left\\{\left\|\frac{\mathrm{d}T_{\Theta}}{\mathrm{d}u}v\right\|:\
\|v\|=1\right\\}\leq\gamma.$ (24)

That is, the operator norm is bounded by $\gamma$. Together the Cauchy-Schwarz
inequality and (24) imply

$\left<u,\frac{\mathrm{d}T_{\Theta}}{\mathrm{d}u}u\right>\leq\left\|\frac{\mathrm{d}T_{\Theta}}{\mathrm{d}u}\right\|\|u\|^{2}\leq\gamma\|u\|^{2},\
\ \ \mbox{for all}\ u\in{\cal U}.$ (25)

Whence the bilinear form $\left<\ \cdot\ ,\ {\cal J}_{\Theta}\ \cdot\ \right>$
is $(1-\gamma)$ coercive, i.e.

$\left<u,{\cal
J}_{\Theta}u\right>=\|u\|^{2}-\left<u,\frac{\mathrm{d}T_{\Theta}}{\mathrm{d}u}u\right>\geq(1-\gamma)\|u\|^{2},\
\ \ \mbox{for all}\ u\in{\cal U}.$ (26)

Step 2. Consider any kernel element $w\in\ker({\cal J}_{\Theta})$. Then (26)
implies

$(1-\gamma)\|w\|^{2}\leq\left<w,{\cal J}_{\Theta}w\right>=\left<w,0\right>=0\
\ \ \Longrightarrow\ \ \ (1-\gamma)\|w\|^{2}\leq 0\ \ \ \Longrightarrow\ \ \
w=0.$ (27)

Consequently, the kernel of ${\cal J}_{\Theta}$ is trivial, i.e.

$\ker({\cal J}_{\Theta})\triangleq\\{u:{\cal J}_{\Theta}u=0\\}=\\{0\\}.$ (28)

Thus, the linear operator ${\cal J}_{\Theta}$ is invertible.

Step 3. By (23) and an elementary result in functional analysis,

$\|{\cal J}_{\Theta}^{\top}{\cal J}_{\Theta}\|=\|{\cal
J}_{\Theta}\|^{2}\leq\left(\|\mathrm{I}\|+\left\|\frac{\mathrm{d}T_{\Theta}}{\mathrm{d}u}\right\|\right)^{2}\leq(1+\gamma)^{2}.$
(29)

Hence

$\|u\|^{2}=\left<u,u\right>=\left<{\cal J}_{\Theta}^{-1}u,({\cal
J}_{\Theta}^{\top}{\cal J}_{\Theta}){\cal
J}_{\Theta}^{-1}u\right>\leq(1+\gamma)^{2}\left\|{\cal
J}_{\Theta}^{-1}u\right\|^{2},\ \ \ \mbox{for all $u\in{\cal U}$.}$ (30)

Combining (26) and (30) reveals

$\dfrac{1-\gamma}{(1+\gamma)^{2}}\left<u,u\right>\leq(1-\gamma)\|{\cal
J}_{\Theta}^{-1}u\|^{2}\leq\left<{\cal J}^{-1}_{\Theta}u,{\cal
J}_{\Theta}({\cal J}_{\Theta}^{-1}u)\right>=\left<{\cal
J}_{\Theta}^{-1}u,u\right>,\ \ \ \mbox{for all $u\in{\cal U}$.}$ (31)

This establishes (20), and we are done. ∎

######  Lemma A.2.

If $A\in{\mathbb{R}}^{t\times t}$ is symmetric with positive eigenvalues,

$\overline{\lambda}\triangleq\dfrac{\lambda_{\mathrm{max}}(A)+\lambda_{\mathrm{min}}(A)}{2}\
\ \ \mbox{and}\ \ \ S\triangleq\overline{\lambda}\mathrm{I}-A,$ (32)

then

$\|S\|=\dfrac{\lambda_{\mathrm{max}}(A)-\lambda_{\mathrm{min}}(A)}{2}.$ (33)

###### Proof.

Since $A$ is symmetric, the spectral theorem asserts it possesses a set of
eigenvectors that form an orthogonal basis for ${\mathbb{R}}^{t}$. This same
basis forms the set of eigenvectors for $\overline{\lambda}\mathrm{I}-A$, with
eigenvalues denoted by $\\{\lambda_{i}\\}_{i=1}^{t}$. So, there exists
orthogonal $P\in{\mathbb{R}}^{t\times t}$ and diagonal $\Lambda$ with entries
given by each of the eigenvalues $\lambda_{i}$ such that

$S=\overline{\lambda}\mathrm{I}-P^{\top}\Lambda
P=P^{\top}\left(\overline{\lambda}\mathrm{I}-\Lambda\right)P.$ (34)

Substituting this equivalence into the definition of the operator norm yields

$\|S\|\triangleq\sup\left\\{\left\|S\xi\right\|:\
\|\xi\|=1\right\\}=\sup\left\\{\|P^{\top}(\overline{\lambda}I-\Lambda)P\xi\|:\
\|\xi\|=1\right\\}.$ (35)

Leveraging the fact $P$ is orthogonal enables the supremum above to be
restated via

$\|S\|=\sup\left\\{\|(\overline{\lambda}I-\Lambda)P\xi\|:\
\|\xi\|=1\right\\}=\sup\left\\{\|(\overline{\lambda}I-\Lambda)\zeta\|:\
\|\zeta\|=1\right\\}.$ (36)

Because $\overline{\lambda}\mathrm{I}-\Lambda$ is diagonal, (36) implies

$\|S\|=\max_{i\in[t]}|\overline{\lambda}-\lambda_{i}|=\dfrac{\lambda_{\mathrm{max}}(A)-\lambda_{\mathrm{min}}(A)}{2},$
(37)

and the proof is complete. ∎

Theorem 3.1. If Assumptions 2.1, 3.1, 3.2, and 3.3 hold for given weights
$\Theta$ and data $d$, then

$p_{\Theta}\triangleq-\left[\frac{\partial\ell}{\partial
y}\left(\frac{\mathrm{d}S_{\Theta}}{\mathrm{d}u}\frac{\partial
T_{\Theta}}{\partial\Theta}+\frac{\partial
S_{\Theta}}{\partial\Theta}\right)\right]_{(u,d)=(u_{d}^{\star},d)}$ (38)

forms a descent direction for $\ell(y_{d},{\cal M}_{\Theta}(d))$ with respect
to $\Theta$.

###### Proof.

To complete the proof, it suffices to show

$\left<\frac{\mathrm{d}\ell}{\mathrm{d}\Theta},p_{\Theta}\right><0,\ \ \
\mbox{for all}\ \ \frac{\mathrm{d}\ell}{\mathrm{d}\Theta}\neq 0.$ (39)

Let any weights $\Theta$ and data $d$ be given, and assume the gradient
$\mathrm{d}\ell/\mathrm{d}\Theta$ is nonzero. We proceed in the following
manner. First we show $p_{\Theta}$ is equivalent to a preconditioned gradient
(Step 1). We then
show $M^{\top}\mathrm{d}\ell/\mathrm{d}\Theta$ is nonzero, with $M$ as in (14)
of Assumption 3.3 (Step 2). These two results are then combined to verify the
descent inequality (39) for the provided $\Theta$ and $d$ (Step 3).

Step 1. Denote the dimension of each component of the gradient
$\mathrm{d}\ell/\mathrm{d}\Theta$ using101010We assumed each space is a real-
valued finite dimensional Hilbert space, making it equivalent to some
Euclidean space. So, it suffices to show everything in Euclidean spaces.

$\frac{\partial T_{\Theta}}{\partial\Theta}\in{\mathbb{R}}^{p\times n},\ \ \ \
{\cal J}_{\Theta}^{-1}\in{\mathbb{R}}^{n\times n},\ \ \ \ \frac{\partial
S_{\Theta}}{\partial\Theta}\in{\mathbb{R}}^{p\times c},\ \ \
\frac{\mathrm{d}S_{\Theta}}{\mathrm{d}u}\in{\mathbb{R}}^{n\times c},\ \ \ \
\frac{\partial\ell}{\partial y}\in{\mathbb{R}}^{c\times 1}.$ (40)

Combining each of these terms yields the gradient expression111111In the main
text, the ordering was used to make clear application of the chain rule, but
here we reorder terms to get consistent dimensions in each matrix operation.

$\frac{\mathrm{d}\ell}{\mathrm{d}\Theta}=\left[\frac{\partial
T_{\Theta}}{\partial\Theta}{\cal
J}_{\Theta}^{-1}\frac{\mathrm{d}S_{\Theta}}{\mathrm{d}u}+\frac{\mathrm{d}S_{\Theta}}{\mathrm{d}\Theta}\right]\frac{\partial\ell}{\partial
y}.$ (41)

By Assumption 3.2, $S_{\Theta}$ and $T_{\Theta}$ depend on separate components
of $\Theta=(\theta_{S},\theta_{T})$. Thus,

$\frac{\mathrm{d}\ell}{\mathrm{d}\Theta}=\left[\begin{array}[]{c}\frac{\partial
S_{\Theta}}{\partial\theta_{S}}\\\\[3.0pt] \frac{\partial
T_{\Theta}}{\partial\theta_{T}}{\cal
J}_{\Theta}^{-1}\frac{\mathrm{d}S_{\Theta}}{\mathrm{d}u}\end{array}\right]\frac{\partial\ell}{\partial
y}=\underbrace{\left[\begin{array}[]{cc}\frac{\partial
S_{\Theta}}{\partial\theta_{S}}&0\\\ 0&\frac{\partial
T_{\Theta}}{\partial\theta_{T}}\end{array}\right]}_{M}\underbrace{\left[\begin{array}[]{cc}\mathrm{I}&0\\\
0&{\cal J}_{\Theta}^{-1}\end{array}\right]}_{\tilde{{\cal
J}}_{\Theta}^{-1}}\underbrace{\left[\begin{array}[]{c}\mathrm{I}\\\
\frac{\mathrm{d}S_{\Theta}}{\mathrm{d}u}\end{array}\right]\frac{\partial\ell}{\partial
y}}_{v},$ (42)

where we define121212Note this choice of $M$ coincides with the matrix $M$ in
Assumption 3.3. $M\in{\mathbb{R}}^{p\times(n+c)}$,
$\tilde{J}_{\Theta}^{-1}\in{\mathbb{R}}^{(n+c)\times(n+c)}$, and
$v\in{\mathbb{R}}^{(n+c)\times 1}$ to be the underbraced quantities. This
enables the gradient to be concisely expressed via the relation

$\frac{\mathrm{d}\ell}{\mathrm{d}\Theta}=M\tilde{{\cal J}}_{\Theta}^{-1}v,$
(43)

and our proposed gradient alternative in (38) is given by

$p_{\Theta}=-Mv.$ (44)

Because $M$ has full column rank (by Assumption 3.3), $M^{+}M=\mathrm{I}$,
enabling us to rewrite $p_{\Theta}$ via

$p_{\Theta}=-M\tilde{{\cal J}}_{\Theta}M^{+}M{\cal
J}_{\Theta}^{-1}v=-{(M\tilde{{\cal
J}}_{\Theta}M^{+})}\frac{\mathrm{d}\ell}{\mathrm{d}\Theta}.$ (45)

Hence $p_{\Theta}$ is a preconditioned gradient (n.b. the preconditioner is
not necessarily symmetric).

Step 2. Set

$w\triangleq
M^{\top}\frac{\mathrm{d}\ell}{\mathrm{d}\Theta}=M^{\top}M\tilde{{\cal
J}}_{\Theta}^{-1}v.$ (46)

The fact that $M$ has full column rank implies it has a trivial kernel. In
particular,

$0\neq\frac{\mathrm{d}\ell}{\mathrm{d}\Theta}=M\tilde{{\cal
J}}_{\Theta}^{-1}v\ \ \ \Longrightarrow\ \ \ 0\neq\tilde{{\cal
J}}_{\Theta}^{-1}v.$ (47)

Again leveraging the full column rank of $M$, we know $M^{\top}M$ is
invertible and, thus, has trivial kernel as well. This fact together with (47)
reveals

$0\neq(M^{\top}M)\tilde{{\cal J}}_{\Theta}^{-1}v=w.$ (48)

Step 3. Inserting the definition of $w$ and $p_{\Theta}$ formulation of (45)
into the scalar product in (39) yields

$\left<\frac{\mathrm{d}\ell}{\mathrm{d}\Theta},p_{\Theta}\right>=-\left<M^{\top}M\tilde{{\cal
J}}_{\Theta}^{-1}v,\tilde{{\cal J}}_{\Theta}M^{+}M\tilde{{\cal
J}}_{\Theta}^{-1}v\right>\\\ =-\left<w,\tilde{{\cal
J}}_{\theta}(M^{\top}M)^{-1}w\right>,$ (49)

noting $M^{+}=(M^{\top}M)^{-1}M^{\top}$. Let $\lambda_{+}$ and $\lambda_{-}$
be the maximum and minimum eigenvalues of $(M^{\top}M)^{-1}$, respectively.
Note $(M^{\top}M)$ is positive definite since the full column rank of $M$
implies

$\left<\xi,M^{\top}M\xi\right>=\|M\xi\|^{2}>0,\ \ \ \mbox{for all nonzero
$\xi\in{\mathbb{R}}^{n+c}$.}$ (50)

Thus, $(M^{\top}M)^{-1}$ is positive definite, making
$\lambda_{+},\lambda_{-}>0$. Let $\overline{\lambda}$ be the average of these
terms, i.e.

$\overline{\lambda}\triangleq\dfrac{\lambda_{+}+\lambda_{-}}{2}.$ (51)

Substituting in this choice of $\overline{\lambda}$ to (49) by adding and
subtracting $\overline{\lambda}\mathrm{I}$ gives the inequality

$-\left<w,\tilde{{\cal
J}}_{\theta}(M^{\top}M)^{-1}w\right>\leq-\overline{\lambda}(1-\gamma)\|w\|^{2}+\left<w,\tilde{{\cal
J}}_{\Theta}(\overline{\lambda}\mathrm{I}-(M^{\top}M)^{-1})w\right>,$ (52)

noting $\tilde{{\cal J}}_{\Theta}$ is $1-\gamma$ coercive because it is the
block diagonal composition of ${\cal J}_{\Theta}$, which is $1-\gamma$
coercive by (19 ) in Lemma A.1, and the identity matrix, which is 1-coercive.
Application of the Cauchy Schwarz inequality to the right hand side of (52)
reveals

$-\left<w,\tilde{{\cal
J}}_{\theta}(M^{\top}M)^{-1}w\right>\leq-\overline{\lambda}(1-\gamma)\|w\|^{2}+\|\tilde{{\cal
J}}_{\Theta}\|\|\overline{\lambda}\mathrm{I}-(M^{\top}M)^{-1})\|\|w\|^{2}.$
(53)

By Lemma A.2,

$\|\overline{\lambda}\mathrm{I}-(M^{\top}M)^{-1}\|=\dfrac{\lambda_{+}-\lambda_{-}}{2}.$
(54)

Similar block diagonal argument as used above to verify $\tilde{{\cal
J}}_{\Theta}$ is coercive can also be applied to bound the operator norm of
$\tilde{{\cal J}}_{\Theta}$. Indeed, (24) implies

$\|{\cal J}_{\Theta}\|\leq 1+\gamma\ \ \ \Longrightarrow\ \ \ \|\tilde{{\cal
J}}_{\Theta}\|\leq 1+\gamma.$ (55)

Hence (49), (53), (54), and (55) together yield

$\left<\frac{\mathrm{d}\ell}{\mathrm{d}\Theta},p_{\Theta}\right>\leq-\dfrac{1}{2}\big{(}(1-\gamma)(\lambda_{+}+\lambda_{-})-(1+\gamma)(\lambda_{+}-\lambda_{-})\big{)}\|w\|^{2}=-2(\lambda_{-}-\gamma\lambda_{+})\|w\|^{2}.$
(56)

The right hand expression in (56) is negative since (48) shows $w\neq 0$ and
the conditioning inequality (15) in Assumption 3.3 implies
$(\lambda_{-}-\gamma\lambda_{+})$ is positive. This verifies (39), completing
the proof. ∎

Corollary 3.1. Given weights $\Theta$ and data $d$, there exists
$\varepsilon>0$ such that if $u^{\varepsilon}\in{\cal U}$ satisfies
$\|u_{d}^{\varepsilon}-u_{d}^{\star}\|\leq\varepsilon$ and the assumptions of
Theorem 3.1 hold, then

$p_{\Theta}^{\varepsilon}\triangleq-\left[\frac{\partial\ell}{\partial
y}\left(\frac{\mathrm{d}S_{\Theta}}{\mathrm{d}u}\frac{\partial
T_{\Theta}}{\partial\Theta}+\frac{\partial
S_{\Theta}}{\partial\Theta}\right)\right]_{(u,d)=(u_{d}^{\varepsilon},d)}$
(57)

is a descent direction for the loss function $\ell(y_{d},{\cal
M}_{\Theta}(u_{d}^{\star},d))$ with respect to $\Theta$.

###### Proof.

For notational convenience, for all $\tilde{u}\in{\cal U}$, define

$p_{\Theta}(\tilde{u})\triangleq-\left[\frac{\partial\ell}{\partial
y}\left(\frac{\mathrm{d}S_{\Theta}}{\mathrm{d}u}\frac{\partial
T_{\Theta}}{\partial\Theta}+\frac{\partial
S_{\Theta}}{\partial\Theta}\right)\right]_{(u,d)=(\tilde{u},d)},$ (58)

noting $p_{\Theta}^{\varepsilon}=p_{\Theta}(u_{d}^{\varepsilon})$. Also define
the quantity

$\nabla\triangleq\frac{\mathrm{d}}{\mathrm{d}\Theta}\left[\ell(y_{d},{\cal
M}_{\Theta}(d))\right].$ (59)

Assuming $\nabla\neq 0$, it suffices to show

$\left<p_{\Theta}^{\varepsilon},\nabla\right><0.$ (60)

By the smoothness of $\ell$, $S_{\Theta}$, and $T_{\Theta}$ (see Assumption
3.1), there exists $\delta>0$ such that

$\|u-u_{d}^{\star}\|\leq\delta\ \ \ \Longrightarrow\ \ \
\left\|p_{\Theta}(u)-p_{\Theta}(u_{d}^{\star})\right\|\leq\dfrac{(\lambda_{-}-\gamma\lambda_{+})\|M^{\top}\nabla\|^{2}}{\|\nabla\|},$
(61)

where $\lambda_{+}$ and $\lambda_{-}$ are the maximum and minimum eigenvalues
of $(M^{\top}M)^{-1}$, respectively. Substituting the inequality (56) in the
proof of Theorem 3.1 into (60) reveals

$\displaystyle\left<p_{\Theta}(u),\nabla\right>$
$\displaystyle=\left<p_{\Theta}(u_{d}^{\star}),\nabla\right>+\left<p_{\Theta}(u)-p_{\Theta}(u_{d}^{\star}),\nabla\right>$
(62a)
$\displaystyle\leq-2(\lambda_{-}-\gamma\lambda_{+})\|M^{\top}\nabla\|^{2}+\left<p_{\Theta}(u)-p_{\Theta}(u_{d}^{\star}),\nabla\right>.$
(62b)

But, the Cauchy Schwarz inequality and (61) enable us to obtain the upper
bound

$|\left<p_{\Theta}(u)-p_{\Theta}(u_{d}^{\star}),\nabla\right>|\leq(\lambda_{-}-\gamma\lambda_{+})\|M^{\top}\nabla\|^{2},\
\ \ \mbox{for all $u\in B(u_{d}^{\star},\delta)$,}$ (63)

where $B(u_{d}^{\star},\delta)$ is the ball of radius $\delta$ centered about
$u_{d}^{\star}$. Combining (62) and (63) yields

$\left<p_{\Theta}(u),\nabla\right>\leq-(\lambda_{-}-\gamma\lambda_{+})\|M^{\top}\nabla\|^{2},\
\ \ \mbox{for all $u\in B(u_{d}^{\star},\delta)$.}$ (64)

In particular, this shows (60) holds when we set $\varepsilon=\delta$. ∎

##  Appendix B Classification Accuracy Plots

Figure 6: MNIST performance using identical models/configurations, but with
two backpropagation schemes: our proposed method (blue) and the standard
Jacobian-based backpropagation in (11) (green), with fixed point tolerance
$\epsilon=10^{-4}$. Our method is 4X faster and yields better test accuracy.

Figure 7: SVHN performance using identical models/configurations, but with two
backpropagation schemes: our proposed method (blue) and the standard Jacobian-
based backpropagation in (11) (green), with fixed point tolerance
$\epsilon=10^{-4}$. Our method is 3X faster and yields better test accuracy.

##  Appendix C Experimental Settings

We present the experimental settings and describe the architecture used for
each dataset. For all models, we use the Adam [kingma2015adam] optimizer for
training.

###  C.1 MNIST

We use two convolutions for the data-space portion of the network $Q(d)$ and
use one fully-connected layer for the latent space with relu activation. For
the fixed point stopping criterion, we stop whenever the fixed point satisfies
$\epsilon<10^{-6}$ or 100 iterations have occurred. We use an initial learning
rate of $10^{-5}$ and decay the learning rate by multiplying it by $0.98$
every 10 epochs. To map from the latent space to inference space,
$S_{\Theta}$, we use a fully connected layer with a relu activation function.
For the loss function, we found mean squared error to be most effective.

###  C.2 SVHN

We use two convolutions for the data-space portion of the network $Q(d)$ and
use one fully-connected layer for the latent space with leaky relu activation.
For the fixed point stopping criterion, we stop whenever the fixed point
satisfies $\epsilon<10^{-4}$ or 100 iterations have occurred. We use an
initial learning rate of $10^{-3}$ and decay the learning rate using cosine
annealing [loshchilov2016sgdr]. To map from the latent space to inference
space, $S_{\Theta}$, we use a fully connected layer with a softmax activation
function. We use cross entropy for the loss function.

###  C.3 CIFAR10

For the unaugmented CIFAR10 results, we use three convolutions for the data-
space portion of the network $Q(d)$ and use two fully-connected layer for the
latent space with leaky relu activation. For the augmented CIFAR10 results, we
use ResNet56 [he2016deep, Idelbayev18a] for the data-space portion of the
network $Q(d)$ and modify its final layer to output to a 45-dimensional latent
space (rather than originally to 10 dimensional space). For the fixed point
stopping criterion, we stop whenever the fixed point satisfies
$\epsilon<10^{-3}$ or 2000 iterations have occurred. We use an initial
learning rate of $5\times 10^{-4}$ and decay the learning rate using cosine
annealing [loshchilov2016sgdr]. To map from the latent space to inference
space, $S_{\Theta}$, we use a fully connected layer with a leaky relu
activation function. We use cross entropy for the loss function.