# A gauge-invariant differential equation for covariant Lyapunov vectors

Doriano Brogioli  Universität Bremen, Energiespeicher– und
Energiewandlersysteme, Bibliothekstraße 1, 28359 Bremen, Germany  Massimo
Marino  Liceo Scientifico Statale “A. Volta”, via Benedetto Marcello 7, 20124
Milan, Italy

###### Abstract

We introduce the concept of “covariant Lyapunov field”, which assigns all the
components of covariant Lyapunov vectors at almost all points of the phase
space of a dynamical system. We show that in ergodic systems such fields can
be characterized as the global solutions of a partial differential equation on
the phase space. Due to the arbitrariness in the choice of a multiplicative
scalar factor for the Lyapunov vector at each point of the phase space, this
partial differential equation presents a gauge invariance that is formally
analogous to that of electrodynamics.

Lyapunov exponents (LE) quantify the rate of divergence of trajectories in a
dynamical system [1]. Given a vector field $\boldsymbol{F}(\boldsymbol{x})$ on
an $n$-dimensional Riemannian manifold $X$, let us consider a trajectory
$\boldsymbol{x}(t)$, satisfying the differential equation

$\frac{\mathrm{d}}{\mathrm{d}t}x^{\mu}(t)=F^{\mu}\left[\boldsymbol{x}(t)\right]$
(1)

for $\mu=1,\dots,n$. Then for a slightly displaced trajectory
$\boldsymbol{x}(t)+\boldsymbol{\delta x}(t)$ we get at first order in
$\boldsymbol{\delta x}(t)$ the equation

$\frac{\mathrm{d}}{\mathrm{d}t}\delta
x^{\mu}(t)=\partial_{\nu}F^{\mu}\left[\boldsymbol{x}(t)\right]\delta
x^{\nu}(t)\,,$ (2)

where $\partial_{j}$ is the partial derivative with respect to the $j$-th
coordinate and the terms are summed over the repeated Greek indices (Einstein
notation).

It is typically observed that the norm of $\boldsymbol{\delta x}(t)$,
asymptotically for $t\to+\infty$, behaves exponentially:

$\lVert\boldsymbol{\delta x}(t)\rVert=e^{\lambda^{+}t+o(t)}\,,$ (3)

where $o(t)/t$ vanishes in the limit $t\to+\infty$. The real parameter

$\lambda^{+}=\lim_{t\to+\infty}\frac{\ln\lVert\boldsymbol{\delta
x}(t)\rVert}{t}\,,$

which is called forward Lyapunov exponent (LE), depends both on the initial
position $\boldsymbol{x}(0)=\boldsymbol{x}_{0}$ and on the initial
displacement vector $\boldsymbol{\delta x}(0)=\boldsymbol{\delta x}_{0}$.
Under quite general hypotheses, it can be proved instead that the value of the
LE does not depend on the particular choice of the metric tensor used on the
manifold $X$ for the calculation of the norm appearing in Eq. (3).

If $\mu$ is a measure on $X$ which is preserved by the flux generated by
$\boldsymbol{F}(\boldsymbol{x})$, Oseledec’s theorem [2] says that Eq. (3)
holds for any $\boldsymbol{\delta x}_{0}$ at almost all points
$\boldsymbol{x}_{0}\in X$ with respect to the measure $\mu$. It is easy to see
that, given $\boldsymbol{x}_{0}$, there can be at most $n$ distinct forward
LEs $\lambda_{j}^{+}$ as a function of $\boldsymbol{\delta x}_{0}$, thanks to
the linearity of Eq. (2) in $\boldsymbol{\delta x}(t)$. For a rigorous
discussion of the conditions for the existence of the LEs we refer the readers
to Ref. [3].

LEs give a deep characterization of dynamical systems and have proved to be an
invaluable tool for the analysis of chaotic systems and attractors, either in
numerical calculations or experimental data [4]. In particular, they provide
the Kolmogorov-Sinai entropy [5] and the Kaplan-Yorke dimension of an
attractor [6], which is an upper bound for the information dimension of the
system.

The whole set of forward LEs $\lambda_{j}^{+}$, for a given
$\boldsymbol{x}_{0}$, can be calculated by suitable numerical procedures [7,
8] which, as a by-product, also return a set of $n$ so-called forward
“orthonormal Lyapunov vectors” (OLVs): each forward LE $\lambda_{j}^{+}$ is
obtained by taking one of the forward OLVs as initial displacement vector
$\boldsymbol{\delta x}_{0}$. As the name suggests, the OLVs form an
orthonormal set with respect to the chosen metric tensor on $X$. However, at
variance with the LEs, these OLVs do depend on the arbitrary choice of such a
metric tensor, and for this reason they do not represent an intrinsic
characterization of the dynamical system.

A different point of view arises when the forward dynamics is compared with
the backward dynamics [1]. In analogy with Eq. (3), an exponential behaviour
of the displacement vector is also observed looking at the evolution back in
time, for $t\to-\infty$:

$\lVert\boldsymbol{\delta x}(t)\rVert=e^{\lambda^{-}|t|+o(t)}\,,$ (4)

where $\lambda^{-}$ represents the backward LE and $o(t)/t$ vanishes for
$t\to-\infty$. As for the forward LEs $\lambda_{j}^{+}$, there can be at most
$n$ distinct backward LEs $\lambda_{j}^{-}$ as a function of
$\boldsymbol{\delta x}_{0}$. In general, there is no relation between the
forward and backward LEs. However, in this work we are considering systems
with a preserved measure $\mu$; then, for almost every initial position
$\boldsymbol{x}_{0}$ with respect to $\mu$, the forward and backward LEs,
$\lambda_{j}^{+}$ and $\lambda_{j}^{-}$, are opposites, and there exist
initial displacement vectors $\boldsymbol{\delta x}_{0}$ giving rise to these
opposite LEs [9]. This implies that

$\lVert\boldsymbol{\delta x}(t)\rVert=e^{\lambda_{j}t+o(t)}\,,$ (5)

where $o(t)/t$ vanishes for both $t\to+\infty$ and $t\to-\infty$ and

$\lambda_{j}:=\lambda_{j}^{+}=-\lambda_{j}^{-}\,.$ (6)

The displacement vectors giving rise to Eq. (5) are called covariant (or
characteristic) Lyapunov vectors with LEs $\lambda_{j}$ defined by Eq. (6).
Consistently with these results we will adopt the following general
definition.

######  Definition 1 (covariant Lyapunov vector).

Let $\boldsymbol{\delta x}(t)$ be the solution of Eq. (2) with initial data
$\boldsymbol{\delta x}(0)=\boldsymbol{v}$, where $\boldsymbol{v}$ is a tangent
vector at a point $\boldsymbol{x}_{0}\in X$. We say that $\boldsymbol{v}$ is a
“covariant Lyapunov vector” (CLV) at the point $\boldsymbol{x}_{0}$ if

$\lim_{t\to-\infty}\frac{\ln\lVert\boldsymbol{\delta
x}(t)\rVert}{t}=\lim_{t\to+\infty}\frac{\ln\lVert\boldsymbol{\delta
x}(t)\rVert}{t}\,.$ (7)

The common value $\lambda$ of the two above limits is called the “Lyapunov
exponent of the CLV $\boldsymbol{v}$”.

According to the above definition, at a given point $\boldsymbol{x}_{0}$ of
the phase space, the CLVs corresponding to each $\lambda_{j}$ form a linear
space of dimension $\nu_{j}$, called the multiplicity of the LE $\lambda_{j}$,
and the sum of the multiplicities of all the LEs equals the dimension $n$ of
the space $X$. The linear space corresponding to each LE $\lambda_{j}$ is
independent of the particular metric which is used on the space $X$. Hence
these spaces provide a splitting of the tangent space at almost every point of
$X$ (the so called “Oseledec splitting”) which represents an intrinsic
characterization of the dynamical system. The possible presence of
multiplicities larger than 1, called degeneration, is often neglected in the
literature, e.g. in Ref. [10]. In the absence of degeneration, one might say
that there is a single CLV $\boldsymbol{v}_{j}$, defined up to an arbitrary
scalar factor, for each LE $\lambda_{j}$, with $1\leq j\leq n$.

Various methods have been recently developed with the aim of calculating the
CLVs; a discussion can be found in Ref. [11]. They are divided into the so-
called “static” [12, 13] and “dynamic” methods [10]. With the aid of such
methods, the CLVs have been evaluated and used as a diagnostic tool in various
systems, e.g. in spatially extended dynamical systems exhibiting chaos [14],
with hyperbolic chaotic dynamics [15], in large chaotic systems consisting in
globally coupled maps [16], in stationary system out of equilibrium [17], and
in the phase synchronization transition of chaotic oscillators [18].

It is easy to see that, if $\boldsymbol{v}=\boldsymbol{\delta x}_{0}$ is a CLV
at a point $\boldsymbol{x}_{0}$ with LE $\lambda$, then the vector
$\boldsymbol{\delta x}(t)$, evolving from $\boldsymbol{\delta x}_{0}$
according to Eqs. (1) and (2), is for any $t$ a CLV at the point
$\boldsymbol{x}(t)$ with the same LE $\lambda$. This property is expressed by
saying that the CLVs are “invariant under the linearized flow” [12], or that
“CLVs are mapped to other CLVs by the linear propagator along trajectories”
[19]. This fact implies that, for a given real number $\lambda$, the set
$D\subseteq X$, of all the points at which $\lambda$ is a LE, is invariant
under the evolution of the system.

If $\nu$ is the multiplicity of the LE $\lambda$, one may arbitrarily choose,
at each point of $D$, $\nu$ CLVs which form a basis of the corresponding
linear subspace. This leads us to introduce the following definition.

######  Definition 2 (covariant Lyapunov field).

Let $\boldsymbol{F}(\boldsymbol{x})$ be a vector field on the Riemannian
manifold $X$, and let $D\subseteq X$ be an invariant set with respect to the
evolution generated by $\boldsymbol{F}$. Let $\boldsymbol{v}(\boldsymbol{x})$
be a vector field, continuous and differentiable on $D$, such that
$\lVert\boldsymbol{v}\rVert$ and $\lVert\boldsymbol{v}\rVert^{-1}$ are both
limited on $D$. If, at every point $\boldsymbol{x}\in D$,
$\boldsymbol{v}(\boldsymbol{x})$ is a CLV with LE $\lambda$ independent of
$\boldsymbol{x}$, then we say that $\boldsymbol{v}$ is a “covariant Lyapunov
field” (CLF) on the domain $D$ with LE $\lambda$.

We then see that a suitable choice of $\nu$ linearly independent CLVs, at each
point of $D$, determines $\nu$ linearly independent CLFs associated with the
LE $\lambda$. The aim of this letter is to show some important properties of
such CLFs. Note that, since the CLVs, corresponding to a given LE, form a
linear space at each point, a bound on $\lVert\boldsymbol{v}\rVert$ and
$\lVert\boldsymbol{v}\rVert^{-1}$ can always be imposed by multiplying, if
necessary, the vector field by a suitable scalar function.

The following lemma shows that any CLF, associated with a nondegenerate LE, is
the solution of a particular partial differential equation of first order
along any trajectory contained in the domain $D$.

######  Lemma 1.

Let the vector field $\boldsymbol{F}(\boldsymbol{x})$ generate a flux on the
Riemannian manifold $X$ according to Eq. (1), and let
$\boldsymbol{v}(\boldsymbol{x})$ be a CLF, on an invariant domain $D\subseteq
X$, corresponding to a nondegenerate LE $\lambda$. Then there exists on $D$ a
scalar function $b$ such that, for any trajectory $\gamma\subseteq D$, defined
by a function $\boldsymbol{x}(t)$ satisfying Eq. (1) for $-\infty<t<+\infty$,
the differential equation

$\frac{\mathrm{d}}{\mathrm{d}t}v^{\mu}(\boldsymbol{x}(t))=\partial_{\nu}F^{\mu}\left[\boldsymbol{x}(t)\right]v^{\nu}(\boldsymbol{x}(t))-b(\boldsymbol{x}(t))v^{\mu}(\boldsymbol{x}(t))$
(8)

holds on $\gamma$, and the time average of $b$ over $\gamma$ is equal to the
LE $\lambda$:

$\lim_{t\to\pm\infty}\frac{\int_{0}^{t}b\left[\boldsymbol{x}(t^{\prime})\right]{\mathrm{d}}t^{\prime}}{t}=\lambda\,.$
(9)

###### Proof.

For a given trajectory $\gamma\subseteq D$, let $\boldsymbol{\delta x}(t)$ be
the solution of Eq. (2) with initial condition $\boldsymbol{\delta
x}(0)=\boldsymbol{v}(\boldsymbol{x}_{0})$, where
$\boldsymbol{x}_{0}=\boldsymbol{x}(0)$. As we have already recalled, the
property of being a CLV is maintained by the linearized flow, so
$\boldsymbol{\delta x}(t)$ is for any $t$ a CLV with LE $\lambda$ at the point
$\boldsymbol{x}(t)$. On the other hand, by hypothesis also
$\boldsymbol{v}(\boldsymbol{x}(t))$ is for any $t$ a CLV with LE $\lambda$.
Since $\lambda$ is assumed to be a nondegenerate LE, the corresponding CLVs
form at all the points of $D$ a 1-dimensional linear space, so there exists a
nonvanishing real function $a(t)$ such that

$\boldsymbol{\delta x}(t)=a(t)\boldsymbol{v}\left[\boldsymbol{x}(t)\right]$
(10)

at all the points of $\gamma$. Since $\boldsymbol{\delta x}(t)$ and
$\boldsymbol{v}\left[\boldsymbol{x}(t)\right]$ are continuous functions of $t$
and $a(0)=1$, we get $a(t)>0$ for any $t$.

By substituting Eq. (10) into Eq. (2) we obtain

$\frac{\mathrm{d}}{\mathrm{d}t}\left\\{a(t)v^{\mu}\left[\boldsymbol{x}(t)\right]\right\\}=\partial_{\nu}F^{\mu}\left[\boldsymbol{x}(t)\right]a(t)v^{\nu}\left[\boldsymbol{x}(t)\right]$

and so

$\frac{\mathrm{d}}{\mathrm{d}t}v^{\mu}(\boldsymbol{x}(t))=\partial_{\nu}F^{\mu}\left[\boldsymbol{x}(t)\right]v^{\nu}\left[\boldsymbol{x}(t)\right]-\frac{\dot{a}(t)}{a(t)}v^{\mu}\left[\boldsymbol{x}(t)\right]\,.$

This shows that Eq. (8) holds with

$b(\boldsymbol{x}(t))=\frac{\dot{a}(t)}{a(t)}=\frac{\mathrm{d}}{\mathrm{d}t}\ln
a(t)\,.$ (11)

The fact that $\boldsymbol{v}(\boldsymbol{x}_{0})$ is a CLV with LE $\lambda$
means that

$\displaystyle\lambda$
$\displaystyle=\lim_{t\to\pm\infty}\frac{\ln\lVert\boldsymbol{\delta
x}(t)\rVert}{t}$ $\displaystyle=\lim_{t\to\pm\infty}\left(\frac{\ln
a(t)}{t}+\frac{\ln\lVert\boldsymbol{v}\left[\boldsymbol{x}(t)\right]\rVert}{t}\right)\,.$
(12)

From Eq. (11) and from $a(0)=1$ we get

$\ln
a(t)=\int_{0}^{t}b\left[\boldsymbol{x}(t^{\prime})\right]{\mathrm{d}}t^{\prime}\,.$
(13)

Moreover, definition 2 implies that both $\lVert\boldsymbol{v}\rVert$ and
$\lVert\boldsymbol{v}\rVert^{-1}$ are limited on $\gamma$, so that

$\lim_{t\to\pm\infty}\frac{\ln\lVert\boldsymbol{v}\left[\boldsymbol{x}(t)\right]\rVert}{t}=0\,.$

From Eq. (Proof.) one then obtains Eq. (9). ∎

For an ergodic system, the set of LEs is the same at almost all points of $X$
with respect to the preserved measure $\mu$ [9]. For such systems one can then
expect that there exist CLFs defined almost everywhere on $X$, and so the
result, which was proved in lemma 1 for individual trajectories, can be
generalized in the following way.

######  Proposition 1.

Let the vector field $\boldsymbol{F}$ generate an ergodic flux on a manifold
$X$, and let $\mu$ be a measure on $X$ preserved by the flux, such that
$\mu(X)<+\infty$. Let $\boldsymbol{v}$ be a CLF, corresponding to a
nondegenerate LE $\lambda$, on an open invariant domain $D\subseteq X$ such
that $\mu(D)=\mu(X)$. Then there exists a scalar field $b$ such that the
partial differential equation

$\left[\boldsymbol{v},\boldsymbol{F}\right]=b\boldsymbol{v}$ (14)

holds on $D$, where $\left[\cdot,\cdot\right]$ is the commutator:

$\left[\boldsymbol{v},\boldsymbol{F}\right]^{\mu}=v^{\nu}\partial_{\nu}F^{\mu}-F^{\nu}\partial_{\nu}v^{\mu}\,.$
(15)

Moreover

$\lambda=\left<b\right>\,,$ (16)

where $\left<b\right>$ denotes the average of $b$ over the manifold $X$:

$\left<b\right>=\frac{1}{\mu(X)}\int_{X}\mathrm{d}\mu(\boldsymbol{x})\,b(\boldsymbol{x})\,.$
(17)

###### Proof.

From lemma 1 it follows that there exists a scalar field $b$ on $D$ such that
Eq. (8) holds along any trajectory contained in $D$. But since the domain $D$
is open, we have

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}v^{\mu}(\boldsymbol{x}(t))$
$\displaystyle=\dot{x}^{\nu}(t)\partial_{\nu}v^{\mu}(\boldsymbol{x}(t))$
$\displaystyle=F^{\nu}\left[\boldsymbol{x}(t)\right]\partial_{\nu}v^{\mu}(\boldsymbol{x}(t))\,.$

By substituting this equation into Eq. (8) we obtain that Eq. (14) holds at
all points of $D$.

In addition, the ergodicity implies that the time average $\lambda$ of the
function $b$ along a trajectory equals the average of $b$ over the phase
space. Hence Eq. (16) follows from Eq. (9). ∎

The following proposition can be considered in some respect as the inverse of
the previous one. It shows in fact that in an ergodic system, under quite
general hypotheses, the fact of satisfying Eq. (14) is a sufficient condition
for a vector field $\boldsymbol{v}$ in order to be a CLF.

######  Proposition 2.

Let the vector field $\boldsymbol{F}$ generate an ergodic flux on a manifold
$X$, and let $\mu$ be a measure on $X$ preserved by the flux, such that
$\mu(X)<+\infty$. Let $\boldsymbol{v}$ and $b$ be respectively a nonvanishing
vector field and a scalar function satisfying Eq. (14) on an open invariant
domain $D\subseteq X$, such that $\mu(D)=\mu(X)$. Let also suppose that both
$\lVert\boldsymbol{v}\rVert$ and $\lVert\boldsymbol{v}\rVert^{-1}$ are limited
on $D$, and that the function $b$ is integrable over $X$, so that we can
define its average according to Eq. (17). Then $\boldsymbol{v}$ is a CLF, with
LE $\lambda=\left<b\right>$, on an invariant domain $D^{\prime}\subseteq D$
such that $\mu(D^{\prime})=\mu(X)$.

###### Proof.

Let us take a point $\boldsymbol{x}_{0}\in D$ and let $\boldsymbol{x}(t)$ be
the corresponding trajectory, i.e. the solution of Eq. (1) with initial
condition $\boldsymbol{x}(0)=\boldsymbol{x}_{0}$. We have

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}v^{\mu}\left[\boldsymbol{x}(t)\right]$
$\displaystyle=F^{\nu}\partial_{\nu}v^{\mu}\big{|}_{\boldsymbol{x}(t)}$
$\displaystyle=\big{[}[F,v]^{\mu}+v^{\nu}\partial_{\nu}F^{\mu}\big{]}_{\boldsymbol{x}(t)}$
$\displaystyle=\big{[}-bv^{\mu}+v^{\nu}\partial_{\nu}F^{\mu}\big{]}_{\boldsymbol{x}(t)}\,,$

where Eq. (14) has been applied in the last equality. Hence, for any
nonvanishing scalar function $a(t)$ we have

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\left\\{a(t)v^{\mu}\left[\boldsymbol{x}(t)\right]\right\\}=$
$\displaystyle\
\dot{a}(t)v^{\mu}\left[\boldsymbol{x}(t)\right]+a(t)\frac{\mathrm{d}}{\mathrm{d}t}v^{\mu}\left[\boldsymbol{x}(t)\right]$
$\displaystyle=$
$\displaystyle\left\\{\dot{a}(t)-a(t)b\left[\boldsymbol{x}(t)\right]\right\\}v^{\mu}\left[\boldsymbol{x}(t)\right]$
$\displaystyle+a(t)v^{\nu}\left[\boldsymbol{x}(t)\right]\partial_{\nu}F^{\mu}\left[\boldsymbol{x}(t)\right]\,.$

We then see that, by putting $\boldsymbol{\delta
x}(t)=a(t)\boldsymbol{v}\left[\boldsymbol{x}(t)\right]$, we obtain the
solution of Eq. (2) with initial condition $\boldsymbol{\delta
x}(0)=\boldsymbol{v}(\boldsymbol{x}_{0})$, provided that $a(0)=1$ and $a(t)$
satisfies Eq. (11). Hence $\ln a(t)$ is given by Eq. (13) and we get

$\displaystyle\lim_{t\to+\infty}\frac{\ln\lVert\boldsymbol{\delta
x}(t)\rVert}{t}$ $\displaystyle=\lim_{t\to+\infty}\left(\frac{\ln
a(t)}{t}+\frac{\ln\lVert\boldsymbol{v}\left[\boldsymbol{x}(t)\right]\rVert}{t}\right)$
$\displaystyle=\lim_{t\to+\infty}\frac{\int_{0}^{t}b\left[\boldsymbol{x}(t^{\prime})\right]{\mathrm{d}}t^{\prime}}{t}\,,$
(18)

where for the last equality we have used the hypothesis that both
$\lVert\boldsymbol{v}\rVert$ and $\lVert\boldsymbol{v}\rVert^{-1}$ are limited
on $D$.

The last member of Eq. (Proof.) represents the time average of the function
$b(\boldsymbol{x})$ on the considered trajectory for positive times. Thanks to
the ergodicity, for almost all the points $\boldsymbol{x}_{0}\in D$ this
average equals the average $\left<b\right>$ of $b$ over the whole manifold
$X$.

By analyzing in a similar way the limit for $t\to-\infty$, we obtain that
there exists a subset $D^{\prime}\subseteq D$, with $\mu(D^{\prime})=\mu(X)$,
such that

$\lim_{t\to-\infty}\frac{\ln\lVert\boldsymbol{\delta
x}(t)\rVert}{t}=\lim_{t\to+\infty}\frac{\ln\lVert\boldsymbol{\delta
x}(t)\rVert}{t}=\langle b\rangle$

for all the points $\boldsymbol{x}_{0}\in D^{\prime}$. According to definition
1, this means that $\boldsymbol{v}(\boldsymbol{x}_{0})$ is a CLV at
$\boldsymbol{x}_{0}$ with LE $\lambda=\langle b\rangle$. The set $D^{\prime}$
is obviously invariant under the evolution of the system so, by applying
definition 2, we conclude that $\boldsymbol{v}$ is a CLF on $D^{\prime}$ with
LE $\lambda=\langle b\rangle$. ∎

Propositions 1 and 2 together imply the remarkable fact that, if the system is
ergodic and $\lambda$ is a nondegenerate LE, then a vector field
$\boldsymbol{v}$ is a CLF with LE $\lambda$ if and only if it satisfies Eq.
(14) almost everywhere on $X$. Note that this is a global condition on the
vector field $\boldsymbol{v}$. It is in fact easy to see that a local solution
of the first order differential equation (14) can be obtained after
arbitrarily assigning the vector $\boldsymbol{v}$ on a $(n-1)$-dimensional
surface transversal to the flux generated by $\boldsymbol{F}$. So the property
of being a CLF only arises if one requires that Eq. (14) be satisfied almost
everywhere on $X$, i.e. in a global sense.

The fact that CLVs are determined up to an arbitrary scalar factor implies
that, if $\boldsymbol{v}(\boldsymbol{x})$ is a CLF, then the same is true also
for the vector field

$\boldsymbol{v}^{\prime}(\boldsymbol{x})=c(\boldsymbol{x})\boldsymbol{v}(\boldsymbol{x})\,,$
(19)

where $c(\boldsymbol{x})$ is an arbitrary nonvanishing smooth scalar function.
If $\boldsymbol{v}$ satisfies Eq. (14), it is immediate to see that
$\boldsymbol{v}^{\prime}$ satisfies the equation

$\left[\boldsymbol{v}^{\prime},\boldsymbol{F}\right]=b^{\prime}\boldsymbol{v}^{\prime}$
(20)

with

$b^{\prime}(\boldsymbol{x})=b(\boldsymbol{x})-F^{\mu}(\boldsymbol{x})\partial_{\mu}\ln|c(\boldsymbol{x})|\,.$
(21)

Since the two Eqs. (14) and (20) are formally identical, we can say that Eq.
(14) is invariant under the local “gauge transformation” expressed by Eqs.
(19) and (21).

Since the function $c$ nowhere vanishes, by continuity it has constant sign
over the domain $D$ of the CLF. Assuming that the sign is positive, we can
write $c(\boldsymbol{x})=e^{\phi(\boldsymbol{x})}$, where $\phi$ is an
arbitrary smooth scalar function. Introducing the Lie derivative
$\mathcal{L}_{\boldsymbol{F}}\varphi:=F^{\mu}\partial_{\mu}\varphi$ of
$\varphi$ with respect to $\boldsymbol{F}$, we can write the transformation
given by Eqs. (19) and (21) in the form

$\left\\{\begin{array}[]{l}\boldsymbol{v}(\boldsymbol{x})\mapsto
e^{\varphi(\boldsymbol{x})}\,\boldsymbol{v}(\boldsymbol{x})\\\
b(\boldsymbol{x})\mapsto
b(\boldsymbol{x})-\mathcal{L}_{\boldsymbol{F}}\varphi(\boldsymbol{x})\,.\end{array}\right.$
(22)

From a mathematical point of view, such a gauge invariance recalls that of
field theories in fundamental physics. For instance, in quantum
electrodynamics, the fact the wavefunction $\psi$ is defined at each space-
time point up to an arbitrary phase factor, implies that Dirac equation

$\gamma^{\mu}\left[i\partial_{\mu}-eA_{\mu}(x)\right]\psi(x)-m\psi(x)=0$

is invariant under the local gauge transformation

$\left\\{\begin{array}[]{l}\psi(x)\mapsto e^{ie\alpha(x)}\,\psi(x)\\\
A^{\mu}(x)\mapsto A^{\mu}(x)-\partial_{\mu}\alpha(x)\,,\end{array}\right.$
(23)

where $x$ stands for the four space-time coordinates and $\alpha(x)$ is an
arbitrary real scalar function [20].

The analogy between Eqs. (22) and (23) is obvious. The transformation on the
four-vector potential $A^{\mu}$, given by Eq. (23), does not alter the value
of the physically relevant electromagnetic field tensor
$F^{\mu\nu}=\partial^{\mu}A^{\nu}-\partial^{\nu}A^{\mu}$. In a similar way,
the transformation on the scalar function $b$, given by Eq. (22), does not
alter the physically relevant value of the LE $\lambda$, since it is easy to
see that $\langle b^{\prime}\rangle=\langle b\rangle$.

If a metric tensor has been defined over the manifold $X$, e.g. the euclidean
tensor in the chosen system of coordinates, then a convenient gauge can be
$\lVert\boldsymbol{v}\rVert=1$, i.e. the request that the CLVs are normalized
at every point. However, since dynamical systems in general are not equipped
with an intrinsically defined metric, no preferential gauge is available in
most cases.

It is worth remarking that the definition 1 of CLV and of the corresponding
LE, similarly to other definitions of Lyapunov vectors and exponents adopted
in the literature, is explicitly based on the existence of a norm of vectors,
as shown by Eq. (7). The same is then true also for the definition 2 of CLF.
Despite this fact, as we have already pointed out, both the property of being
a CLV, and the value of the LE, are actually independent of the choice of a
particular metric tensor on the space $X$.

It is therefore interesting to note that propositions 1 and 2 provide the
possibility of an alternative definition of CLF, and of the corresponding LE,
which does not mention at all the existence of a norm. One could in fact
define as CLF any vector field satisfying Eq. (14), and define its LE as
$\lambda=\langle b\rangle$. In the case of nondegenerate CLFs in ergodic
systems, such a definition would be equivalent to definition 2. In the case of
a LE with degeneracy $\nu>1$, one could conjecture, under suitable hypotheses,
the existence of $\nu$ linearly independent vector fields, each one satisfying
an equation of the form of Eq. (14).

The two authors equally contributed to this work.

## References

  * [1] Arkady Pikovsky and Antonio Politi.  Lyapunov Exponents: A Tool to Explore Complex Dynamics.  Cambridge University Press, Cambridge, 2016.
  * [2] V. I. Oseledec.  A multiplicative ergodic theorem. Lyapunov characteristic numbers for dynamical systems.  Trans. Moscow Math. Soc., 19, 1968.
  * [3] J.-P. Eckmann and D. Ruelle.  Ergodic theory of chaos and strange attractors.  Rev. Mod. Phys., 57, 1985.
  * [4] Alan Wolf, Jack B. Swift, Harry L. Swinney, and John A. Vastano.  Determining lyapunov exponents from a time series.  Physica D: Nonlinear Phenomena, 16(3):285–317, 1985.
  * [5] Y. B. Pesin.  Characteristic Lyapunov exponents and smooth ergodic theory.  Russian Mathematical Surveys, 32(4):55, 1977.
  * [6] Paul Frederickson, James L. Kaplan, Ellen D. Yorke, and James A. Yorke.  The Liapunov dimension of strange attractors.  Journal of Differential Equations, 49(2):185–207, 1983.
  * [7] G. Benettin, L. Galgani, A. Giorgilli, and J. M. Strelcyn.  Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. Part 1: Theory.  Meccanica, 15:9–20, 1980.
  * [8] Ippei Shimada and Tomomasa Nagashima.  A numerical approach to ergodic problem of dissipative dynamical systems.  Progress of Theoretical Physics, 61(6):1605–1616, 1979.
  * [9] D. Ruelle.  Ergodic theory of differentiable dynamical systems.  Publ. Math. de l’I.H.É.S., 50, 1979.
  * [10] F. Ginelli, P. Poggi, A. Turchi, H. Chaté, R. Livi, and A. Politi.  Characterizing dynamics with covariant Lyapunov vectors.  Phys. Rev. Lett., 99:130601, 2007.
  * [11] F. Ginelli, H. Chaté, R. Livi, and A. Politi.  Covariant Lyapunov vectors.  J. Phys. A: Mathematical and theoretical, 46(25), 2013.
  * [12] C. L. Wolfe and R. M. Samelson.  An efficient method for recovering Lyapunov vectors from singular vectors.  Tellus A: Dynamic Meteorology and Oceanography, 59(3):355–366, 2007\.
  * [13] P. V. Kuptsov and U. Parlitz.  Theory and computation of covariant Lyapunov vectors.  J. of Nonlinear Science, 22(5), 2012.
  * [14] Ivan G. Szendro, Diego Pazó, Miguel A. Rodríguez, and Juan M. López.  Spatiotemporal structure of lyapunov vectors in chaotic coupled-map lattices.  Phys. Rev. E, 76:025202, 2007.
  * [15] Pavel V. Kuptsov and Sergey P. Kuznetsov.  Violation of hyperbolicity in a diffusive medium with local hyperbolic attractor.  Phys. Rev. E, 80:016205, 2009.
  * [16] Kazumasa A. Takeuchi, Francesco Ginelli, and Hugues Chaté.  Lyapunov analysis captures the collective dynamics of large chaotic systems.  Phys. Rev. Lett., 103:154103, 2009.
  * [17] Hadrien Bosetti, Harald A. Posch, Christoph Dellago, and William G. Hoover.  Time-reversal symmetry and covariant Lyapunov vectors for simple particle models in and out of thermal equilibrium.  Phys. Rev. E, 82:046218, 2010.
  * [18] Katsuya Ouchi, Takehiko Horita, and Tomoji Yamada.  Characterizing the phase synchronization transition of chaotic oscillators.  Phys. Rev. E, 83:046202, 2011.
  * [19] F. Noethen.  Computing covariant lyapunov vectors – a convergence analysis of Ginelli’s algorithm.  Ph. D. thesis, University of Hamburg, 2019.
  * [20] Claude Itzykson and Jean-Bernard Zuber.  Quantum Field Theory.  McGraw-Hill, 1980.


11institutetext: Department of Engineering Design, Indian Institute of
Technology Madras, India
11email<EMAIL_ADDRESS>
# Structurally aware bidirectional unpaired image to image translation between
CT and MR

Vismay Agrawal  Avinash Kori  Vikas Kumar Anand  Ganapathy Krishnamurthi

###### Abstract

Magnetic Resonance (MR) Imaging and Computed Tomography (CT) are the primary
diagnostic imaging modalities quite frequently used for surgical planning and
analysis. A general problem with medical imaging is that the acquisition
process is quite expensive and time-consuming. Deep learning techniques like
generative adversarial networks (GANs) can help us to leverage the possibility
of an image to image translation between multiple imaging modalities, which in
turn helps in saving time and cost. These techniques will help to conduct
surgical planning under CT with the feedback of MRI information. While
previous studies have shown paired and unpaired image synthesis from MR to CT,
image synthesis from CT to MR still remains a challenge, since it involves the
addition of extra tissue information. In this manuscript, we have implemented
two different variations of Generative Adversarial Networks exploiting the
cycling consistency and structural similarity between both CT and MR image
modalities on a pelvis dataset, thus facilitating a bidirectional exchange of
content and style between these image modalities. The proposed GANs translate
the input medical images by different mechanisms, and hence generated images
not only appears realistic but also performs well across various comparison
metrics, and these images have also been cross verified with a radiologist.
The radiologist verification has shown that slight variations in generated MR
and CT images may not be exactly the same as their true counterpart but it can
be used for medical purposes.

###### Keywords:

MRI CT GAN CNN

##  1 Introduction

Computer Tomography (CT) and Magnetic Resonance (MR) Imaging are the two most
frequently used imaging modalities for medical diagnosis and surgical
planning. While CT provides contrast information for bone studies, MR imaging
provides information for soft tissue contrast. In this manuscript we propose a
Generative Adversarial Networks (GANs) [4] based technique to convert CT
images to realistic looking MR images. As the physics of acquisition in the
case of MR and CT differ significantly, developing techniques for direct
transformation of images from MR space to CT space is not a trivial problem.
But recent advancements in data-driven modeling, specifically GANs have
provided realistic looking results on image to image translation by separating
content and style from an image. We make use of similar techniques to separate
style and content from MR and CT images and interchange them. To overcome the
challenges in obtaining paired data, we conducted our experiments along the
lines of cycleGAN [14] based approach, where the models were trained on
unpaired data to maintain the cyclic consistency between different imaging
modalities. The use of unpaired data drastically increases the amount of
input-output pairs used in translational networks. It also assures that the
network is not overfitting itself on the given data but is learning the visual
properties of different image modalities. We have also played with loss
functions and have shown that using structural similarity in loss function can
generate realistically looking images which can be used in looking at bone
structure and fractures in a better way. These findings were also cross
verified by a radiologist.

Previously, Conditional GANs [6] have been used with U-net [10] based
architecture as a generator and convolutional PatchGAN classifier as a
discriminator to learn the mapping from input images to output images. In some
instances, CNNs and GANs have also been implemented in MR to CT translation
[8] [9]. cycleGANs have been used to perform a bidirectional unpaired image to
image translation between Cardiac CT and MR data [2] and it has been shown
that the generated data can be used to increase the accuracy of the network in
segmentation tasks. Wolterink et al. [12] has shown that cycleGAN trained on
unpaired data was able to outperform the model trained on paired data for MR
to CT translation. Jin et al. [7] study was focused on CT to MR translation
which consisted of dual cycle-consistency loss using paired and unpaired
training data, but due to the use of paired data in the training process, the
method can’t be used in the situation with a lack of paired data. With all
these, we can see that although many research work has been done on brain and
cardiac medical scans, there aren’t many successful works with pelvis scans.
Also, we have observed that most of the research was focused on obtaining CT
images from their MR counterpart, with little or no details about bi-
directional translations.

##  2 Materials and Methods

###  2.1 Data

Our data consisted of human pelvis unpaired (taken from multiple patients) MR
and CT 3D volumes. T2 weighted images were obtained from local hospital
(Apollo Speciality Hospital, Chennai, TamilNadu, India) with appropriate
ethical clearance which were acquired using a 1.5 Tesla system (Achieva,
Philips Medical System), while CT volumes we extracted from Liver Tumour
Segmentation (LiTS) challenge database [3], under the guidance of expert
radiologists. In total, we have about 55 MR and CT volumes, average with about
30 slices in MR and about 80 slices in CT volumes. We extracted axial slices
from these MR and CT volumes to train the network. The data split used in
training and testing is shown in table 1.

Table 1: Dataset Details | Training Set | Test Set | Volume Dimension |
---|---|---|---|---
CT Volume | 50 | 5 | (512, 512, 80) |
MR Volume | 50 | 5 | (768, 768, 30) |

###  2.2 Pre-processing of data

As the number of slices was different in MR and CT, to maintain an equal
number of slices we resampled all volumes to have 80 slices on an average.
Each image slice was then normalized using the min-max normalization method.
The obtained image was resized to the dimension (286, 286) using bicubic
interpolation and subsequently, it was cropped to a dimension (256, 256) from
a random location. Data augmentation was performed by flipping and rotating
the input image at random.

###  2.3 Transformer Network

Our pipeline consisted of two identical Generator networks for CT to MR and MR
to CT image to image translation, i.e. $G_{CT}:I_{MR}\rightarrow{}I_{CT}$ and
$G_{MR}:I_{CT}\rightarrow{}I_{MR}$. During training, one of the networks
performs the required translation and the other one tries to undo the same. We
have used ResNet [5] as a generator network with 9 resblocks. For each
generator, we have a corresponding discriminator network, $D_{CT}$ and
$D_{MR}$, which aims to discriminate between translated and original, CT and
MR images respectively. The pipeline of our network is given in Fig. 1

Figure 1: Schematics of our method used inter sequence MRI images

###  2.4 Loss function

In the two proposed networks, the first one (cycleGAN) uses loss function as
proposed in [13] and in the second network (cycleGAN-SSIM) we have included
the structural similarity loss function [11] on top of all the other losses.
The loss functions can be divided into multiple parts as stated below.

####  2.4.1 GAN Loss

$\displaystyle\mathcal{L}_{GAN}(D_{MR},G_{MR},I_{CT})=\ $
$\displaystyle\mathbb{E}_{I_{CT}\sim
Data_{CT}}[(D_{MR}(G_{MR}(I_{CT}))-1)^{2}]$
$\displaystyle+\mathbb{E}_{I_{MR}\sim
Data_{MR}}[(D_{CT}(G_{CT}(I_{MR}))-1)^{2}]$

The GAN loss aims to minimize the possibility of the translated image being
recognized as a fake image by the discriminator.

####  2.4.2 Cycle Loss

$\displaystyle\mathcal{L}_{cycle}=\ $ $\displaystyle\mathbb{E}_{I_{MR}\sim
Data_{MR}}[||G_{MR}(G_{CT}(I_{MR}))-I_{MR})||_{1}]$
$\displaystyle+\mathbb{E}_{I_{CT}\sim
Data_{CT}}[||G_{CT}(G_{MR}(I_{CT}))-I_{CT})||_{1}]$

The cycle loss function ensures the cyclic consistency between the input and
generated images.

####  2.4.3 Identity Loss

$\displaystyle\mathcal{L}_{identity}=\ $ $\displaystyle\mathbb{E}_{I_{CT}\sim
Data_{CT}}[||(G_{CT}(I_{CT}))-I_{CT})||_{1}]$
$\displaystyle+\mathbb{E}_{I_{MR}\sim
Data_{MR}}[||(G_{MR}(I_{MR}))-I_{MR})||_{1}]$

The root of the identity loss lies in the fact that the generator should be
aware of the input image modality, and the generator shouldn’t translate the
input image if it corresponds to same imaging modality as the output of the
generator.

####  2.4.4 Structural Similarity (SSIM) Loss

$\displaystyle\mu_{MR/CT}$ $\displaystyle=\mathbb{E}_{I_{MR/CT}\sim
Data_{MR/CT}}[G_{CT/MR}(I_{MR/CT})]$ $\displaystyle\sigma_{MR/CT}$
$\displaystyle=\mathbb{E}_{I_{MR/CT}\sim
Data_{MR/CT}}[(G_{CT/MR}(I_{MR/CT}))^{2}]$ $\displaystyle\sigma_{MR,CT}$
$\displaystyle=\mathbb{E}_{I_{CT}\sim Data_{CT},I_{MR}\sim
Data_{MR}}[G_{MR}(I_{CT})G_{CT}(I_{MR})]$ $\displaystyle\mathcal{L}_{SSIM}$
$\displaystyle=1-\
\frac{(\mu_{CT}\mu_{MR}+c_{1})(2\sigma_{MR,CT}+c_{2})}{(\mu_{MR}^{2}+\mu_{CT}^{2}+c_{1})(\sigma_{MR}^{2}+\sigma_{CT}^{2}+c_{2})}$

Here, we chose $c_{1}$ and $c_{2}$ to be 0.0001 and 0.009 respectively. SSIM
helps in maintaining the structural integrity between inter-domain images.

####  2.4.5 Generator Net Loss

$\displaystyle\mathcal{L}_{G_{net}}=\ $
$\displaystyle\mathcal{L}_{GAN}+\lambda_{cyc}*\mathcal{L}_{cycle}+\lambda_{id}*\mathcal{L}_{identity}+\lambda_{ssim}*\mathcal{L}_{ssim}$

Net generator loss is defined as a superposition of all the aforementioned
losses namely, GAN loss, cycle loss, identity and SSIM loss.

####  2.4.6 Discriminator Loss

$\displaystyle\mathcal{L}_{dis_{MR/CT}}=$ $\displaystyle\
\mathbb{E}_{I_{CT/MR}\sim Data_{CT/MR}}[\ (D_{MR/CT}(I_{CT/MR})-1)^{2}\ ]$
$\displaystyle+\mathbb{E}_{I_{CT/MR}\sim Data_{CT/MR}}[\
(D_{MR/CT}(G_{MR/CT}(I_{CT/MR}))^{2}\ ]$

As our aim is to predict label 1 for a real imaging modality and label 0 for a
translated one, the loss function is the distance between the prediction and
the real label.

###  2.5 Comparison Metrics

Quantitative analysis of the generated images was performed using FID index
(Fréchet Inception Distance) [1], MI (mutual information) [1] and SSIM index
[1].

####  2.5.1 FID

In case of FID, pretrained densenet121 ($FID(.)$) is considered as base
network, and feature embedding for real MRI/CT and corresponding fake MRI/CT
is obtained, the distance between these two embeddings are observed. If the
distribution of generator matched with the distribution of real image the
distance between feature embeddings will reduce and FID score will increase.

$\displaystyle{FID}_{CT}=$ $\displaystyle\ \mathbb{E}_{I_{CT}\sim
Data_{CT},I_{MR}\sim Data_{MR}}<FID(G_{CT}(I_{MR})),FID(I_{CT})>$

where $<.>$ indicates inner product between two vectors. FID for MR can be
calculated similarly.

####  2.5.2 SSIM

To quantitatively measure the structural integrity of the generated images we
made use structural similarity as a metric. SSIM between generated CT/MR and
input MR/CT provides us with an information how structurally close the given
pair of generated MR and CT are.

$\displaystyle{SSIM}$
$\displaystyle=\frac{(\mu_{CT}\mu_{MR}+c_{1})(2\sigma_{MR,CT}+c_{2})}{(\mu_{MR}^{2}+\mu_{CT}^{2}+c_{1})(\sigma_{MR}^{2}+\sigma_{CT}^{2}+c_{2})}$

####  2.5.3 MI

To measure the textural similarity between generated MR/CT with actual MR/CT
we made use of mutual information as a metric.

$MI(G,R)=\mathbf{E}(P_{GR}(G,R))\times
log(\frac{P_{GR}(G,R)}{P_{G}(G)P_{R}(R)})$

where $P_{GR}(G,R)$ denotes joint distribution between generated and real
distribution, $\mathbf{E}$ denotes expectation value and $P_{G}(G),P_{R}(R)$
denotes marginals.

####  2.5.4 Pixel wise accuracy (pixacc)

To gauge the performance of transformation network in generating recovered MR/
CT using real MR/ CT images, cosine similarity between them was used as a
metric, as stated below.

$pixacc_{MR/CT}=\mathbb{E}_{I_{CT}\sim Data_{CT},I_{MR}\sim
Data_{MR}}\frac{I_{MR/CT}\cdot
G_{MR/CT}(I_{CT/MR})}{||I_{MR/CT}||\cdot||G_{MR/CT}(I_{CT/MR})||}$

##  3 Results and Discussion

The performance metrics were being calculated for 200 slices and the average
results are shown in Table 2. The loss plots between cycleGAN and cycleGAN-
SSIM, given in Fig 2 clearly show the lower convergence in case of using
structural similarity loss for both generator and discriminator, we can
conclude that the network was able to exploit the structural consistencies
between different image modalities to boost its learning rate.

Table 2: Mean comparison metrics  | CT to MR Translation | MR to CT Translation
---|---|---
Model | FID | SSIM | MI | pixacc | FID | SSIM | MI | pixacc
cycleGAN | 0.193 | 0.408 | 0.273 | 0.986 | 0.177 | 0.562 | 0.332 | 0.994
cycleGAN-SSIM | 0.200 | 0.416 | 0.290 | 0.988 | 0.184 | 0.562 | 0.336 | 0.995

(a) Total Generator loss

(b) Total Discriminator loss

Figure 2: comparison of models based on loss convergence

We can see that although all the results were higher in the case of the
network trained with SSIM loss function, there aren’t any significant
differences in their scores. The use of structural similarity loss function
was preserving information about the edges of the generated images, especially
in the case of CT to MR translation (Fig 3), wherein visual inspection shows
enhanced bone details, and relatively higher structural similarity and mutual
information scores for cycleGAN-SSIM model. The MR to CT translation (Fig 4)
had similar SSIM scores for both the models, this denotes that cycleGAN was
learning to structural information in CT images. In such a case addition of
the SSIM loss function was just acting as a catalyst to increase the learning
rate of the model. High pixacc scores reflect the fact that both the models
generate realistically looking images and such synthesised images are suitable
for use in medical diagnosis. MR to CT translation had higher pixacc scores
than CT to MR translation because MR has more tissue contrasts and such
information is missing in CT, thus it isn’t possible for the network to
completely recover the lost information in latter transformation.

According to the reviews by the radiologist, in comparison to cycleGAN-SSIM
model, the images generated by cycleGAN were more similar to real and CT MR
images. Radiologist pointed out that MR images generated by cycleGAN model
would be better for routine diagnosis, whereas images generated by cycleGAN-
SSIM model would be better for special cases of diagnosis. In cycleGAN-SSIM
model, due to preservation in textural details of input CT image content, its
MR counterpart contained contrast feature of MR image with blended with
structural information of CT images. Such images can be used in parallel with
CT images for diagnosis of bone fractures and also in imaging the lung to
identify tiny nodules and calcification.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

Figure 3: CT to MR conversion. In the above image column (i) corresponds to
real CT, (ii) generated MR using cycleGAN (iii) recovery CT using cycleGAN
(iv) generated MR with cycleGAN-SSMI (v) recovery CT with cycleGAN-SSMI

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

Figure 4: MR to CT conversion. In the above image column (i) corresponds to
real MR, (ii) generated CT using cycleGAN (iii) recovery MR using cycleGAN
(iv) generated CT with cycleGAN-SSMI (v) recovery MR with cycleGAN-SSMI

##  4 Conclusion and Future work

In this work, we have shown that the cycleGAN can be used in medical image
translation tasks. Moreover, the use of structural similarity loss with
cycleGAN can boost the learning rate of GAN, images generated by such model
are more structurally sound and hence gives a better diagnosis of fractures.
In future, we plan to conduct similar experiments with entire 3D volumes, i.e
3D style and content separation from MRI and CT volumes. We would also like to
extend this work to identify tiny nodule and calcification in the lungs, as
pointed by the radiologist.

## References

  * [1] Borji, A.: Pros and cons of gan evaluation measures. Computer Vision and Image Understanding 179, 41–65 (2019)
  * [2] Chartsias, A., Joyce, T., Dharmakumar, R., Tsaftaris, S.A.: Adversarial image synthesis for unpaired multi-modal cardiac data. In: International Workshop on Simulation and Synthesis in Medical Imaging. pp. 3–13. Springer (2017)
  * [3] Christ, P.F., Elshaer, M.E.A., Ettlinger, F., Tatavarty, S., Bickel, M., Bilic, P., Rempfler, M., Armbruster, M., Hofmann, F., D’Anastasi, M., et al.: Automatic liver and lesion segmentation in ct using cascaded fully convolutional neural networks and 3d conditional random fields. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 415–423. Springer (2016)
  * [4] Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., Bengio, Y.: Generative adversarial nets. In: Advances in neural information processing systems. pp. 2672–2680 (2014)
  * [5] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 770–778 (2016)
  * [6] Isola, P., Zhu, J.Y., Zhou, T., Efros, A.A.: Image-to-image translation with conditional adversarial networks. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 1125–1134 (2017)
  * [7] Jin, C.B., Kim, H., Liu, M., Jung, W., Joo, S., Park, E., Ahn, Y.S., Han, I.H., Lee, J.I., Cui, X.: Deep ct to mr synthesis using paired and unpaired data. Sensors 19(10), 2361 (2019)
  * [8] Kaiser, B., Albarqouni, S.: Mri to ct translation with gans. arXiv preprint arXiv:1901.05259 (2019)
  * [9] Nie, D., Trullo, R., Lian, J., Petitjean, C., Ruan, S., Wang, Q., Shen, D.: Medical image synthesis with context-aware generative adversarial networks. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 417–425. Springer (2017)
  * [10] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: MICCAI (2015)
  * [11] Wang, Z., Bovik, A.C., Sheikh, H.R., Simoncelli, E.P., et al.: Image quality assessment: from error visibility to structural similarity. IEEE transactions on image processing 13(4), 600–612 (2004)
  * [12] Wolterink, J.M., Dinkla, A.M., Savenije, M.H., Seevinck, P.R., van den Berg, C.A., Išgum, I.: Deep mr to ct synthesis using unpaired data. In: International Workshop on Simulation and Synthesis in Medical Imaging. pp. 14–23. Springer (2017)
  * [13] Yang, H., Sun, J., Carass, A., Zhao, C., Lee, J., Xu, Z., Prince, J.: Unpaired Brain MR-to-CT Synthesis using a Structure-Constrained CycleGAN (9 2018), http://arxiv.org/abs/1809.04536
  * [14] Zhu, J.Y., Park, T., Isola, P., Efros, A.A.: Unpaired image-to-image translation using cycle-consistent adversarial networks. In: Proceedings of the IEEE international conference on computer vision. pp. 2223–2232 (2017)


# Entropy Transformer Networks: A Learning Approach via
Tangent Bundle Data Manifold

Pourya Shamsolmoali1, Masoumeh Zareapoor2
1East China Normal University 2Shanghai Jiao Tong University

###### Abstract

This paper focuses on an accurate and fast interpolation approach for image
transformation employed in the design of CNN architectures. Standard Spatial
Transformer Networks (STNs) use bilinear or linear interpolation as their
interpolation, with unrealistic assumptions about the underlying data
distributions, which leads to poor performance under scale variations.
Moreover, STNs do not preserve the norm of gradients in propagation due to
their dependency on sparse neighboring pixels. To address this problem, a
novel Entropy STN (ESTN) is proposed that interpolates on the data manifold
distributions. In particular, random samples are generated for each pixel in
association with the tangent space of the data manifold and construct a linear
approximation of their intensity values with an entropy regularizer to compute
the transformer parameters. A simple yet effective technique is also proposed
to normalize the non-zero values of the convolution operation, to fine-tune
the layers for gradients’ norm-regularization during training. Experiments on
challenging benchmarks show that the proposed ESTN can improve predictive
accuracy over a range of computer vision tasks, including image
reconstruction, and classification, while reducing the computational cost.

Figure 1: Data transformation and resulting distributions. Left: Data manifold
$\cal{M}$ (yellow curve), context points (blue dots), interpolate points (red
dots), Green star is the model observation. Stationary parts in target density
(green lines). Right: Plot of the transformed density (blue curve), and
histogram of the projected samples (red).

##  1 Introduction

Determining discriminant features is challenging in various computer vision
tasks, such as object detection [32] and image retrieval [21]. To detect data
relations we need an image descriptors to efficiently encode image features
that are invariant to geometric image transformations. For example, to detect
objects in real-world images, we have learned that the appearance of an object
in an image depends on different types of deformations, such as rotations and
scale changes. However, the features of an object should remain consistent
(invariant) to these transformations. Generally, CNNs do not make use of prior
learning to deal with the transformation invariance of objects across
different domains. As a substitute, CNNs use data augmentation to integrate
such priors during training. To some extent, data augmentation improves CNN
performance, but provides no assurance that transformation invariance will be
achieved for the data that has never been seen throughout the training. In
[7], the authors investigate the role of rotational equivariance in CNNs.

In addition to data augmentation, group equivariant CNNs [4], [5] have been
proposed to constantly address particular transformations of the input. For
example, in [8], a CNN is integrated with equivariant self-attention layers,
to separate the feature aggregation process into a kernel generator, and
decouple the spatial and additional geometric dimensions in the computing.
Moreover, Jaderberg et al. [11] have proposed the STN, a discernible module
that uses an attention mechanism to perform spatial manipulation and
transformation on the image by a CNN. In contrast to pooling layers, in which
the receptive fields are constant, the STN is an active process that uses the
attention mechanism to transform the feature map of the input image by
producing a transformation for the image. In [29], a rotation and translation
approach is proposed, which generalizes poorly to the composition of feature
transformations, because this model does not process the local information in
complex scenarios.

Indeed, the above models [11], [29] are based on bilinear interpolation for
calculating the values of pixels. Bilinear interpolation is local, which is
one of its weaknesses, as it only leverages the four nearest pixels of the
query. Since bilinear interpolation does not consider high-level
transformation, the performance of deep neural networks are reduced when the
image size changes. Among the techniques used for predicting missing pixels,
linear interpolation has attracted widespread attention due to its simplicity.
However, this technique cannot perform accurately because of redundant noise.
The limitation of the bilinear interpolation was addressed in [12] by
proposing Linearized-STNs (L-STNs). This approach has better performance
compared to [11], however, in L-STNs random noise is generated during
training, over-fitting often occurs, and linear regression is inefficient.

STNs are generally adopted in CNNs as pre-alignment modules for different
tasks [17]. Current STN approaches cannot deal with extensive downsampling to
increase the network capacity. To address this problem, we propose a Entropy
STN (ESTN), where the interpolation method allows to increase the capacity of
STN for learning invariance against different type of transformations. First,
a local tangent space approximation [31] is performed to estimate the data
tangent space by using the variational autoencoders (VAEs) [14], and enforce
the interpolation module to generate auxiliary samples in the neighborhood of
the query pixel. Adopting an entropy regularization allows exploring more
samples in the search space along the data manifold and improving the
convergence properties. This allows learning the geometric correlations that
lie at the intersection of the input image and the geometric manifolds. The
space of geometric manipulations is based on the transformed image (see Fig.
1), and samples are formed in the transformed domain, where their spatial
locations are determined by the transformation. Therefore, gradient features
are transformation-aware to improve network performance by identifying the
important regions of an image and detecting the changes in accordance with the
transformation.

To reconstruct an image, we rely on the content information in the local
tangent space, followed by smoothing of the orthogonal direction factors in
the local tangent space. Smoothing the orthogonal direction factors allows to
maintain the local geometries of the surface that are linearly approximated in
the tangent space. The ESTN improves the learning ability of deep neural
networks, which in turn can be used to promote any equivariance translation in
CNN architectures. To incorporate the proposed ESTN into any CNN architecture
and have smooth training, we introduce a norm-regularization method by
retaining the same variance of the weighted gradients in the layers. Our
method does not require singular value decomposition to process the non-zero
values of a CNN on the underlying data manifold. With the above modifications,
our model achieves higher accuracy with lower computation costs on different
tasks.

##  2 Related Work

Equivariant CNNs. Recently, there has been a strong interest from the machine
learning community in the design of CNN models that are equivariant in terms
of image transformation. Equivariance in terms of transformation is possible
via transforming CNN filters or manipulating feature maps [28, 4]. In [34], an
equivariant CNN to space and scaling is proposed, which is equivariance for
the regular representation of scaling and translation. Moreover, [23] proposed
an efficient equivariant network by reducing the computation of moving frames
at the input stage rather than repeating computations at all layers. In fact,
these approaches support limited transformations because the computational
cost increases linearly with the cardinality ratio of the transformation. The
equivariance approaches with regard to uninterrupted transformation are
divided into two classes: (1) those that increase the input in a guiding
foundation [5], and (2) those that assess convolutions through a comprehensive
geodetic system [9]. [33] introduces an equivariant hypergraph module that can
be adopted in various neural networks. [24] also recently proposes efficient
equivariant graph neural networks that can easily scale to higher-dimensional
spaces. In addition, a permutation-equivariant network is proposed in [13] to
capture both global and channel-wise contextual information. In [26], to build
an equivariant network, orthogonal moments of the function are used as an
effective means for encoding global invariance with respect to rotation and
translation in fully-connected layers of the network.

A new class of equivariant CNNs was proposed in [6], based on data manifold
learning and intrinsic geometry. In [27], to address the problems of learning
representations, the authors proposed a manifold neural network model that
learns from a level-set function of the constraint. However, these methods can
not realize the relative feature-space relationship between two transformed
encoded inputs. For example, the relative feature space relationship between
two transformed images.

Spatial Transformer Networks (STNs). is a method to integrate the learned
image deformation into a CNN architecture. A ST layer consists of a sub-
network to predict parameters with a deformation function. In [17], Inverse
Compositional STNs (IC-STNs) was proposed as an alignment framework by
adapting the standard Lucas-Kanade (LK) algorithm into STNs. In this model,
the image coordinate $x\in R^{2}$ and the input intensity at the same
coordinate $U_{(x)}\in R^{C}$ ($C$ denotes the number of channels) are used to
estimate the transformation parameter $\theta$ as $\mathcal{T}_{\theta}(x)$
and the transformed output $V_{(\tilde{x})}$ at the target coordinate
$\tilde{x}$.

$V_{(\tilde{x})}=\sum_{x}U(\mathcal{T}_{\theta}(x))K(\tilde{x},\mathcal{T}_{\theta}(x)),$
(1)

in which $K(.,.)$ is the convolution’s kernel to reflect the impact of this
process on each pixel. In this model, the kernel influences the whole image
and creates gradient values for all the pixels. This process involves back-
propagation throughout the pixels in the input image and to all the pixels in
the transformed one, which is expensive. The $K$ is adjusted in bilinear
interpolation, thus, $K(x,y)=0$ where $x$ and $y$ are not neighbours. For this
reason, the gradient just flows over sub-pixels gradient which reduce the
model performance under down-sampling. Moreover, the gradient of the sub-
pixels cannot address the sweeping changes that raised upon the changes of
transformation parameter $\theta$, because it is impossible to capture minor
transformations by direct neighbouring pixels. Therefore, the intensity values
with non-zero gradients from $K$ need to be selected in such a manner that is
invariant to deformations, in order to accommodate how the image would be
change under deformations. Our ESTN has been designed to address different
variations of geometry transformations by interpolating on the data manifold
distribution.

##  3 Entropy Spatial Transformer Network

In this section, we illustrate our proposed ESTN. We first discuss the
manifold learning, before introducing our model.

###  3.1 Interpolation on Data Manifold

In our model, we use the VAE, which is a strong generative model that contains
two modules: (1) an encoder, which is a recognition network to learn the
latent manifold representation of the data in the input space; (2) the
decoder, which learns how to reconstruct the data from the latent space. A
prior distribution is defined for the latent manifold representations, and a
mapping function is used to generate a surface in input space. In general,
each $x$ has an underlying manifold with a lower dimension represented as
${\cal{M}}\cong R^{C}$. We approximate Jacobian propagation to extend the
input data perturbed in the Jacobian order. In terms of regression, this
strategy can be implemented by generating auxiliary samples $z$ over the
underlying manifold $\cal{M}$. Therefore, due to the smoothness of the data on
the manifold and the smooth inductive bias of neural networks, realistic
results can be obtained by perturbing each input pixel with a Gaussian random
vector and marginalizing the samples $z$ in the context of linear regression.
In this work, for each pixel, we generate auxiliary noise $R_{z}$ within
$\cal{M}$ and the key point is how we ingeniously estimate $\cal{M}$ and
design $R_{z}$ to compel the intended smoothness. Therefore, linear regression
is used to find $w$ for minimizing $\|z-xw\|^{2}\in R^{(\hat{c})}$ ($xw$ is an
inner product so $C=\hat{c}$ and is a 2D matrix). To retain the input
dimension, we use Bernoulli probability $(p)$, therefore, the input $x$
represented as $(r\circ x)$ in which $r\in\\{0,1\\}$ is a random matrix. To
marginalize $\mu$ the samples, we formulate the following function

$\min_{r}\mu_{r\sim(p)}[\|z-(r\circ x)\|^{2}].$ (2)

where it returns a set of possible alignments compatible with $x$ drawn from
the manifold distribution and this allows to efficiently obtain both local and
contextual structures. To form $R_{\cal M}$ on the basis of the Jacobian, we
use

$R_{\cal{M}}=\mu_{x\sim dis(x)}\sum_{v\in T_{x}\cal{M}}\|J_{x}.v\|,$ (3)

where, $dis(x)$ is the true data distribution, $T_{x}\cal{M}$ denotes the
tangent space of $\cal{M}$ at $x$, $J$ is the Jacobin. These parameters are
used to regularize the manifold’s smoothness in view of the norm of its
Jacobian and the perception of the tangent space. This specifies a space
$m_{x}(s_{1}),...,m_{x}(s_{c})\in T_{x}\cal{M}$, in which $s_{i}$ denotes a
standard space of $R^{C}$. The linear approximation around $x$ and close to
$T_{x}\cal{M}$ creates a broader context to the local transformation, and
enhances the robustness of the network. Therefore, to compute the
transformation parameter $\mathcal{T}_{\theta}$ of a point $x_{i}$, the local
linear interpolation $L_{(x)}$ on the tangent subspace is calculated as

$L_{(x)}=\\{U(\mathcal{T}_{\theta}({x_{i}}))+l_{i}(\mathcal{T}_{\theta}(x)-\mathcal{T}_{\theta}(x_{i}))\\}.$
(4)

where $l_{i}$ is a matrix that stands for the local linearization over the
underlying data manifold that we look for. To accurately perform
interpolation, we should define a local chart around $x$ which is compatible
with the nature of the data manifold. Indeed, ${\mathcal{T}}_{\theta}$ needs
to be very sensitive to changes, for example, from point $x_{i}$ to another
transformed point $x_{j}$, but insensitive to any other changes. To obtain
$l_{i}$, the local linearization should be performed at the transformed
coordinate ${\mathcal{T}}_{\theta}({\cal M}_{x_{i}})$, while
${\mathcal{T}}_{\theta}({\cal M}_{x})$ is invariant. Thus, $l_{i}$ in (4)
represents the gradient of $L_{(x)}$ in terms of $x$ on the tangent manifold
and can be expressed by an orthogonal basis. To acquire $l_{i}$, firstly, we
perform interpolation near the point $x_{i}$ on the underlying data manifold
$T_{x_{i}}\cal M$, $\\{x_{i}^{k}\sim N(x_{i},\sigma)\\}$ where
$\\{k\in(1,2,…,K)\\}$, in which $N(x_{i},\sigma)$ is the data distribution
derived from $x_{i}$ with standard deviation $\sigma$. With the weighted
least-squares method [16], we can produce sample points in the $N\times N$
neighborhood of $x_{i}$ and estimate the nearest points. In the weighted
least-squares, we set a weight to each point $x$ in the $N\times N$
neighborhood of $x_{i}$ in accordance with the distance between $x$ and
$x_{i}$. The weight on $T_{x_{i}}\cal M$ is computed by

Figure 2: Interpolation methods. (a) Sample generation by linear
interpolation. (b) Sample generation based on data manifold.

$W_{i}x=\frac{1+\frac{N-1}{2}\sqrt{\frac{N-1}{2}}-\sqrt{(x-x_{i})^{2}}}{\frac{N-1}{2}\sqrt{\frac{N-1}{2}}}$
(5)

We streamline the notation $L({\mathcal{T}}_{\theta}({\cal M}_{x_{i}^{k}}))$
as $L_{i}^{k}$, and ${\mathcal{T}}_{\theta}({\cal M}_{x_{i}^{k}})$ as
$\tilde{x}_{i}^{k}$, where
$\tilde{x}_{i}^{k}=[\tilde{u}_{i}^{k},\tilde{v}_{i}^{k}]$, is based on the LK
[2], therefore we create two data matrices for the sampled points as follows

$y_{i}=x_{i}l_{i}|l_{i}\in R^{C},$ (6)
$y_{i}=\begin{bmatrix}L_{i}^{1}-L_{i}^{0}&L_{i}^{2}-L_{i}^{0}&...&L_{i}^{k-1}-L_{i}^{0}\end{bmatrix}^{ER}$
(7)
$x_{i}=\begin{bmatrix}\tilde{u}_{i}^{1}-\tilde{u}_{i}^{0}&\tilde{u}_{i}^{2}-\tilde{u}_{i}^{0}&&\tilde{u}_{i}^{k-1}-\tilde{u}_{i}^{0}\\\
\tilde{v}_{i}^{1}-\tilde{v}_{i}^{0}&\tilde{v}_{i}^{2}-\tilde{v}_{i}^{0}&...&\tilde{v}_{i}^{k-1}-\tilde{v}_{i}^{0}\\\
1&1&&1\end{bmatrix}^{ER}$ (8)

To estimate the displacement in the pixels, we analyse the changes in pixels
by using the intensity gradients of the image in all the pixels. Therefore,
$l_{i}$ can be computed by weighted least-squares $W_{i}$ on data manifold
$T_{x_{i}}\cal M$ as

$l_{i}=(x_{i}^{ER}W_{i}x_{i}+\xi A)^{-1}x_{i}^{ER}W_{i}y_{i}$ (9)

where $A$ denotes a $s\times s$ matching matrix, and $\xi$ is a distribution
with mean vector 0. To solve the nonconvex problem of Eq. (9), generally the
Tikhonov method is use as a regularizer [12]. However, the performance of this
method is sensitive to the choice of parameter $\gamma$. In our model,
projected gradient descent is used, whereas both the gradient and the
projections relied on the basis of the Tikhonov metric. By taking
$\frac{1}{\gamma}$ as the learning rate, the projected gradient is equivalent
to solve the non-convex problem with the entropy regularizer (ER) [1]:

$V(\tilde{x})=\sum_{u}U(\gamma{\mathcal{T}}_{\theta}(x))K(\tilde{x},{\mathcal{T}}_{\theta}(x))+\gamma
H({\mathcal{T}}_{\theta})$ (10)

in which
$H({\mathcal{T}}_{\theta})=\sum_{i,j}{\mathcal{T}}_{ij}\log{\mathcal{T}}_{ij}$,
here ${\mathcal{T}}_{ij}$ is an element of ${\mathcal{T}}_{\theta}$ and
represents the probability that sample point $i$ matches $j$. By replacing the
Tikhonov estimation with the ER, we derive the weight of the regularizer which
is a significant hyperparameter to improve the convergence. Fig. 2, shows the
performance of our proposed interpolation model in comparison with linear
interpolation. In our model, the generated samples are in line with the data
samples (context) and the data manifold.

###  3.2 Network Regularization

In this section, we define the proposed transformation layer and the essential
properties. To execute a transformation, the sampler needs to take a group of
sampling points along with $U_{(x)}$ to generate the output feature map
$V_{(x)}$. To perform such transformation, we project the proposed ESTN layer
kernel into the group of norm preserving kernels by defining its
differentiable values. We will discuss how to set the convolution layer to
maintain the norm via singular value regularization without dependency on the
decomposition of the values.

More precisely, in a convolution layer, let $k$, $c$, and $d$ be the kernel
size, input channels, and output channels respectively, therefore the gradient
is formulated as

$\nabla_{u}=\hat{W}\nabla_{v},$ (11)

in which $u\in R^{c}$ is a c-dimensional vector, $\nabla_{u}$ denotes the
input’s gradient, $v\in R^{d}$ is the $d$-dimensional vector, and $\nabla_{v}$
denotes the output’s gradient of the convolution. $\hat{W}$ is an $c\times
d$-dimensional matrix which performs backpropagation for the convolution
layer. This transformation illustrated as

$\begin{array}[]{lr}\nabla_{u}=\sum_{i=1}^{c}\psi_{i}p_{i}<\nabla_{v},q_{i}>,\\\
\nabla_{v}=\sum_{i=1}^{c}q_{i}<\nabla_{v},q_{i}>,\end{array}$ (12)

in which $\psi_{i}$ is a singular value of $\hat{W}$ and $p_{i}$, $q_{i}$ are
corresponding left and right singular vectors, respectively. Therefore, the
estimated values of the norms of gradients are calculated as

$\begin{array}[]{lr}E[\|\nabla_{u}\|_{2}^{2}]=\sum_{i=1}^{c}\psi_{i}^{2}E[c({\cal
M})|<\nabla_{v},q_{i}>|^{2}],\\\
E[\|\nabla_{v}\|_{2}^{2}]=\sum_{i=1}^{d}E[d({\cal
M})|<\nabla_{v},q_{i}>|^{2}],\\\ \end{array}$ (13)

in which, for the equal data range, we have: $p_{i}^{T}p_{j}=q_{i}^{T}q_{j}=1$
if $i=j$, otherwise 0. Therefore, to preserve the gradients’ norm, we propose
to have $E[\nabla_{u}]=E[\nabla_{v}]$, by assigning non-zero singular values
to $\psi$. It can be obtained by

$\psi^{2}=\frac{\sum_{i=1}^{d}E[|<\nabla_{v},q_{i}>|^{2}]}{\sum_{i,\psi_{i}}E[|<\nabla_{v},q_{i}>|^{2}]}$
(14)

Figure 3: Marginal Log-likelihood and qualitative results on MNIST, and
Fashion MNIST. The images stand for the context points (pixels quantity)
(first row), ESTN target predictions (second row), CVAE target predictions
(third row), and $\beta$-VAE target predictions. Each column represents a
specific percentage of the context distribution.

in which the denominator aggregation in a singular vector $q_{i}$ are based on
the non-zero singular values, $\psi_{i}\neq 0$. The ratio in Eq. (14) is the
ratio of the estimated gradient of $\nabla_{v}$ divided by the total gradient
that never fall in the null-space, or kernel of the matrix $\hat{W}$, so we
can approximate this assumption as $\frac{d}{\min(d,c)}$ . Based on this
assumption, $\frac{\min(d,c)}{d}$ of the gradient of $\nabla_{v}$ will fall in
the $\min(d,c)$-dimensional subspace, for basis $\\{q_{i}|\psi_{i}\neq 0\\}$.
Therefore, we can regularize the norm by changing the singular values ($\neq
0$) to $\sqrt{\frac{d}{\min(d,c)}}$. However, implementing this values
$\frac{d}{min(d,c)}$ is computationally expensive, due to the square root of
matrix and if not being well implemented, it can delay the training process.
To mitigate this issue, we use an iterative algorithm to compute the square
root of matrix with the help of matrix multiplication. Inspired by [10], this
iterative process is presented as:

$\begin{array}[]{lr}Y_{k+1}=\frac{1}{2}Y_{k}(3I-Z_{k}Y_{k}),\\\
Z_{k+1}=\frac{1}{2}(3I-Z_{k}Y_{k})Z_{k},\end{array}$ (15)

where $k$ is the iteration number. It has been shown that, $Y_{k}$ and $Z_{k}$
converge to the $A^{\frac{1}{2}}$ and $A^{-\frac{1}{2}}$ respectively.
Therefore, as the iterations only entail matrix multiplication, we can
efficiently implement them on GPUs.

Table 1: Marginal likelihood evaluation on five datasets. Model | Param | MNIST | F-MNIST | SVHN | CelebA | ZSMM
---|---|---|---|---|---|---
VAE | 210K |  1.08 $\pm$ 0.04 |  2.29 $\pm$ 0.02 |  3.76 $\pm$ 0.03 |  2.92 $\pm$ 0.04 |  0.48 $\pm$ 0.02
$\beta$-VAE [3] | 415K |  1.14 $\pm$ 0.03 |  2.36 $\pm$ 0.01 |  3.88 $\pm$ 0.02 |  3.08 $\pm$ 0.02 |  0.61 $\pm$ 0.01
CVAE [22] | 433K |  1.23 $\pm$ 0.02 |  2.41 $\pm$ 0.02 |  3.91 $\pm$ 0.04 |  3.13 $\pm$ 0.02 |  0.74 $\pm$ 0.03
ESTN (ours) | 219K | 1.28 $\pm$ 0.01 | 2.47 $\pm$ 0.00 | 3.95 $\pm$ 0.01 | 3.22 $\pm$ 0.01 | 1.12 $\pm$ 0.01
Figure 4: Performance evaluation. From each dataset a random image is selected
and different quantity of noise is generated. The first rows from left to
right contain the noisy images (5%, 10%, 20%, 30%, 50% of random pixels). The
second rows show the performance of ESTN in image completion. Table 2:
Classification errors (%) of different models with various numbers of data on
Fashion MNIST, SVHN, and CelebA.

Method | 1k labels | 3k labels | 5k labels
---|---|---|---
CNN |  14.47 $\pm$ 0.23 |  12.63 $\pm$ 0.24 |  9.52 $\pm$ 0.15
VAE |  14.21 $\pm$ 0.30 |  12.49$\pm$ 0.18 |  9.46 $\pm$ 0.11
STN [11] |  11.84$\pm$ 0.21 |  9.37$\pm$ 0.05 |  6.76 $\pm$ 0.06
GCNN [26] |  9.63 $\pm$ 0.10 |  7.08 $\pm$ 0.19 |  4.49 $\pm$ 0.15
L-STN [12] |  8.92 $\pm$ 0.12 |  6.83 $\pm$ 1.24 |  4.24 $\pm$ 1.08
SE3Mov [23] |  8.07 $\pm$ 0.24 |  4.87 $\pm$ 0.27 |  2.83 $\pm$ 0.15
ScDCFNet [34] |  7.98 $\pm$ 0.09 |  4.72 $\pm$ 0.12 |  2.65 $\pm$ 0.09
ESTN (ours) | 7.21$\pm$ 0.06 | 4.43 $\pm$ 0.8 | 2.19 $\pm$ 0.16

Method | SVHN 10k labels | CelebA 20k labels
---|---|---
CNN |  16.45 $\pm$ 0.17 |  22.24 $\pm$ 0.12
VAE |  16.81 $\pm$ 0.24 |  21.67 $\pm$ 0.35
STN [11] |  14.97 $\pm$ 0.16 |  18.82 $\pm$ 0.22
GCNN [26] |  12.85 $\pm$ 0.26 |  16.98 $\pm$ 0.35
L-STN [12] |  11.45 $\pm$ 0.14 |  15.31 $\pm$ 0.21
SE3Mov [23] |  9.85 $\pm$ 0.27 |  13.67 $\pm$ 0.24
ScDCFNet [34] |  9.56 $\pm$ 0.09 |  13.63 $\pm$ 0.14
ESTN (ours) | 9.14 $\pm$ 0.24 | 13.24 $\pm$ 0.16

Table 3: The architecture of CNNs for experiments on Fashion MNIST, SVHN and CelebA. All Conv. and fully connected (FC) layers are followed by ReLU activations and batch normalization. F-MNIST network | SVHN network | CelebA network
---|---|---
$32\times 32$ image
$3\times 3$ conv. 32 |  $3\times 3$ conv. 64 |  $3\times 3$ conv. 128
$3\times 3$ conv. 32 |  $3\times 3$ conv. 64 |  $3\times 3$ conv. 128
EST layer
$3\times 3$ conv. 32 |  $3\times 3$ conv. 64 |  $3\times 3$ conv. 128
$2\times 2$ max-pool, stride 2 dropout
$3\times 3$ conv. 64 |  $3\times 3$ conv. 128 |  $3\times 3$ conv. 256
$3\times 3$ conv. 64 |  $3\times 3$ conv. 128 |  $3\times 3$ conv. 256
EST layer
$3\times 3$ conv. 64 |  $3\times 3$ conv. 128 |  $3\times 3$ conv. 256
$2\times 2$ max-pool, stride 2 dropout
$3\times 3$ conv. 64 |  $3\times 3$ conv.128 |  $3\times 3$ conv. 512
$1\times 1$ conv. 64 |  $1\times 1$ conv. 128 |  $1\times 1$ conv. 256
EST layer
$1\times 1$ conv. 64 |  $1\times 1$ conv. 128 |  $1\times 1$ conv. 128
global average pool, $6\times 6\rightarrow 1\times 1$
FC 64$\rightarrow$10 |  FC 128$\rightarrow$ 10 |  FC 128 $\rightarrow$ 10

##  4 Experiments

In this section, the performance of our ESTN model is evaluated on image
reconstruction and classification, and compared with several state-of-the-art
methods. Image reconstruction can be seen as a prediction of pixel intensities
to determine a target pixel’s location which is sensitive to the detected
contextual points (pixels). In our experiments, the contextual set may be
diverse but the target set consists of all the image pixels. For image
reconstruction, the same model is utilized for 1-D function regression
(excluding the setting of the last layer as 3-D information for RGB images).
The models are evaluated on five public datasets: MNIST [15], Fashion-MNIST
[30], SVHN [20], zero shot multi-MNIST (ZSMM) [19], and CelebA [18]. In all
our experiments, we use a VAE architecture based on [25]. The encoder contains
4 layers (convolution sizes = $[32,32,64,2\times 8]$, strides = [1, 2, 2, 1],
and kernel sizes = [5, 4, 3, 5]). ReLU activation function is used in the
first three layers, whereas the final layer is linear. The decoder also
contains 4 layers (convolution sizes = [32, 16, 16, 1], strides = [1, 2, 2,
1], and kernel sizes = [4, 5, 5, 3]). ReLU activation function is used in the
first three layers, whereas the sigmoid activation function is used in the
final layer, and batch normalization is used on all the layers except the last
ones. To train our network, Adam is used to optimize the weights, with a
learning rate of 0.0001, weight decay = 0.001, and momentum = 0.9. Our ESTN
contains two ST layers that are trained together with the encoder and the
decoder, and the loss is determined as the difference between the decoder’s
output and the target. We implement the proposed ESTN in Pytorch, and run it
on a Tesla P40 GPU.

Predicting pixel values can be viewed as a regression problem that maps a 2D
pixel position to its pixel intensity. In general, each image represents a
single realization of the process that was sampled on a fixed two-dimensional
grid. We use the regular train/test split with up to 250 contextual points
(image pixels) to train the ESTN on the datasets. In Table 1, we report the
mean and standard deviations of log-likelihood of ESTN compare to different
VAE-based models such as Conditional VAE (CVAE) and $\beta$-VAE. CVAE [22] is
a modification of VAE to reconstruct an image conditioned on the given
attributes and $\beta$-VAE [3] is another variation of the VAE that is often
used as an unsupervised approach for the representation learning. Results
indicate that ESTN performs considerably better than CVAE and $\beta$-VAE.
Fig. 3 shows a visual representation of the marginal likelihoods of the images
with a random number of contextual pixels, along with the distribution results
of the ESTN, CVAE and $\beta$-VAE. Results show that ESTN has accurate
contextual reconstruction performance compared to the other baselines.
Moreover, Fig. 4 shows the performance of ESTN on qualitative samples with
different noise set-ups. These results show that even when the target pixels
are filled, ESTN generates smooth results due to the effective incorporation
of transformer layers and manifold interpolation.

Figure 5: Model performance, efficiency and inference latency comparison on
Fashion MNIST. Performance is measured with batch size 1 on the same machine.

###  4.1 Fashion MNIST

For the classification task, we applied ESTN to learn the underlying data
manifold. To undertake experiments on the Fashion MNIST, we use three sets of
the labeled data, i.e., 1k, 3k and 5k, respectively and preserve 150 data
points for validation from the initial training set. The network architecture
for classification is reported in Table 3. We set the batch size to 32. The
model is trained for 5k epochs. Adam optimizer is adopted with the initial
learning rate of 0.001, and linearly decayed on every 1k epochs. The training
time was about 20 hours. The hyperparameters are tuned in the scale of the
tangent perturbation where we use rotations, and scale changes with respect to
the vertical-horizontal directions. We adopt Gaussian noise with standard
deviation 1. As pre-processing, we normalized the pixels value into [-1, 1],
standard deviation is set to $\pi/4$ for rotations, $\sqrt{2}$ for scale
changes either in vertical or horizontal direction, and 0.2 for translations.
0.1 applied for all the other hyperparameters. In the VAEs network, for
estimating the underlying manifold, six convolution layers and two fully
connected layers are used, where the decoder is symmetric to the encoder,
replacing the convolution layers with the deconvolution ones.

Table 2 reports the results. ESTN has lower classification errors as compared
to STN and L-STN, in various numbers of the labeled data. The performance of
ESTN is relying on the prediction of the underlying data manifold. Therefore,
ESTN contributes better results compared to other learning-based equivariant
CNNs.

###  4.2 SVHN and CelebA

We performed two sets of experiments to prove the success of ESTN in
classification, SVHN with 10k labeled data, and CelebA with 20k labeled data.
As pre-processing, we normalized the pixels value into [-1, 1]. The networks’
structure is provided in Table 3. We set the batch size to 32. The SVHN’s
network is trained for 8k epochs, and the CelebA’s network is trained for 12k
epochs. The initial learning rate is set to 0.001, and linearly decayed on
every 2k epochs. The other hyperparameters are similar to the Fashion MNIST’s
setting. Table 2 shows the results of our experiments on SVHN and CelebA. The
evaluation results indicate the superior performance of ESTN compared to the
other state-of-the-art approaches. Particularly, in comparison with L-STN that
used linear interpolation, ESTN shows significant improvements with very
similar network settings. In addition, ESTN has higher accuracy compared to
the recent developed methods such as ScDCFNet [34] and SE3Mov [23].

Table 4: Ablation experiments on three datasets. Model | MNIST | F-MNIST | SVHN
---|---|---|---
ESTN (ours) | 1.28 $\pm$ 0.01 | 2.47 $\pm$ 0.00 | 3.95 $\pm$ 0.01
ESTN w/o norm. |  1.20 $\pm$ 0.03 |  2.38 $\pm$ 0.02 |  3.84 $\pm$ 0.02
ESTN w/o ER |  1.22 $\pm$ 0.02 |  2.41 $\pm$ 0.01 |  3.89 $\pm$ 0.00

###  4.3 Ablation study

We performed the ablation study to show the impact of each module on the
overall performance of ESTN. The evaluation results are shown in Tables 4. In
our proposed model, both the entropy regularizer, and the norm-regularization
methods are extremely helpful to improve the performance, otherwise the
network fails to impose the manifold invariance or cannot be robust against
the off-manifold noise. Our experiments illustrate that sampling based on the
underlying data manifold can significantly improves the performance of ESTN.
ESTN performs satisfactorily with a wide range of manifold coordinate charts.
In our experiments, CNN produces more detailed images than VAEs, on the other
hand, VAEs actively cooperates with ESTN to produce more realistic samples. It
has been observed that sampling based on tangent bundle data manifold helps
ESTN to explore high-dimensional data by characterizing the learning task.

Fig. 5 (a) shows, ESTN has much more rapid decreases in reconstruction error
and lower values at convergence. Compare to the other baselines ESTN has
around $10^{-2}$ lower reconstruction errors. In Fig. 5 (b), we evaluate the
efficiency of ESTN in comparison with that of the other models. The results
indicate that ESTN can considerably enhance the system efficiency. For
instance, ESTN surpasses CVAE and $\beta$-VAE . This suggests that we can have
better performance by interpolating based on the underlying data manifold, and
capture the attention to the significance of norm-regularization by modifying
a small part of the network. Fig. 5 (c) demonstrate the GPU latencey-accuracy
curves for the models, where ESTN achieves higher accuracy in shorter time
with much fewer FLOPS against the other models. Our ESTN is computational
cheaper, for example, up to 2.1$\times$ faster than CVAE.

##  5 Conclusion

In this paper, we have proposed a new variant of spatial transformer networks
based on sampling on the tangent space to explore more samples in the search
space along the data manifold to improve the gradients and preserve the norm
for differentiable image handling. We have demonstrated that ESTN have been
successful on a wide range of challenges and is highly scalable. In terms of
the arithmetic operations in the vector space, our tangent space manifold
convolution is sensible. In fact, the proposed ESTN can be used as an example
to improve the gradient estimation with better performance on image
reconstruction and classification tasks. Moreover, since the proposed model
preserves the norm of gradients, it can be used in highly deep networks and
reduces the optimization bottlenecks of deep networks.

## References

  * [1] Jean-David Benamou, Guillaume Carlier, Marco Cuturi, Luca Nenna, and Gabriel Peyré.  Iterative bregman projections for regularized transportation problems.  SIAM J. on Sci. Comput., 37(2):A1111–A1138, 2015.
  * [2] Che-Han Chang, Chun-Nan Chou, and Edward Y Chang.  Clkn: Cascaded lucas-kanade networks for image alignment.  In CVPR, pages 2213–2221, 2017.
  * [3] Tian Qi Chen, Xuechen Li, Roger B Grosse, and David K Duvenaud.  Isolating sources of disentanglement in variational autoencoders.  In Advances in Neural Information Processing Systems (NIPS), pages 2610–2620, 2018.
  * [4] Taco Cohen and Max Welling.  Group equivariant convolutional networks.  In ICML, pages 2990–2999, 2016.
  * [5] Taco S Cohen, Mario Geiger, Jonas Köhler, and Max Welling.  Spherical cnns.  ICLR, 2018.
  * [6] Taco S Cohen, Maurice Weiler, Berkay Kicanaoglu, and Max Welling.  Gauge equivariant convolutional networks and the icosahedral cnn.  ICML, 2019.
  * [7] Jan Gerken, Oscar Carlsson, Hampus Linander, Fredrik Ohlsson, Christoffer Petersson, and Daniel Persson.  Equivariance versus augmentation for spherical images.  In ICML, pages 7404–7421, 2022.
  * [8] Lingshen He, Yuxuan Chen, Yiming Dong, Yisen Wang, Zhouchen Lin, et al.  Efficient equivariant network.  NIPS, 34:5290–5302, 2021.
  * [9] Joao F Henriques and Andrea Vedaldi.  Warped convolutions: Efficient invariance to spatial transformations.  In ICML, pages 1461–1469, 2017.
  * [10] Nicholas J Higham.  Stable iterations for the matrix square root.  Numerical Algorithms, 15(2):227–242, 1997.
  * [11] Max Jaderberg, Karen Simonyan, Andrew Zisserman, et al.  Spatial transformer networks.  In NIPS, pages 2017–2025, 2015.
  * [12] Wei Jiang, Weiwei Sun, Andrea Tagliasacchi, Eduard Trulls, and Kwang Moo Yi.  Linearized multi-sampling for differentiable image transformation.  In ICCV, pages 2988–2997, 2019.
  * [13] Jinwoo Kim, Saeyoon Oh, Sungjun Cho, and Seunghoon Hong.  Equivariant hypergraph neural networks.  In ECCV, pages 86–103, 2022.
  * [14] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114, 2013.
  * [15] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [16] Erbo Li, Hanlin Mo, Dong Xu, and Hua Li.  Image projective invariants.  IEEE TPAMI, 41(5):1144–1157, 2018.
  * [17] Chen-Hsuan Lin and Simon Lucey.  Inverse compositional spatial transformer networks.  In CVPR, pages 2568–2576, 2017.
  * [18] Ziwei Liu, Ping Luo, Xiaogang Wang, and Xiaoou Tang.  Large-scale celebfaces attributes (celeba) dataset.  Retrieved August, 15:2018, 2018.
  * [19] Jiang Lu, Jin Li, Ziang Yan, and Changshui Zhang.  Zero-shot learning by generating pseudo feature representations.  arXiv preprint arXiv:1703.06389, 2017.
  * [20] Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Bo Wu, and Andrew Y Ng.  Reading digits in natural images with unsupervised feature learning.  2011\.
  * [21] Hyeonwoo Noh, Andre Araujo, Jack Sim, Tobias Weyand, and Bohyung Han.  Large-scale image retrieval with attentive deep local features.  In ICCV, pages 3456–3465, 2017.
  * [22] Siddharth Ramchandran, Gleb Tikhonov, Otto Lönnroth, Pekka Tiikkainen, and Harri Lähdesmäki.  Learning conditional variational autoencoders with missing covariates.  arXiv preprint arXiv:2203.01218, 2022.
  * [23] Mateus Sangalli, Samy Blusseau, Santiago Velasco-Forero, and Jesus Angulo.  Moving frame net: Se (3)-equivariant network for volumes.  ICML, 2022.
  * [24] Vıctor Garcia Satorras, Emiel Hoogeboom, and Max Welling.  E (n) equivariant graph neural networks.  In ICML, pages 9323–9332. PMLR, 2021.
  * [25] Lukas Schott, Jonas Rauber, Matthias Bethge, and Wieland Brendel.  Towards the first adversarially robust neural network model on mnist.  In Proceedings of the International Conference on Learning Representations (ICLR), 2018.
  * [26] Jaspreet Singh and Chandan Singh.  Learning invariant representations for equivariant neural networks using orthogonal moments.  In IJCNN, pages 1–8, 2022.
  * [27] Giovanni Sutanto, Isabel M Rayas Fernández, Peter Englert, Ragesh K Ramachandran, and Gaurav S Sukhatme.  Learning equality constraints for motion planning on manifolds.  CoRL, 2020.
  * [28] Rui Wang, Robin Walters, and Rose Yu.  Approximately equivariant networks for imperfectly symmetric dynamics.  ICML, 2022.
  * [29] Daniel E Worrall, Stephan J Garbin, Daniyar Turmukhambetov, and Gabriel J Brostow.  Harmonic networks: Deep translation and rotation equivariance.  In CVPR, pages 5028–5037, 2017.
  * [30] Han Xiao, Kashif Rasul, and Roland Vollgraf.  Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms.  arXiv preprint arXiv:1708.07747, 2017.
  * [31] Jian Yang, Fuxin Li, and Jue Wang.  A better scaled local tangent space alignment algorithm.  In IJCNN, volume 2, pages 1006–1011, 2005.
  * [32] Shifeng Zhang, Longyin Wen, Xiao Bian, Zhen Lei, and Stan Z Li.  Single-shot refinement neural network for object detection.  In CVPR, pages 4203–4212, 2018.
  * [33] Zhen Zhong, Guobao Xiao, Linxin Zheng, Yan Lu, and Jiayi Ma.  T-net: Effective permutation-equivariant network for two-view correspondence learning.  In ICCV, pages 1950–1959, 2021.
  * [34] Wei Zhu, Qiang Qiu, Robert Calderbank, Guillermo Sapiro, and Xiuyuan Cheng.  Scaling-translation-equivariant networks with decomposed convolutional filters.  JMLR, 23(68):1–45, 2022.