A Spectral Regularizer for Unsupervised Disentanglement

Aditya Rameshopenai,formerNYU
Youngduck Choiyale
Yann LeCunnyu

formerNYUContinuation of work done at New York University
yaleDepartment of Computer Science, Yale University
nyuDepartment of Computer Science, New York University

Aditya<EMAIL_ADDRESS>Machine Learning, ICML, Disentangled Representations, Generative Modeling

A generative model with a disentangled representation allows for control over independent aspects of
the output. Learning disentangled representations has been a recent topic of great interest, but
it remains poorly understood. We show that even for GANs that do not possess disentangled
representations, one can find curved trajectories in latent space over which local disentanglement
occurs. These trajectories are found by iteratively following the leading right-singular vectors of
the Jacobian of the generator with respect to its input. Based on this insight, we describe an
efficient regularizer that aligns these vectors with the coordinate axes, and show that it can be
used to induce disentangled representations in GANs, in a completely unsupervised manner.

§ INTRODUCTION

Rapid progress has been made in the development of generative models capable of producing realistic
samples from distributions over complex, high-resolution natural images (Karras et al., 2017,
Brock et al., 2018). Despite the success of these models, it is unclear how one can achieve control over
the data-generation process when labeled instances are not available. Perturbing an
individual component of a latent variable typically results in an unpredictable change to the
output. When the model has a disentangled representation, these changes become interpretable, and
each component of the latent variable affects a distinct attribute. So far, the problem of learning
disentangled reprersentations remains poorly understood. Most approaches for doing this have focused
on VAEs (Kingma & Welling, 2013, Rezende et al., 2014), which produce blurry samples in practice. The few
approaches that have been developed for GANs, such as InfoGAN [Chen et al., 2016], have had comparably
limited success.

Trajectory obtained by following the leading eigenvector of $M_z(z_0)$ from a fixed
embedding $z_0 \in \reals^{256}$ in the latent space of a GAN trained on the LSUN Bedrooms dataset
at ${256 \times 256}$ resolution. The camera angle is varied smoothly, while the content of the
bedroom is preserved. The parameters used for Algorithm <ref> are given
by ${\alpha, \rho, N} = 1.5 \cdot 10^{-2}, 0.99, 2000$; we show iterates $0, 20, \ldots, 180$.  The
architecture and training procedure are as described in [Mescheder et al., 2018]. Videos of the
trajectories for the first four eigenvectors can be viewed at this URL:

Learning disentangled representations is of interest, because they allow us gauge the ability of
generative models to benefit downstream tasks. For example, a robotics application may require a
vision model that reliably determines the orientation of an object, subject to changes in viewpoint,
lighting conditions, and intra-class variation in appearance. A generative model with independent
control over the full extent of each of these attributes may help to achieve the level of robustness
that is required for these applications. Moreover, the ability to synthesize concepts in ways that
are not present in the original data (e.g., rotation of an object to a position not encoutered
during training) is a useful benchmark for reasoning and out-of-distribution generalization. Often
times, the goal of training a generative model is to observe that this type of generalization

Following [Higgins et al., 2016], we say a generative model $G : Z \to X$, where ${Z \subset \reals^m}$
and ${X \subset \reals^n}$, possesses a disentangled representation when it satisfies two important
properties. The first property is independence of the components of the latent variable.
Roughly speaking, this means that a perturbation to any component of a latent variable should result
in a change to the output that is distinct, in some sense, from the changes resulting from
perturbations to the other components. For example, if the first component of the latent variable
controls hair length, then the other components should not influence this attribute. The second
property is interpretability of the changes in $X$ resulting from these perturbations.
Suppose that the process of sampling from the distribution $\pdata$ modeled by $G$ can be realized
by a simulator that is parameterized over a list of independent, scalar-valued factors or
attributes. These attributes might correspond to concepts such as lighting, azimuth, gender, age,
and so on. This property is met when $G(z + \alpha e_i)$ results in a change along exactly one of
these attributes, for each $i \in [1, m]$, where $e_i$ is the $i$th standard basis vector
and $\alpha > 0$. Quantifying the extent to which this property holds is generally challenging.

Classical dimensionality reduction techniques, such as PCA and metric MDS, can be used to obtain a
latent representation satisfying the first property. For instance, the former ensures that changes
to the output resulting from changes to distinct components of the latent variable are orthogonal,
which is a particularly restrictive form of distinctness. These techniques are guaranteed to
faithfully represent the underlying structure only when the data occupy a linear subspace. For many
applications of interest, such as image modeling, the data manifold could be a highly curved and
twisted surface that does not satisfy this assumption. The Swiss Roll dataset is a simple example
for which this is the case: two points with small euclidean distance could have large geodesic
distance on the data manifold. Since both techniques measure similarity between each pair of points
using euclidean distance, distant outputs in $X$ could get mapped to nearby embeddings in $Z$. Thus,
the second property often fails to hold in practice.

Deep generative models have the potential to learn representations that satisfy both of the required
properties for disentanglement. Our work focuses on the case in which $G$ is a GAN generator.
Let ${\jg(z_0) \in \reals^{n \times m}}$ be the Jacobian of $G$ evaluated at the latent
variable ${z_0 \in Z}$. Our main contributions are as follows:

* We show that the leading right-singular vectors of $\jg(z_0)$ can be used to obtain a local
disentangled representation about a neighborhood of $z_0$;

* We show that following the path determined by a leading right-singular vector in the latent
space of a GAN yields a trajectory along which local changes to the output are interpretable; and

* We formulate a regularizer that induces disentangled representations in GANs, by aligning the
top $k$ right-singular vectors ${v_1, \ldots, v_k}$ with the first $k$ coordinate directions, ${e_1,
\ldots, e_k}$.

§ RELATED WORK

To date, the two most successful approaches for unsupervised learning of disentangled
representations are the $\beta$-VAE [Higgins et al., 2016] and InfoGAN [Chen et al., 2016]. The former proposes
increasing the weight $\beta$ for the KL-divergence term in the objective of the VAE (Kingma & Welling, 2013,
Rezende et al., 2014), which is normally set to one. This weight controls the tradeoff between
minimizing reconstruction error, and minimizing the KL-divergence between the approximate posterior
and the prior of the decoder. The authors observe that as the KL-divergence is reduced, more
coordinate directions in the latent space of the decoder end up corresponding to disentangled
factors. To measure the extent to which this occurs, they describe a disentanglement metric score,
which we use in Section <ref>. Follow-up work [Burgess et al., 2018]
analyzes $\beta$-VAE from the perspective of information channel capacity, and describes principled
way of controlling the tradeoff between reconstruction error and disentanglement. Several variants
of the $\beta$-VAE have also since been developed (Kim & Mnih, 2018, Chen et al., 2018, Esmaeili et al., 2018).

InfoGAN Chen et al., 2016 augments the GAN objective [Goodfellow et al., 2014] with a term that maximizes the mutual
information between the generated samples and small subset of latent variables. This is done by
means of an auxiliary classifier that is trained along with the generator and the discriminator.
Adversarial training offers the potential for a high degree of realism that is difficult to achieve
with models like VAEs, which are based on reconstruction error. However, InfoGAN finds fewer
disentangled factors than VAE-based approaches on datasets such as CelebA [Liu et al., 2015] and
3DChairs [Aubry et al., 2014], and offers limited control over each latent factor (e.g., maximum
rotation angle for azimuth). The mutual information regularizer also detriments sample quality.
Hence, the development of an unsupervised method for learning disentangled representations with
GANs, while meeting or exceeding the quality of those found by VAE-based approaches, remains an open

Our work makes an important step in this direction. In contrast to previous approaches, which are
based on information theory, we leverage the spectral properties of the Jacobian of the generator.
This new perspective not only allows us to induce high-quality disentangled representations in GANs,
but also to find curved paths over which disentanglement occurs, for GANs that do not possess
disentangled representations.

§ IDENTIFYING LOCAL DISENTANGLED FACTORS

We begin by describing how the singular value decomposition of $\jg(z_0)$ can be used to define a
local generative model about $z_0$ that satisfies the first property of disentangled
representations. In particular, we show that the left-singular vectors of $\jg(z_0)$ form an
orthonormal set of directions from $G(z_0)$ along which the magnitudes of the instantaneous changes
in $G$ are maximized. Since perturbations from $z_0$ along the right-singular vectors result in
changes from $G(z_0)$ along the left-singular vectors, the right-singular vectors result in distinct
changes to the output of $G$. This simple relationship allows us to define a local generative model
that satisfies the independence property.

First, we show that the perturbations from $z_0$ along the right-singular vectors maximize the
magnitude of the instantaneous change in $G(z_0)$. Given an arbitrary vector ${v \in \reals^m}$, the
directional derivative

\begin{equation}
\lim_{\epsilon \to 0} \frac{G(z_0 + \epsilon v) - G(z_0)}{\epsilon} = \jg(z_0) v.
\end{equation}

measures the instantaneous change in $G$ resulting from a perturbation along $v$ from $z_0$. The
magnitude of this change is given by

\begin{align}
\| \jg(z_0) v \|
&= (v^t \jg(z_0)^t \jg(z_0) v)^{1 / 2} \notag\\
&\eqqcolon (v^t M_z(z_0) v)^{1 / 2}
\eqqcolon n_z(v, z_0),
\end{align}

which is a seminorm involving the positive semidefinite matrix ${M_z(z_0) \coloneqq \jg(z_0)^t
\jg(z_0) \in \reals^{m \times m}}$. The unit-norm perturbation from $z_0$ that maximizes $n_z(v,
z_0)$ is given by

\begin{equation}
\coloneqq \max_{v \in \mathbb{S}^{m - 1}} n_z(v, z_0)
= \max_{v \in \mathbb{S}^{m - 1}} v^t \jg(z_0)^t \jg(z_0) v,
\end{equation}

where $\mathbb{S}^{m - 1} \subset \reals^m$ is the unit sphere. This is the first eigenvector
of $M_z(z_0)$, which coincides with the first right-singular vector of $\jg(z_0)$. It follows from
the singular value decomposition of $\jg(z_0)$ that the first left-singular vector is given by ${u_1
= \sigma_1^{-1} \jg(z_0) v_1}$, where $\sigma_1$ is the first singular value. Hence, a perturbation
from $z_0$ along $v_1$ maximizes the magnitude of the instantaneous change in $G(z_0)$, and this
change occurs along $u_1$.

Next, we consider the unit-norm perturbation orthogonal to $v_1$ that maximizes $n_z(v, z_0)$. It is
given by

\begin{equation}
v_2 \coloneqq \max_{v \in \mathbb{S}^{d - 1} \;\cap\;
\vspan(v_1)^\perp} v^t \jg(z_0)^t \jg(z_0) v.
\end{equation}

This is the second eigenvector of $M_z(v, z_0)$, which coincides with the second right-singular
vector of $\jg(z_0)$. As before, we get ${u_2 = \sigma_2^{-1} \jg(z_0) v_2}$, where $\sigma_2$ is
the second singular value. So a perturbation from $z_0$ along $v_2$ results in an instantaneous
change in $G(z_0)$ along $u_2$. Continuing in this way, we consider the $k$th unit-norm perturbation
orthogonal to ${v_1, \ldots, v_{k - 1}}$ that maximizes $n_z(v, z_0)$, for each $k \in [2, r]$,
where $r \coloneqq \rank(M_z(z_0))$. This shows that the right-singular vectors of $\jg(z_0)$
maximize the magnitude of the instantaneous change in magnitude of $G(z_0)$, and these changes occur
along the corresponding left-singular vectors.

Now, we use the right-singular vectors of $\jg(z_0)$ to define a local generative model about $z_0$.
Consider the function ${\bar{G}_{z_0} : \reals^r \to \reals^n}$,

\[
\bar{G}_{z_0} : \alpha \mapsto G\left( z_0 + \sum_{i \in [1, r]} \alpha_i v_i \right).
\]

The components of $\alpha$ control perturbations along orthonormal directions, and these directions
also result in orthonormal changes to $G(z_0)$. Hence, $\bar{G}_{z_0}$ satisfies the first property
for a generative model to possess a disentangled representation, but only about a neighborhood
of $z_0$. Figure <ref> in Appendix <ref>
investigates whether $\bar{G}_{z_0}$ also satisfies the second property: interpretability of changes
to individual components of $\alpha$. We can see that perturbations along the leading
eigenvectors $M_z(z_0)$, especially the principal eigenvector, often result in the most drastic
changes. These changes are interpretable, and tend to make modifications to isolated attributes of
the face. To see this in more detail, we consider the top two rows of subfigure (g). Movement along
the first two eigenvectors changes hair length and facial orientation; movement along the third
eigenvector decreases the length of the bangs; movement along the fourth and fifth eigenvectors
changes background color; and movement along the sixth and seventh eigenvectors changes hair color.

§ FINDING QUASI-DISENTANGLED PATHS

Trajectories obtained by following the first three eigenvectors of $M_z(z_0)$ from a fixed
embedding $z_0 \in \reals^3$ in the latent spaces of two GANs with identical architectures, trained
on the dSprites dataset [Matthey et al., 2017]. Subfigure (a) shows plots of the
trajectories ${\{\gamma_k\}_{k \in [1, 3]}}$, and subfigure (b) shows the outputs of $G$ at iterates $0,
100, \ldots, 2000$ of each trajectory $\gamma_k$, for $k \in [1, 3]$, from top to bottom,
respectively. Subfigures (c) and (d) show the same information for an alignment-regularized
GAN (${k, \lambda = 3, 0.1}$); the trajectories are now axis-aligned, as expected. The parameters
used for Algorithm <ref> are given by ${\alpha, \rho, N} = 5 \cdot 10^{-3},
0.99, 2000$. Details regarding the model architecture are given in
Appendix <ref>.

Generative models known to possess disentangled representations, such
as $\beta$-VAEs [Higgins et al., 2016], allow for continuous manipulation of attributes via perturbations
to individual coordinates. Starting from a latent variable $z_0 \in Z$, we can move along the
path ${\gamma : t \mapsto z_0 + t e_i}$ in order to vary a single attribute of $G(z_0)$, while
keeping the others held fixed. GANs are not known to learn disentangled representations
automatically. Nonetheless, the previous section shows that the local generative
model $\bar{G}_{z_0}$ does possess a disentangled representation, but only about a neighborhood of
the base point $z_0$. We explore whether it is possible to extend this local model to obtain
disentanglement along a continuous path from $z_0$. To do this, we construct a trajectory ${\gamma_k
: \reals \to \reals^n}$, ${t \mapsto G(\gamma_k(t))}$ by repeatedly following the $k$th leading
eigenvector of $M_z(\gamma_k(t))$, where $\gamma(0) \coloneqq z_0$. The procedure used to do this is
given by Algorithm <ref>.

We first test the procedure on a toy example for which it is possible to explicitly plot the
trajectories. We use the dSprites dataset [Matthey et al., 2017], which consists of $64 \times 64$ white
shapes on black backgrounds. Each shape is completely determined by six attributes: symbol (square,
ellipse, or heart), scale (6 values), rotation angle from $0$ to $2 \pi$ (40 values), and $x$- and
$y$-positions (30 values each). We trained a GAN on this dataset using a latent variable size of
three, fewer than the six latent factors that determine each shape.
Figure <ref> shows that outputs of the generator along these trajectories
vary locally along only one or two attributes at a time. Along the first trajectory $\gamma_1$, the
generator first decreases the scale of the square, then morphs it into a heart, increases the scale
again, and finally begins to rotate it. Similar comments apply to the other two trajectories,
$\gamma_2$ and $\gamma_3$.

Next, we test the procedure on a GAN trained on the CelebA dataset [Liu et al., 2015] at $64 \times 64$
resolution. Figure <ref>(a) shows the trajectories $\gamma_1$ starting at
four fixed embeddings ${z_1, \ldots, z_4}$. Although the association between the ordinal $k$ of the
eigenvector $v_k$ and the attribute of the image being varied is not consistent throughout latent
space, local changes still tend to occur along only one or two attributes at a time.
Figure <ref>(b) shows the trajectories ${\gamma_1, \gamma_2, \gamma_3}$
and $\gamma_5$, all starting from the same fixed embedding $z_5$. As is the case for the dSprites
dataset, we can see that trajectories $\gamma_k$ for distinct $k$ tend to effect changes to $G(z_0)$
along distinct attributes. These results suggest that, along trajectories from $z_0$ determined by a
leading eigenvector of $M_z(z_0)$, changes in the output tend to occur along isolated attributes.

[t]Procedure to trace path determined by $k$th leading

[1]
${\mv : \reals^m \times \reals^m \to \reals^d}$, ${z, v \mapsto M_z(z) v}$ is a function
that computes matrix-vector products with the implicitly-defined matrix $M_z(z) \in \reals^{m \times
$z \in \reals^m$ is the embedding from which to begin the trajectory.
$k \in [1, m]$ is ordinal of the eigenvector to trace, with $k = 1$ corresponding to the
leading eigenvector.
$\alpha > 0$ is the step size.
$\rho \in [0, 1)$ is the decay factor.
$N \geq 1$ is the required number of steps in the trajectory.

TraceEigenpath$\mv, z, k, \alpha, \rho, N$
$z_0 \gets z$
$i \in [1, N]$
$M_z(z) \gets \Call{EvaluateNormalJacobian}{\mv, z}$ Details in Appendix <ref>
$V, D, V^t \gets \Call{SVD}{M_z(z)}$
$w_i \gets v_k$ Take the $k$th eigenvector
$i \geq 2$
$\inner{w_{i - 1}}{w_i} < 0$
$w_i \gets -w_i$ Prevent backtracking by ensuring that
${\angle(w_{i - 1}, w_i) \leq \pi / 2}$
$w_i \gets \rho w_{i - 1} + (1 - \rho) w_i$ Apply decay to smoothen
the trajectory
$z_i \gets z_{i - 1} - \alpha w_i$
$\{z_0, \ldots, z_N\}$

[][Trajectories corresponding to the principal eigenvector for embeddings $z_1, \ldots,
z_4$, from top to bottom, respectively. For embeddings $z_1, z_2$, and $z_3$, we show iterates $0,
20, \ldots, 180$, and for embedding $z_4$, we show iterates $0, 10, \ldots, 90$. The first
trajectory varies azimuth; the second, gender; the third, hair length; and the fourth, hair

[][Trajectories corresponding to the first, second, third, and fifth leading eigenvectors
of embedding $z_5$, from top to bottom, respectively; we show iterates $0, 20, \ldots, 200$. The
first trajectory primarily varies azimuth; the second, gender; the third, age; and the fourth, hair
color and presence of facial hair.]
Trajectories found by following leading eigenvectors from five fixed embeddings $z_1,
\ldots, z_5$, using Algorithm <ref> (${\alpha = 1.5 \cdot 10^{-2}}, {\rho =
0}$). Details regarding the model architecture are given in
Appendix <ref>.

§ ALIGNING THE LOCAL DISENTANGLED FACTORS

The $\beta$-VAE [Higgins et al., 2016] is known to learn disentangled representations, so that traveling
along paths of the form

\begin{equation}
\gamma : t \mapsto z_0 + t e_j,
\label{eq:aligned_path}
\end{equation}

for certain coordinate directions $e_j$, produces changes to isolated attributes of $G(z_0)$. In the
previous section, we saw that such paths still exist for GANs, but they are not simply straight
lines oriented along the coordinate axes (see Figure <ref>). We develop an
efficient regularizer that encourages these paths to take the form of
Equation <ref>, based on alignment of the top $k$ eigenvectors of $M_z(z_0)$ with
the first $k$ coordinate directions ${e_1, \ldots, e_k}$. Before proceeding, we describe a useful
visualization technique to help measure the extent to which this happens. Let ${M_z(z_0) = V D V^t}$
be the eigendecomposition of $M_z(z_0)$, where $V$ is an orthogonal matrix whose columns are
eigenvectors, and $D$ the diagonal matrix of nonnegative eigenvalues, sorted in descending order.
Now we define ${\tilde{V} : z \mapsto V}$ to be the function that maps $z$ to the corresponding
eigenvector matrix $V$, and let

\begin{equation}
F \coloneqq \E_{z \sim p_z} \tilde{V}(z) \circ \tilde{V}(z),
\label{eq:heatmap_matrix}
\end{equation}

where $p_z$ is the prior over $Z$, and `$\circ$' denotes Hadamard product. If the $k$th column is
close to a one-hot vector, with values close to zero everywhere except at entry $j$, then we know
that on average, the $k$th leading eigenvector of $M_z(z)$ is aligned with $e_j$. A heatmap
generated from this matrix therefore allows us to gauge the extent to which each eigenvector $v_k$
is aligned with the coordinate direction $e_k$. Figure <ref> shows that this
does not happen automatically for a GAN, even when it is trained with a small latent variable size.
Interestingly, eigenvector alignment does occur for $\beta$-VAEs.
Appendix <ref> explores this connection in more detail.

We begin by considering the case where we only seek to align the leading eigenvector ${v_1 \in
\reals^m}$ with the first coordinate direction ${e_1 \in \reals^m}$. A simple way to do this is to
obtain an estimate $\hat{v}_1$ for $v_1$ using $T$ power iterations, and then
renormalize ${\hat{v}_1 \in \reals^m}$ to a unit vector. We can then maximize the value of the first
component of the elementwise squared vector ${\hat{v}_1 \circ \hat{v}_1}$, and minimize values of
the remaining components. Using the mask ${s_1 \coloneqq (-1, 1, \ldots, 1) \in \reals^m}$, we
define the regularizer ${R_1 : \reals^m \to \reals}$,

\begin{equation}
R_1(z) \coloneqq \sum_{i \in [1, m]} (s_1 \circ \hat{v}_1 \circ \hat{v}_1)_i.
\end{equation}

Since $\hat{v}_1$ is constrained to unit norm, this regularizer is bounded. It can be incorporated
into the loss for the generator using an appropriate penalty weight ${\lambda > 0}$.

Next, we consider the case where we would like to align the first two leading eigenvectors ${v_1,
v_2 \in \reals^m}$ with the first two coordinate directions ${e_1, e_2 \in \reals^m}$. One potential
approach is to first compute an estimate $\hat{v}_1$ to $v_1$ using $T$ power iterations, as before.
Then, we could apply a modified version of the power method to obtain an estimate $\hat{v}_2$
for $v_2$, in which we project the result of each power iteration onto $\vspan(\hat{v}_1)^\perp$
using the projection ${P_1 \coloneqq I - \hat{v}_1 \hat{v}_1^t}$. There are two problems with this
approach. Firstly, it could be inaccurate: unless $\|v_1 - \hat{v}_1\| < \tau$ for sufficiently
small $\tau > 0$, which may require a large number of power iterations, $P_1$ will not be an
accurate projection onto $\vspan(v_1)^\perp$. Error in approximating $v_1$ would then jeopardize the
approximation to $v_2$. Second, the approach is inefficient. We can only run the power method to
estimate $v_2$ after the we have already obtained an estimate for $v_1$. If we use $T$ power
iterations to estimate each eigenvector, then estimating the first $k$ eigenvectors will require a
total of $kT$ power iterations. This is too slow to be practical.

Our specific application of the power method enables an optimization that allows for the $k$ power
iterations to be run simultaneously. Once we apply the regularizer to the GAN training procedure,
the first eigenvector $v_1$ will quickly align with the first coordinate direction $e_1$. We
therefore assume that ${v_1 = e_1}$. This assumption would imply that ${\vspan(v_1)^\perp =
\vspan(e_2, \ldots, e_m)}$, so applying $P_1$ would amount to zeroing out the first component
of $\hat{v}_2$ after each power iteration. Since $P_1$ would no longer depend on $\hat{v}_1$, we can
run the power iterations for $v_1$ and $v_2$ in parallel. To formally describe this, we let ${c_p
\coloneqq 1 \in \reals^p}$ be the constant vector of ones, and let

\begin{equation}
M_2 \coloneqq \begin{pmatrix*} c_m & \begin{matrix*} 0 \\ c_{m - 1} \end{matrix*} \end{pmatrix*}.
\label{eq:eigvec_mask_2}
\end{equation}

Given a matrix ${\hat{V}_t^{(2)} \in \reals^{m \times 2}}$ whose columns are the current estimates
for $v_1$ and $v_2$, respectively, we can describe the power iterations for $v_1$ and $v_2$ using
the recurrence

\begin{equation}
\hat{V}_{t + 1}^{(2)} \coloneqq M_2 \circ (M_z(z_0) \hat{V}_t^{(2)}).
\label{eq:power_iter_2}
\end{equation}

Now, let $\hat{V}^{(2)} \in \reals^{m \times 2}$ be the final estimate for $v_1$ and $v_2$. To
implement the regularizer, we let $s_2 \coloneqq (1, -1, 1, \ldots, 1) \in \reals^m$, $S_2 \coloneqq
\begin{pmatrix*} s_1 & s_2 \end{pmatrix*} \in \reals^{m \times 2}$, and define

\begin{equation}
R_2(z) \coloneqq \sum_{\substack{i \in [1, m] \\ j \in [1, 2]}}
(S_2 \circ \hat{V}^{(2)} \circ \hat{V}^{(2)})_{i, j}.
\label{eq:reg_2}
\end{equation}

[t]Procedure to estimate the top $k$ eigenpairs of $M_z(z)$.

[1]
$\mv : \reals^m \to \reals^m$ is a function that computes matrix-vector products with the
implicitly-defined matrix $M_z(z) \in \reals^{m \times m}$.
$V \in \reals^{m \times k}$ is a matrix whose columns are the initial estimates for the
$T \geq 1$ is the required number of power iterations.

$\epsilon \gets 10^{-8}$ Guards against division by numbers close to zero

EstimateLeadingEigenpairs$\mv, V, T$
Let $M_k \in \reals^{m \times k}$ be given by Equation <ref>
$V_0 \gets M_k \circ V$ `$\circ$' denotes Hadamard product

$i \in [1, T]$
$V_i \gets M \circ \mv(V_{i - 1})$
$\Lambda_i \gets \diag(\Call{ColumnNorms}{V_i})$
$V_i \gets V_i\, (\Lambda_i + \epsilon I)^{-1}$ Renormalize columns

$\Lambda_T, V_T$

ColumnNorms$A$ $A \in \reals^{p \times q}$
$(\|a_1\|, \ldots, \|a_q\|)$ $a_i \in \reals^p$ is the $i$th column
of $A$

[t]Procedure to evaluate the alignment penalty.

[1]
$k \in [1, m]$ is number of leading eigenvectors to align with $e_1, \ldots, e_k$.
$\mv, T$ are as defined in Algorithm <ref>.

EvaluateAlignmentRegularizer$k, \mv, T$
Let $S_k$ be given by Equation <ref>

$\alpha \gets 2 / (k (k + 1))$
$A \gets \diag(\alpha \cdot (k, k - 1, \ldots, 1))$
$S_k \gets S_k A$ Reweight columns to prioritize alignment of leading eigenvectors

$V_0 \gets \Call{RandomRademacher}{m, k}$
$\hat{\Lambda}, \hat{V} \gets \Call{EstimateLeadingEigenpairs}{\mv, V_0, t}$
$\Call{Sum}{S_k \circ \hat{V} \circ \hat{V}}$

RandomRademacher$p, q$
$A \in \reals^{p \times q}$, where $a_{ij} = 1$ with probability $1 / 2$ and $-1$
with probability $1 / 2$

It is straightforward to generalize this approach to the case where we seek to align the first $k$
eigenvectors ${v_1, \ldots, v_k}$ with ${e_1, \ldots, e_k}$. For each eigenvector $v_i$, with ${i
\in [2, k]}$, we assume that eigenvectors ${v_1, \ldots v_{i - 1}}$ are already aligned with ${e_1,
\ldots, e_{i - 1}}$. The projections onto $\vspan(e_1)^\perp, \,\vspan(e_1, e_2)^\perp, \ldots,
\,\vspan(e_1, \ldots, e_{i - 1})^\perp$ can be implemented using columns ${2, 3, \ldots, i}$,
respectively, of the mask $M_k \in \reals^{m \times k}$. This mask, which is a generalization of the
one defined by Equation <ref>, is given by ${\reals^{m \times k} \ni M_k \coloneqq}$

\begin{equation}
\smallmath{\begin{pmatrix*}
c_k & c_k - e_1 & c_k - (e_1 + e_2) & \cdots & c_k - (e_1 + \cdots + e_{k - 1}) \\
c_{m - k} & c_{m - k} & c_{m - k} & \cdots & c_{m - k}
\end{pmatrix*}}. \label{eq:eigvec_mask}
\end{equation}

The resulting procedure to estimate the leading $k$ eigenvectors is described by
Algorithm <ref>. Figure <ref> shows that the runtime
of Algorithm <ref> scales linearly with respect to the number of
eigenvectors $k$ and the number of power iterations $T$. Next, we generalize
Equation <ref>, in order to describe how to evaluate the regularizer $R_k$. Let ${s_p \in
\reals^m}$ be given by $(s_p)_i = -1$ if $i = p$ and $1$ otherwise, and define

\begin{equation}
S_k \coloneqq \smallmath{\begin{pmatrix*}[l]
-1                & \phantom{-}1      & \phantom{-}1      & \cdots & \phantom{-}1      \\
\phantom{-}1      & -1                & \phantom{-}1      & \cdots & \phantom{-}1      \\
\phantom{-}1      & \phantom{-}1      & -1                & \cdots & \phantom{-}1      \\
\phantom{-}\vdots & \phantom{-}\vdots & \phantom{-}\vdots & \ddots & \phantom{-}\vdots \\
\phantom{-}1      & \phantom{-}1      & \phantom{-}1      & \cdots & -1                \\
\phantom{-}\vdots & \phantom{-}\vdots & \phantom{-}\vdots & \ddots & \phantom{-}\vdots \\
\phantom{-}1      & \phantom{-}1      & \phantom{-}1      & \cdots & \phantom{-}1
\end{pmatrix*}}
= \begin{pmatrix*} s_1 & \cdots & s_k \end{pmatrix*} \in \reals^{m \times k}.
\label{eq:reg_mask}
\end{equation}

Algorithm <ref> shows how $S_k$ is used with the result of
Algorithm <ref> to evaluate $R_k$. Finally,
Figure <ref>(c) shows that incorporating the alignment regularizer into the
generator loss successfully aligns the trajectories produced by Algorithm <ref>
with the coordinate axes.

Median relative cost per generator update as a function of the number of eigenvectors $k$
to align, and the number of power iterations $T$ to use for Algorithm <ref>. The
cost is measured relative to the time required to make one RMSProp [Tieleman & Hinton, 2012] update to the
parameters of a DCGAN generator (approximately $20.419$ ms) with base feature map count $64$ and
batch size $32$. Complete details regarding the model architecture are given in
Appendix <ref>. Medians were computed using the update times for the first $500$
iterations; the median absolute deviations are small enough as to be indistinguishable from the
medians on the plot.

Comparison of the top-left $8 \times 8$ corner of the
matrix $F$ (Equation <ref>) for alignment-regularized GANs (${k = 8, T = 8,
\lambda = 0.1}$), trained with reweighting (left) and without reweighting (right) of the columns of
the matrix $S_k$ defined by Equation <ref>. Both GANs were trained on the CelebA
dataset [Liu et al., 2015] at $64 \times 64$ resolution, with latent variable size 128. See
Appendix <ref> for complete details regarding the model architecture and training

If the alignment regularizer is implemented exactly as described by Equation <ref>, it will
fail to have the intended effect. The reason for this has to do with the assumption behind the
optimization used to run the $k$ power iterations in parallel. Before attempting to align $v_i$
with $e_i$, we assume that that ${v_1, \ldots, v_{i - 1}}$ are already aligned with ${e_1, \ldots,
e_{i - 1}}$. When this assumption fails to hold, the projections computed using the columns of $M_k$
will no longer be valid. Figure <ref>(a) shows that the matrix $F$ does not
have a diagonal in its top-left corner, which is what we would expect to see if ${v_1, \ldots, v_k}$
were aligned with ${e_1, \ldots, e_k}$. Fortunately, there is a simple fix that remedies the
situation. We would like to encourage the optimizer to prioritize alignment of $v_i$ with $e_i$ over
alignment of ${v_{i + 1}, \ldots, v_k}$ with ${e_{i + 1}, \ldots, e_k}$, for all $i \in [1, k - 1]$.
A simple way to do this is to multiply the $i$th column of $M_k$ by a weight of ${(k - i + 1)
\alpha}$. We choose $\alpha$ based on the condition that these weights sum to one, i.e.,

\begin{equation}
\sum_{i \in [1, k]} i \alpha = 1,
\quad\text{implying that}\quad
\alpha = \frac{2}{k (k + 1)}.
\end{equation}

This reweighting scheme is implemented in lines 5–7 of Algorithm <ref>.
Figure <ref>(b) confirms that this modification induces the desired structure
in the top-left corner of $F$.

§ RESULTS

Model     2@cDisentanglement Score

Ground truth        2@l100

Raw pixels          45.75     0.8

PCA                 84.90     0.4

ICA                 42.03     10.6

DC-IGN              99.3     0.1

InfoGAN$^1$         73.5      0.9

VAE (untrained)     44.14     2.5

VAE                 61.58     0.5

$\beta$-VAE         99.23     0.1

Ours                92.34     0.4

Comparison of disentanglement metric scores for our method to those reported
in [Higgins et al., 2016]. We use an alignment-regularized GAN (${k, \lambda = 6, 0.6}$) with latent
variable size $10$. Details regarding the model architecture and training procedure are given in
Appendix <ref>.

We first make a quantitative comparison between our approach and previous methods that have been
used to obtain disentangled representations. This requires us to measure the extent to which the
second property of disentangled representations – namely, interpretability of changes resulting
from perturbations to individual coordinates of a latent variable – holds. Suppose that we had
knowledge of the ground truth latent factors for the dataset, along with a simulator that can
synthesize new outputs given assignments to these latent factors. Then we could generate pairs of
outputs, such that the outputs in each pair differ only along a single attribute. Let ${(x_0, x_1)}$
be one such pair, and ${z_0 \coloneqq G^{-1}(x_0)}$ and ${z_1 \coloneqq G^{-1}(x_1)}$ the
corresponding latent variables for the generator ${G : \reals^m \to \reals^n}$. If $G$ satisfies the
second property, then we would expect $z_0$ and $z_1$ to be approximately equal along all
components, except the one corresponding to the attribute that was varied. Hence, ${|z_0 -
z_1|}$ should be a one-hot vector in expectation. [Higgins et al., 2016] propose an evaluation metric based
on this idea. It involves training a linear classifier to predict the index of the latent factor
that was varied in order to generate each pair. At each step of training for the classifier, we
sample a batch of input-target pairs in accordance with the procedure described in [Higgins et al., 2016],
and update the classifier using the cross-entropy loss.

Application of the evaluation metric to GANs is complicated by the fact that a direct procedure to
invert the generator is usually not available. Models such as the $\beta$-VAE consist of an encoder
that effectively functions as an inverse for the generative model, so this is not an issue. For the
purpose of evaluation, we also train an encoder to invert the fixed generator, after GAN training
has finished. This additional training procedure uses a standard autoencoding loss, and the details
are specified in Appendix <ref>.
Table <ref> compares our approach to the ones evaluated
in [Higgins et al., 2016] on the dSprites dataset [Matthey et al., 2017], which we describe in
Section <ref>. As stated in Appendix <ref>, we added noise to both
the real and fake inputs of the discriminator in order to stabilize GAN training on this dataset,
for which the pixels are binary-valued. Incorporating this modification into the training procedure
for InfoGAN may improve the score reported by [Higgins et al., 2016]. Our approach achieves a score
competitive to that of DC-IGN, which makes use of supervised information. We note that the encoder
trained to invert the generator does not always succeed in finding a latent variable that yields an
accurate reconstruction of the input. This may lead to a reduction in the disentanglement score, a
problem that does not affect the other approaches for which an accurate inference procedure is

Next, we make a qualitative comparison between our approach and $\beta$-VAE. We train a series
of GANs on the CelebA dataset [Liu et al., 2015], with $k \in \{{8, 16, 32}\}$ for
Algorithm <ref> and varying values for the alignment regularizer weight $\lambda$.
We also train a series of $\beta$-VAEs with ${\beta \in \{1, 2, 4, 8, 16\}}$. The results for our
method with ${k, \lambda = 32, 0.5}$ are shown in Figure <ref>, and the results
for $\beta$-VAE with $\beta = 8$ in Figure <ref>.
Figure <ref> shows the top-left corners of the matrix $F$ given by
Equation <ref>, for three configurations of our approach with $k \in \{{8, 16,
32\}}$. A penalty weight of $\lambda = 0.5$ was not sufficient to result in a clear diagonal
structure in the top-left corner of the matrix $F$ for the configuration with $k = 32$.
Nonetheless, this configuration resulted in the largest number of disentangled factors, and the
results are shown in Figure <ref>. The best results for the configurations
with ${k \in \{8, 16\}}$ are shown in Appendix <ref>. We found that
the $\beta$-VAE configuration with ${\beta = 8}$ resulted in the largest number of disentangled
factors, while still maintaining sufficiently low reconstruction error so as to keep the sample
quality acceptable. These results are shown in Figure <ref>. In addition to
better sample quality, our approach is able to learn concepts such as different kinds of hair styles
(coordinates 9 and 10 in Figure <ref>) that are not modeled by the $\beta$-VAE.

Even with relatively large values for the penalty weight, the heatmaps shown in
Figure <ref> contain nonzero entries above and below the diagonal. This can
result in some degree of leakage of the attribute controlled by one coordinate into the next. To see
this in more detail, we examine the disentangled factors found by the configuration with $k = 16$,
which are shown in Figure <ref>. Coordinates ${\{1, 2, 3\}}$ all involve change
in hair darkness; coordinates ${\{4, 5, 6\}}$ all involve change in gender; coordinates $\{{6, 7, 8,
9}\}$ all involve smiling; and coordinates $\{{12, 13\}}$ all involve change in the location of the
hair partition. This is a limitation of our current approach, and modifications to the
implementation of the alignment regularizer in Algorithm <ref> may help to
mitigate this leakage. We plan to investigate such improvements in future work.

[][Top-left ${8 \times 8}$ corner of $F$ for GAN with ${k, \lambda = 8, 2.3}$.]

[][Top-left ${16 \times 16}$ corner of $F$ for GAN with ${k, \lambda = 16, 1.7}$.]

[][Top-left ${21 \times 21}$ corner of $F$ for GAN with ${k, \lambda = 32, 0.5}$.]

Comparison of the top-left corners of the matrix $F$ (Equation <ref>) for
alignment-regularized GANs with ${k \in \{8, 16, 32\}}$. All GANs were trained with latent variable
size $128$. Details regarding the model architecture and training procedure are given in
Appendix <ref>.

Coordinate     Description     Example

1      Background darkness            ⊷0.jpg

2      Azimuth                        ⊷1.jpg

3      Bangs                          ⊷2.jpg

4      Gender, bangs, and smiling     ⊷3.jpg

5      Smiling                        ⊷4_v2.jpg

6      Hair color                     ⊷5.jpg

7      Hair color and hair style      ⊷6_v2.jpg

8      Lighting color and bangs       ⊷7.jpg

9      Hair color and hair style      ⊷8_v2.jpg

10     Hair color and hair style      ⊷9_v2.jpg

11     Jawline                        ⊷10.jpg

12     Smiling and bangs              ⊷11.jpg

13     Background color               ⊷12.jpg

14     Age and hairline               ⊷13.jpg

15     Age                            ⊷14_v2.jpg

17     Location of hair partition     ⊷15.jpg

18     Lighting and skin tone         ⊷16.jpg

20     Mouth open                     ⊷19.jpg

21     Mouth open                     ⊷20.jpg

Disentanglement results for alignment-regularized GAN (${k, \lambda = 32, 0.5}$) with
latent variable size $128$ on the CelebA dataset, at ${64 \times 64}$ resolution. Details regarding
the model architecture and training procedure are given in
Appendix <ref>. Results with additional samples for each coordinate can be found at
this URL: <https://drive.google.com/open?id=1E2TneLSAdyFYgN4GUzYKnHDYABo-G8RK>.

Coordinate     Description     Example

1      Azimuth                             ⊷0.jpg

2      Background color                    ⊷1.jpg

3      Location of hair partition          ⊷2.jpg

4      Skin color                          ⊷3.jpg

5      Hair direction                      ⊷4.jpg

8      Background color                    ⊷7.jpg

12     Hair length                         ⊷11.jpg

14     Jawline                             ⊷13.jpg

19     Skin tone                           ⊷18.jpg

20     Hair style                          ⊷19.jpg

22     Sunglasses                          ⊷21.jpg

23     Gender                              ⊷22.jpg

25     Background darkness                 ⊷24.jpg

26     Hairline and skin tone              ⊷25.jpg

28     Lighting direction and rotation     ⊷27.jpg

32     Hair size, color, and smiling       ⊷31.jpg

Disentanglement results for $\beta$-VAE (${\beta = 8}$) with latent variable size $32$ on
the CelebA dataset, at ${64 \times 64}$ resolution. Details regarding the model architecture and
training procedure are given in Appendix <ref>. Results with additional samples for
each coordinate can be found at this URL:

§ CONCLUSION

Our work approaches the problem of learning disentangled representations from the perspective of
eigenvector alignment. We develop a novel regularizer which, when incorporated into the GAN
objective, induces disentangled representations of quality comparable to those obtained by VAE-based
approaches (Higgins et al., 2016, Kim & Mnih, 2018, Chen et al., 2018, Esmaeili et al., 2018), without the need to
introduce auxiliary models into the training procedure [Chen et al., 2016]. This approach is also not
specific to the GAN framework [Goodfellow et al., 2014], and could potentially be applied to autoregressive
models, such as those used to generate text and audio. We believe this is an important direction for
future investigation. So far, two different perspectives for viewing disentanglement have been
proposed: maximizing mutual information and eigenvector alignment. An investigation into the
relationship between the two could further our understanding of what disentanglement is and why it

§ ACKNOWLEDGEMENTS

We would like to thank for Ryan-Chris Moreno Vasquez and Emily Denton for early discussions that led
us to think carefully about the paths determined by the leading eigenvectors of $M_z(z_0)$. We are
also grateful for the suggestions provided by Mikael Henaff, Alec Radford, Yura Burda, and Harri
Edwards, which improved the quality of our presentation.

[Abadi et al., 2016]
Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., Dean, J., Devin, M.,
Ghemawat, S., Irving, G., Isard, M., et al.
Tensorflow: A system for large-scale machine learning.
In OSDI, volume 16, pp.  265–283, 2016.

[Aubry et al., 2014]
Aubry, M., Maturana, D., Efros, A. A., Russell, B. C., and Sivic, J.
Seeing 3d chairs: exemplar part-based 2d-3d alignment using a large
dataset of cad models.
In Proceedings of the IEEE conference on computer vision and
pattern recognition, pp.  3762–3769, 2014.

[Baydin et al., 2017]
Baydin, A. G., Pearlmutter, B. A., Radul, A. A., and Siskind, J. M.
Automatic differentiation in machine learning: a survey.
Journal of Machine Learning Research, 180
(153):0 1–153, 2017.

[Brock et al., 2018]
Brock, A., Donahue, J., and Simonyan, K.
Large scale gan training for high fidelity natural image synthesis.
arXiv preprint arXiv:1809.11096, 2018.

[Burgess et al., 2018]
Burgess, C. P., Higgins, I., Pal, A., Matthey, L., Watters, N., Desjardins, G.,
and Lerchner, A.
Understanding disentangling in $\beta$-vae.
arXiv preprint arXiv:1804.03599, 2018.

[Chen et al., 2018]
Chen, T. Q., Li, X., Grosse, R., and Duvenaud, D.
Isolating sources of disentanglement in variational autoencoders.
arXiv preprint arXiv:1802.04942, 2018.

[Chen et al., 2016]
Chen, X., Duan, Y., Houthooft, R., Schulman, J., Sutskever, I., and Abbeel, P.
Infogan: Interpretable representation learning by information
maximizing generative adversarial nets.
In Advances in neural information processing systems, pp.  2172–2180, 2016.

[Esmaeili et al., 2018]
Esmaeili, B., Wu, H., Jain, S., Bozkurt, A., Siddharth, N., Paige, B., Brooks,
D. H., Dy, J., and van de Meent, J.-W.
Structured disentangled representations.
stat, 1050:0 29, 2018.

[Goodfellow et al., 2014]
Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair,
S., Courville, A., and Bengio, Y.
Generative adversarial nets.
In Advances in neural information processing systems, pp.  2672–2680, 2014.

[Higgins et al., 2016]
Higgins, I., Matthey, L., Pal, A., Burgess, C., Glorot, X., Botvinick, M.,
Mohamed, S., and Lerchner, A.
beta-vae: Learning basic visual concepts with a constrained
variational framework.

[Karras et al., 2017]
Karras, T., Aila, T., Laine, S., and Lehtinen, J.
Progressive growing of gans for improved quality, stability, and
arXiv preprint arXiv:1710.10196, 2017.

[Kim & Mnih, 2018]
Kim, H. and Mnih, A.
Disentangling by factorising.
arXiv preprint arXiv:1802.05983, 2018.

[Kingma & Ba, 2014]
Kingma, D. P. and Ba, J.
Adam: A method for stochastic optimization.
arXiv preprint arXiv:1412.6980, 2014.

[Kingma & Welling, 2013]
Kingma, D. P. and Welling, M.
Auto-encoding variational bayes.
arXiv preprint arXiv:1312.6114, 2013.

[Liu et al., 2015]
Liu, Z., Luo, P., Wang, X., and Tang, X.
Deep learning face attributes in the wild.
In Proceedings of International Conference on Computer Vision
(ICCV), 2015.

[Matthey et al., 2017]
Matthey, L., Higgins, I., Hassabis, D., and Lerchner, A.
dsprites: Disentanglement testing sprites dataset.
https://github.com/deepmind/dsprites-dataset/, 2017.

[Mescheder et al., 2018]
Mescheder, L., Geiger, A., and Nowozin, S.
Which training methods for gans do actually converge?
In International Conference on Machine Learning, pp.  3478–3487, 2018.

[Miyato et al., 2018]
Miyato, T., Kataoka, T., Koyama, M., and Yoshida, Y.
Spectral normalization for generative adversarial networks.
arXiv preprint arXiv:1802.05957, 2018.

[Radford et al., 2015]
Radford, A., Metz, L., and Chintala, S.
Unsupervised representation learning with deep convolutional
generative adversarial networks.
arXiv preprint arXiv:1511.06434, 2015.

[Rezende et al., 2014]
Rezende, D. J., Mohamed, S., and Wierstra, D.
Stochastic backpropagation and approximate inference in deep
generative models.
arXiv preprint arXiv:1401.4082, 2014.

[Salimans & Kingma, 2016]
Salimans, T. and Kingma, D. P.
Weight normalization: A simple reparameterization to accelerate
training of deep neural networks.
In Advances in Neural Information Processing Systems, pp.  901–909, 2016.

[Tieleman & Hinton, 2012]
Tieleman, T. and Hinton, G.
Lecture 6.5—RmsProp: Divide the gradient by a running average of
its recent magnitude.
COURSERA: Neural Networks for Machine Learning, 2012.

[Townsend, 2017]
Townsend, J.
A new trick for calculating jacobian vector products.
<https://j-towns.github.io/2017/06/12/A-new-trick.html>, 2017.

[Xiang & Li, 2017]
Xiang, S. and Li, H.
On the effects of batch and weight normalization in generative
adversarial networks.
stat, 1050:0 22, 2017.

§ FORWARD- AND REVERSE-MODE AUTOMATIC DIFFERENTIATION

The entries of $\jg$ are not explicitly stored in memory: it is an implicit-defined matrix.
As such, information about $\jg$ must be accessed via automatic differentiation (AD), of which there
are two kinds: forward-mode and reverse-mode. We briefly describe them here, and refer the reader
to [Baydin et al., 2017] for a more comprehensive survey.

Suppose that we are given a differentiable function ${f : \reals^m \to \reals^n}$ corresponding to a
feedforward neural network with $L$ layers, and let ${G \coloneqq (V, E)}$ be its representation as
a computation graph. For each ${k \in [1, L]}$, let $R_k$ be the set of nodes corresponding to
the $k$th layer, and let ${r_k : \reals^{\alpha_{k - 1}} \to \reals^{\alpha_k}}$ be the function
computed by the nodes in $R_k$, with ${\alpha_0 \coloneqq m}$ and ${\alpha_L \coloneqq n}$. Then ${f
= r_L \circ r_{L - 1} \circ \cdots \circ r_2 \circ r_1}$. Finally, let ${\bar{x} \in \reals^m}$
denote the variable for the input to $f$, and ${\bar{b}_{k - 1} \in \reals^{\alpha_{k - 1}}}$ the
variable for the input to $r_k$, with ${\bar{b}_0 \coloneqq \bar{x}}$. Given a vector ${v \in
\reals^m}$, forward-mode AD computes $Jv$. It works according to the recurrence

\begin{align*}
v_k \coloneqq \left( D_{\bar{x}} r_k \big|_{\bar{x} = x} \right) v
&= \left( D_{\bar{b}_{k - 1}} r_k \big|_{\bar{b}_{k - 1} = b_{k - 1}} \right)
\left( D_{\bar{x}}\; r_{k - 1} \circ \ldots \circ r_1 \big|_{\bar{x} = x} \right) v \\
&= \left( D_{\bar{b}_{k - 1}} r_k \big|_{\bar{b}_{k - 1} = b_{k - 1}} \right) v_{k - 1},
\end{align*}

where $k \in [1, L]$ and $v_0 \coloneqq v$. After step $k$, we will have obtained the product of the
Jacobian of $r_k$ with respect to $x$, and $v$; after step $L$, we will have obtained the desired
product $Jv$.

Forward-mode AD computes $Jv$ far more efficiently than the approach of first evaluating $J$, and
subsequently multiplying by $v$. Suppose for simplicity that ${\alpha_k = n}$ for all ${k \in [0,
L]}$, so that ${J \in \reals^{n \times n}}$. To compute $J$ independently, we would use the chain
rule, which gives

\[
D_{\bar{x}} f \big|_{\bar{x} = x} =
\left( D_{\bar{b}_{L - 1}} r_L \big|_{\bar{b}_{L - 1} = b_{L - 1}} \right)
\left( D_{\bar{b}_{L - 2}} r_{L - 1} \big|_{\bar{b}_{L - 2} = b_{p - 2}} \right)
\cdots \\
\left( D_{\bar{b}_0} r_1 \big|_{\bar{b}_0 = x} \right).
\]

Each layer $r_k$, for ${k \in [1, L]}$, must compute its $n \times n$ Jacobian, and each layer after
the first must multiply its Jacobian with the $n \times n$ Jacobian of the preceding layer with
respect to $x$. This process is described by the recurrence

\[
D_{\bar{x}} r_k \big|_{\bar{x} = x} =
\left( D_{\bar{b}_{k - 1}} r_k \big|_{\bar{b}_{k - 1} = b_{k - 1}} \right)
\left( D_{\bar{x}} r_{k - 1} \big|_{\bar{x} = x} \right).
\]

Assuming that $\Theta(1)$ operations are required to compute each element of the Jacobian, the total
number of operations required is proportional to

\[
N \coloneqq n^2 + (L - 1) (n^2 + n^3) \in \Theta(n^3),
\]

despite the fact that $J$ has only $n^2$ elements. On the other hand, the $k$th step of the
forward-mode recurrence requires $\Theta(n)$ operations, so computing $Jv$ only
requires $\Theta(Ln)$ operations. By running forward-mode with ${v = e_i}$ for each $i \in [1, n]$,
we can form $J$ column-by-column in only $\Theta(Ln^2)$ operations. Henceforth, we will measure
operation count in terms of invocations to the AD engine, rather than by counting elementary

Of the two kinds of AD, it is reverse-mode that finds the most use in machine learning. Given a
vector ${w \in \reals^n}$, it computes $J^t w$. Most applications involve minimizing a scalar-valued
loss function with respect to a vector of parameters, which corresponds to the case where $n = 1$.
This special case is otherwise known as backpropagation. Unlike forward-mode, which begins at the
first layer and ends at the last, reverse-mode begins at the last layer and ends at the first. It
works according to the recurrence

\begin{align*}
w_k^t \coloneqq w^t \left( D_{\bar{b}_{k - 1}} r_k \big|_{\bar{b}_{k - 1} = b_{k - 1}} \right)
&= w^t \left( D_{\bar{b}_k}\; r_L \circ \cdots \circ r_{k + 1} \big|_{\bar{b}_k = b_k} \right)
\left( D_{\bar{b}_{k - 1}} r_k \big|_{\bar{b}_{k - 1} = b_{k - 1}} \right) \\
&= w_{k + 1}^t \left( D_{\bar{b}_{k - 1}} r_k \big|_{\bar{b}_{k - 1} = b_{k - 1}} \right), \\
\end{align*}

where $k$ ranges from $L$ to $1$, and $w_L \coloneqq w$. After $L$ steps, we will have obtained the
desired product $J^t w$. Machine learning frameworks typically expose the full interface to the
reverse-mode AD engine, rather than specializing to the case of backpropagation. E.g., in
TensorFlow [Abadi et al., 2016], one can change the value of $w$ for tf.gradients from a vector of
ones to a custom value specified by the grad_ys parameter.

[style=pythonStyle, caption=TensorFlow implementation of Jacobian-vector
operations., captionpos=b, abovecaptionskip=, label=lst:jacobian_ops]
import tensorflow as tf

def forward_gradients(ys, xs, d_xs):
v = tf.placeholder_with_default(tf.ones_like(ys), shape=ys.get_shape())
g = tf.gradients(ys, xs, grad_ys=v)
return tf.gradients(g, v, grad_ys=d_xs)

def j_v(ys, xs, vs):
return forward_gradients(ys, xs, vs)

def jt_v(ys, xs, vs):
return tf.gradients(ys, xs, vs)

def jt_j_v(ys, xs, vs):
jv = j_v(ys, xs, vs)
return tf.gradients(ys, xs, jv)

def j_jt_v(ys, xs, vs):
jt_v_ = jt_v(ys, xs, vs)
return forward_gradients(ys, xs, jt_v_)

If desired, we can compute the entire Jacobian using either type of AD. By running forward-mode
with $v = e_i$ for each $i \in [1, m]$, we can form $J$ column-by-column, using $m$ total
invocations to AD. Similarly, by running reverse-mode with $w = e_i$ for each $i \in [1, n]$, we can
form $J$ row-by-row using $n$ total invocations to AD. If $n > m$, the former is typically faster;
otherwise, the latter is preferable. We note that while AD offers a relatively efficient approach
for evaluating $J$, doing so at each iteration of optimization becomes impractical.

Many popular machine learning frameworks (e.g., TensorFlow) do not implement forward-mode natively,
since it is seldom used for machine learning. Surprisingly, this is not a limitation:
reverse-mode can be used to implement forward-mode. Given a differentiable function ${f :
\reals^m \to \reals^n}$, we can compute $w^t J$ for a given vector $w \in \reals^n$, using
reverse-mode AD. Treating the input to $f$ as a constant, we can regard $w^t J$ as a function ${g :
\reals^n \to \reals^m}$, $w \mapsto w^t J$. The derivative of $g$ with respect to $w$ is given
by $J^t$, and so another application of reverse-mode AD allows us to compute ${v^t J^t = (Jv)^t}$.
Hence, reverse-mode AD can be used to implement forward-mode AD. This trick was first described
by [Townsend, 2017]. We provide a TensorFlow implementation of the procedures to compute $Jv,
J^t v, J^t J v$, and $J J^t v$ in Listing <ref>.

§ EIGENVECTOR ALIGNMENT FOR $\BETA$-VAES

Comparison between the average squared eigenvector matrices $F$
(Equation <ref>) for a GAN generator with latent variable size $16$ (left) and the
decoder from a $\beta$-VAE (${\beta = 8}$) with latent variable size $32$ (right), both trained on
the CelebA dataset. Complete architecture and training details are provided in
Appendix <ref>. The matrix $F$ for the VAE decoder contains several columns that are
close to one-hot vectors. By contrast, the same matrix for the GAN generator exhibits little
structure, despite the fact that it is trained with a small latent variable

Investigation of whether the matrix $F$ can be used to determine which directions in latent
space correspond to disentangled factors for a $\beta$-VAE. For each $i \in [1, 32]$, we compute the
maximum value along the $i$th row of matrix $F$ shown in Figure <ref>(b).
Then, we plot a point indicating whether or not disentanglement occurs along this direction, as
determined by visual inspection. The directions that do result in disentanglement are shown in
Figure <ref>.

The results from Section <ref> suggest that alignment of the eigenvectors
of $M_z(z_0)$ with the coordinate axes might be sufficient to induce disentanglement in the latent
representation of a generative model. The $\beta$-VAE is known to learn such a representation
when $\beta$ is increased, so that the KL-divergence between the approximate posterior and the prior
of the decoder is made sufficiently small. Suppose that a $\beta$-VAE exhibits disentanglement along
the $j$th coordinate direction in latent space. Then paths of the form ${\gamma : t \mapsto z_0 + t
e_j}$ will produce changes to $G(z_0)$ along isolated factors of variation. If such a path coincides
with the trajectory found by Algorithm <ref> for the $k$th leading eigenvector
throughout latent space, then this eigenvector must be aligned with the $j$th coordinate axis.
Hence, we would expect the $k$th column of the matrix $F$ given by Equation <ref>
to be a one-hot vector close to $e_j$. Figure <ref> shows that, in fact,
several columns of $F$ have high similarity to coordinate directions.

Next, we investigate whether a column of $F$ having high similarity with a coordinate
direction $e_j$ implies that disentanglement actually occurs along ${\gamma : t \mapsto z_0 + t
e_j}$. Figure <ref> shows that with the exception of three
coordinates having similarity greater than 0.35, this turns out not to be the case. In other words,
several of the eigenvectors naturally align with coordinate directions during training, but
disentanglement does not reliably occur along all of them. More work needs to be done in order to
better understand the relationship between eigenvector alignment and disentangled representations.

§ GENERATOR ARCHITECTURES AND TRAINING PROCEDURE

Figure     DCGAN Base Feature Map Count     Latent Variable Size     Notes

<ref>(c)              64      32      1

<ref>(e)              64      64      1

<ref>(g)              64      128     1

<ref>                       64      128     1

<ref>                       128     128     1

<ref>(f)              64      256     1

<ref>(h)              64      512     1

<ref>                  64      3       1, 3

<ref>      64      10      1, 3

<ref>(a)                64      16      1

<ref>     64      32      2

GAN and $\beta$-VAE architectures used for all figures and tables. Note 1: weight
normalization [Salimans & Kingma, 2016] was applied both to the generator and the discriminator, with the
scale $g$ fixed to one for the discriminator. Note 2: spectral weight
normalization [Miyato et al., 2018] was applied both to the encoder and the decoder, with
learned scales for both. Note 3: gaussian noise with standard deviation $0.6$ was added to both real
and fake inputs to the discriminator.

All models in this work are based on the DCGAN [Radford et al., 2015] architecture.
Table <ref> describes the architectures of the generator and discriminator (in the
case of GANs) and the encoder and decoder (in the case of VAEs) associated with each figure and
table. All models use the translated PReLU activation function [Xiang & Li, 2017]; the ReLU leaks are
learned, and clipped to the interval $[0, 1]$ after each parameter update. The only data
preprocessing we applied was to scale the pixel values to the interval $[0, 1]$. The GANs were
trained using the original, non-saturating GAN loss described in [Goodfellow et al., 2014] with multivariate normal
prior. Both the generator and the discriminator were trained using RMSProp [Tieleman & Hinton, 2012] with step
size $10^{-4}$, decay factor $0.9$, and ${\epsilon = 10^{-6}}$. Each parameter update was made using
a batch size of $32$, and the models were trained for a total of 750000 updates. The VAEs were
trained using a gaussian likelihood model for the decoder, and the log-diagonal covariance
parameterization for the encoder. The fixed, per-pixel standard deviation of the decoder was chosen
such that, disregarding constant terms, the log-likelihood corresponds to the reconstruction error,
normalized by the product of the number of channels and the number of pixels. To optimize the
evidence lower bound, we used Adam [Kingma & Ba, 2014] with ${\beta_1, \beta_2, \epsilon = 0.5, 0.99,
10^{-8}}$. Each parameter update was made using a batch size of $32$, and the models were trained
for a total of 1500000 updates. The KL-divergence weight $\beta$ was increased linearly
from $0$ to the final value over the first 1000000 updates. We applied the function $x \mapsto
2 (\tanh(x) + 1 / 2)$ to the outputs of both the GAN generators and the VAE decoders; this ensures
that the output pixel values are in the interval $[0, 1]$. To predict the mean and log-diagonal
covariance with the VAE encoders, we apply the translated ReLU activation function [Xiang & Li, 2017] to
the final convolutional features, followed by two fully-connected layers, one for each statistic.

§ IMPLEMENTATION DETAILS FOR DISENTANGLEMENT METRIC SCORE

[t]Procedure to generate a batch for disentanglement metric

$n_\text{inst} \geq 1$ is the number of aligned pairs to use, in order to create each
instance in the batch.
$n_\text{batch} \geq 1$ is the batch size.
$\enc : \reals^n \to \reals^m$ is the encoder functioning as the inverse of the
generator $G$.

$c \gets (3, 6, 40, 30, 30)$ Ranges for the five attributes that determine each

SampleShape$i, v$ Sample shape with attribute $i$ fixed to value $v$
$j \in [1, 5]$
$u_j \gets \Call{RandomUniform}{1, c_j}$
$u \gets (u_1, u_2, u_3, u_4, u_5)$
$u_i \gets v$

MakeInstance$n_\text{inst}, \enc$
$i \gets \Call{RandomUniform}{2, 5}$ Sample index of attribute to fix
$v \gets \Call{RandomUniform}{1, c_i}$ Sample value for fixed attribute
$z \gets 0 \in \reals^m$
$k \in [1, n_\text{inst}]$
$z_0 \gets \enc(\Call{SampleShape}{i, v})$
$z_1 \gets \enc(\Call{SampleShape}{i, v})$
$z \gets z + |z_0 - z_1|$
$z / n_\text{inst}, i$

MakeBatch$n_\text{inst}, n_\text{batch}, \enc$
$\text{inputs} \gets \varnothing$
$\text{targets} \gets \varnothing$
$k \in [1, n_\text{batch}]$
$x, y \gets \Call{MakeInstance}{n_\text{inst}, \enc}$
$\text{inputs} \gets \text{inputs} \cup \{x\}$
$\text{targets} \gets \text{targets} \cup \{y\}$
$\text{inputs}, \text{targets}$

Algorithm <ref> describes the procedure used to generate the batches to train
and evaluate the classifier for the disentanglement metric [Higgins et al., 2016]. We use ${n_\text{inst},
n_\text{batch} = 64, 32}$, and update the classifier using SGD with nesterov accelerated
gradient (step size $10^{-2}$, momentum $0.99$). We train the classifier using a total
of 10000 parameter updates, and evaluate its performance using 5000 instances, as
reported by [Higgins et al., 2016]. To invert the pretrained GAN generator, we train a VAE decoder with
twice the base feature map count for the generator, using 30000 parameter updates. The details
for this training procedure are identical to those for regular VAE training described in
Appendix <ref>, except that the generator is not updated, and the KL-divergence
weight $\beta$ is set to zero. In other words, we use a standard autoencoding loss with a fixed

§ ADDITIONAL GAN DISENTANGLEMENT RESULTS

Coordinate     Description     Example

1     Background darkness and hair color     ⊷0.jpg

2     Azimuth, lighting, and hair color      ⊷1.jpg

3     Hairline and hair color                ⊷2.jpg

4     Azimuth                                ⊷3.jpg

5     Shadow                                 ⊷4.jpg

6     Smiling, age, skin tone, gender        ⊷5.jpg

7     Smiling, age, jawline                  ⊷6.jpg

8     Jawline and hairstyle                  ⊷7.jpg

Disentanglement results for alignment-regularized GAN (${k, \lambda = 8,
2.3}$). Details regarding the model architecture and training procedure are given in
Appendix <ref>. Results with additional samples for each coordinate can be found at
this URL: <https://drive.google.com/open?id=1Wnd5jIxopFsBRlylMUN2HfWhICXGiU3u>.

Coordinate     Description     Example

1      Background and hair darkness                       ⊷0.jpg

2      Azimuth and hair darkness                          ⊷1.jpg

3      Hair length, hair darkness, and lighting           ⊷2.jpg

4      Smiling and gender                                 ⊷3.jpg

5      Hair length, hair darkness, and gender             ⊷4.jpg

6      Smiling, bangs, and gender                         ⊷5.jpg

7      Smiling and hairstyle                              ⊷6.jpg

8      Smiling, jawline, and glaring expression           ⊷7.jpg

9      Smiling, bangs, and mouth open                     ⊷8.jpg

10     Hairline                                           ⊷9.jpg

11     Raised eyebrows and skin tone                      ⊷10.jpg

12     Raised eyebrows and location of hair partition     ⊷11_v2.jpg

13     Raised eyebrows and location of hair partition     ⊷12.jpg

14     Lighting color                                     ⊷13.jpg

Disentanglement results for alignment-regularized GAN (${k, \lambda = 16,
1.7}$). Details regarding the model architecture and training procedure are given in
Appendix <ref>. Results with additional samples for each corodinate can be found at
this URL: <https://drive.google.com/open?id=1BM-P-hMF7sV_0smNFD_iTAKpq6FeNsCo>.

§ SUPPLEMENTARY FIGURES

[][Log spectra for CelebA models with $n_f = 64$ and $m$ varied.]

[][Log spectra for LSUN Bedroom models with $n_f = 256$ and $m$ varied.]

[][CelebA ($m = 32, n_f = 64, \epsilon = 0.40$), embeddings $z_{10}, z_{11}$.]

[][LSUN Bedroom ($m = 128, n_f = 64, \epsilon = 0.80$), embeddings $z_{5}, z_{12}$.]

[][CelebA ($m = 64, n_f = 64, \epsilon = 0.65$), embeddings $z_6, z_{12}$.]

[][LSUN Bedroom ($m = 256, n_f = 64, \epsilon = 0.80$), embeddings $z_6, z_8$.]

[][CelebA ($m = 128, n_f = 64, \epsilon = 0.80$), embeddings $z_7, z_9$.]

[][LSUN Bedroom ($m = 512, n_f = 64, \epsilon = 0.80$), embeddings $z_8, z_{10}$.]

Effect of perturbing a latent variable $z$ along the leading eigenvectors of $M_z(z)$.
Subfigures (a) and (b) show the top eigenvalues of $M_z(z_i)$, where $\{z_1, \ldots, z_{12}\}$ are
fixed embeddings; small eigenvalues are omitted. Subfigures (c)–(d) compare the effects of
perturbations along random directions to perturbations along leading eigenvectors. Each subfigure
consists of two stacked two-row grids. The leftmost images of each grid are both identical and equal
to $G(z_i)$, for some $i \in [1, 12]$. The first row shows the effect of perturbing $z_i$ along
13 directions sampled uniformly from the sphere of radius $\epsilon$, while the second row shows the
effect of perturbing $z_i$ along the first 13 leading eigenvectors of $M_z(z_i)$, by the same
distance $\epsilon$. Details regarding the model architectures are given in
Appendix <ref>.


Background

Mixing a fluid in some domain $X$ by chaotic advection, essentially stirring,
typically concerns the study of time-$T$ periodic incompressible laminar flows
$\mathbf{v}(\textbf{x},t)=\mathbf{v}(\textbf{x},t+T)$, $\textbf{x}\in X$. The
dynamical features of these flows are described by the trajectories of fluid
particles within the system, in particular their positions after each time
period $T$. This defines a map $f:X\rightarrow X$ which sends the initial
position of a particle to its position after time $T$. The long term behaviour
of the flow $\mathbf{v}$ on $X$ is then described by repeated iterations of
the map $f$, denoted by $f^{n}$, and the incompressibility condition on
$\mathbf{v}$ tells us that $f$ preserves the Lebesgue measure $\mu$ on $X$.
Ergodic theory provides the mathematical framework for understanding the long
term behaviour of iterating such maps, most importantly it gives a precise
definition for what it means for a map to be mixing: A $\mu$-preserving
invertible map $f:X\rightarrow X$ is _mixing_ if $\limn\mu(f^{n}(B)\cap
A)=\mu(A)\mu(B)$ for all measurable $A,B\subset X$, $\mu(X)=1$. A hallmark of
chaotic advection in two dimensions is an iterative ‘stretching and folding’
type action. Starting with a blob of fluid $B$, this action transforms
$f^{n}(B)$ into a thin fluid filament which is repeatedly folded to spread
across the entire domain, or any target set $A$. A natural model for this
behaviour is to compose shear maps, which stretch while preserving $\mu$, on
the torus $\tor$, whose periodic boundaries interweave the long fluid
filaments. A canonical example is the Cat Map $H:\tor\rightarrow\tor$
(arnold_ergodic_1968). Parameterise $\tor$ by
$(x,y)\in(\mathbb{R}/\mathbb{Z})^{2}$, then $H=G\circ F$ composes horizontal
and vertical shears, written in matrix form as

$\pmatrix{x}\\\ y\mapsto\underbrace{\pmatrix{1}&0\\\
11}_{DG}\underbrace{\pmatrix{1}&1\\\ 01}_{DF}\pmatrix{x}\\\
y\text{mod1}=\underbrace{\pmatrix{1}&1\\\ 12}_{M}\pmatrix{x}\\\ y\text{mod1},$

where $DF$, $DG$ denote the Jacobians of the maps $F$, $G$. Since $H$ can be
defined using a single hyperbolic matrix $M$, the map is _uniformly
hyperbolic_ , with the same magnitude and directions of expansion, contraction
across the entire domain. While this allows for a straightforward proof of the
mixing property, mixing in a realistic fluid flow is typically non-uniform due
to the influence of walls and non-linear velocity profiles. The key barriers
to mixing in this setting are _elliptic islands_ , invariant subsets within
which particle paths trace out closed curves. From a fluids perspective these
curves form material lines in the flow which particles cannot penetrate
(except by diffusion), leading to poor mixing (ottino_kinematics_1989). An
illustration of an island pair is given later in Figure 1. Previous studies
which try to minimise island structures in realistic mixing flows include
[franjione_symmetry_1992], looking at eggbeater and duct flows, and
[hertzsch_dna_2007], looking at mixing in DNA microarrays. Several mappings
exist which incorporate realistic flow phenomena and still allow for a proof
of the mixing property. Linked Twist Maps (burton_ergodicity_1980), hereafter
LTMs, compose monotonic shears which act on annuli of the torus, leaving a
region invariant which models a boundary within the domain. Mixing properties
can be shown in the co-rotating case, where the shears act in the positive $x$
and $y$ directions (wojtkowski_linked_1980). In the counter-rotating case,
where the direction of one shear is reversed, island structures develop and
can only be broken up by taking strong shears. This highlights a potential
challenge for mixing by non-monotonic shears, which inherently exhibit this
counter rotating quality. Indeed, there are few examples of non-monotonic
toral maps with proven mixing properties. cerbelli_continuous_2005’s map
(hereafter the CG Map), studied in [cerbelli_continuous_2005],
[mackay_cerbelli_2006], and [cerbelli_characterization_2008], incorporates
non-monotonicity into the first shear by taking

$F(x,y)=\cases{(}x+2y,y)\text{mod1}&\text{for}y\leq\frac{1}{2},\\\
(x+2(1-y),y)\text{mod1}\text{for}y\geq\frac{1}{2},$

and leaves $G$ unchanged. While this introduces a non-hyperbolic derivative
matrix over half the domain, a unique geometric feature of the map ensures
that this does not compromise long term stretching behaviour. Indeed,
hyperbolic and mixing properties of the CG Map can be proven by quite direct
means, with analysis of the unstable foliation revealing that fluid filaments
get stretched and folded in a very regimented fashion. This is not the case
for generic non-monotonic shears $F$, as shown in
[myers_hill_continuous_2021], where shears of the form

$F_{\eta}(x,y)=\cases{\left}(x+\frac{1}{1-\eta}y,y\right)\text{mod1}&\text{for}y\leq
1-\eta,\\\ \par\left(x+\frac{1}{\eta}(1-y),y\right)\text{mod1}\text{for}y\geq
1-\eta,\\\ $ (1)

are considered with $G$ left unchanged. An illustration of this shear is given
in Figure 1; note that parameter values $\eta=0$ and $\eta=1/2$ give the Cat
Map and CG Map respectively. Mixing properties over subsets of the parameter
space $0<\eta<1/2$ are shown using a scheme from [katok_invariant_1986],
Theorem LABEL:thm:katok-strelcyn in the present work, which gives
comparatively easy to verify conditions under which non-uniformly hyperbolic
systems111In particular non-uniformly hyperbolic systems with singularities,
which are inherent to piecewise linear maps. are mixing, compared to arguing
by direct means. Non-uniform hyperbolicity ensures the existence (almost
everywhere) of local stable and unstable manifolds which, roughly speaking,
describe the characteristic local flow direction in backwards and forwards
time respectively (see for example beigie_invariant_1994). The way in which
blobs of fluid are stretched and spread across the domain is then described
well by the images of these local manifolds, and mixing properties follow from
intersection conditions on these images. Broadly speaking, the key challenge
to showing these conditions for non-monotonic systems is the sign alternating
property as described in [cerbelli_continuous_2005]. While hyperbolicity
ensures that the images of local manifolds grow exponentially in length, non-
linear shears like $F_{\eta}$ fold these images back on themselves,
potentially inhibiting their spread across the domain (also commented on in
przytycki_ergodicity_1983). This challenge can often be overcome by
establishing sufficiently strong stretching behaviour over one or several
iterates, a method we will employ here to prove mixing properties of maps
composing two non-monotonic shears. While taking $G$ similar to $F_{\eta}$
(see Figure 1) means that the images of local manifolds will fold back on
themselves more often, in turn the gradients defining $G$ will become steeper,
giving stronger stretching behaviour. Alongside knowing whether a map is
mixing, it is desirable to know the rate at which we approach a mixed state.
That is, what is the rate of decay of $\left|\mu(f^{n}(B)\cap
A)-\mu(A)\mu(B)\right|$ with $n$? Taking indicator functions
$\mathbbm{1}_{A}(x)=1$ if $x\in A$, 0 otherwise, and defining the correlation
function

$C_{n}(\phi,\psi,f,\mu)=\int\left(\phi\circ
f^{n}\right)\psi\d{\mu}-\int\phi\d{\mu}\int\psi\d{\mu}$ (2)

for observables $\phi,\psi:\tor\rightarrow\mathbb{R}$, this amounts to
studying the decay of $|C_{n}(\mathbbm{1}_{B},\mathbbm{1}_{A},f^{-1},\mu)|$.
For uniformly hyperbolic maps, proving exponential correlation decay rate is
quite straightforward using a _Young tower_ construction, shown explicitly for
the Cat Map in [chernov_decay_2000], based on the theory developed in
[young_statistical_1998], [young_recurrence_1999]. The theory has been
extended to hyperbolic systems with singularities (chernov_decay_1999)
culminating in schemes from [chernov_billiards_2005], which give conditions
under which systems enjoy exponential or polynomial decay of correlations,
provided hyperbolicity is sufficiently strong. While aimed at application to
billiards systems, the scheme is readily applicable to models of fluid mixing,
most recently in [springham_polynomial_2014] where the mixing rate for a wide
class of linked twist maps is shown to be at worst polynomial.

##  Statement of results

[][] (0,0) rectangle (4,4); (0,0)–(4,3)–(0,4); [¡-¿] (4.1,4)–(4.1,3);
[scale=1.5] at (4.3,3.5) $\eta$; [-¿] (0.4,2.8)–(2,2.8); [scale=1.7] at
(2,4.3) $F$; [][] (0,0) rectangle (4,4); (0,0)–(3,4) – (4,0); [scale=1.5] at
(3.5,4.4) $\xi$; [¡-¿] (3,4.1)–(4,4.1); [-¿] (2.8,0.4)–(2.8,2); [scale=1.7] at
(2,4.3) $G$;

Figure : A family of area preserving maps $H_{(\xi,\eta)}=G_{\xi}\circ
F_{\eta}$ parameterised by $0<\eta,\xi<1$.

Let $H_{(\xi,\eta)}:\tor\rightarrow\tor$ be the composition of two shears
$F_{\eta}$ and $G_{\xi}$, where $F_{\eta}$ is given by (1) and $G_{\xi}$ maps

$(x,y)\mapsto\cases{\left}(x,y+\frac{1}{1-\xi}x\right)\text{mod1}&\text{for}x\leq
1-\xi,\\\ \par\left(x,y+\frac{1}{\xi}(1-x)\right)\text{mod1}\text{for}x\geq
1-\xi.\\\ $

= 0.2; [scale=1.2] at (0,-0.22) [gray!40] (0,4-4*) rectangle (4-4*,4);
[gray!20] (4-4*,0) rectangle (4,4-4*); (0,0) rectangle (4,4); (4-4*,0) –
(4-4*,4); (0,4-4*) – (4,4-4*); [scale=0.8] at (4-4*,-0.2) $1-\eta$;
[scale=1.5] at (2-0.5*4*,2-0.5*4*) $R_{1}$; [scale=1.5] at (4-0.5*4*,2-0.5*4*)
$R_{2}$; [scale=1.5] at (2-0.5*4*,4-0.5*4*) $R_{3}$; [scale=1.5] at
(4-0.5*4*,4-0.5*4*) $R_{4}$; ; [scale=1.2] at (-5.5,0.025) [fill=gray!40] (0,
4*(1-+*)) – (4*(1-),4) – (0,4) – (0, 4*(1-+*)); [fill=gray!40] (4,4) –
(0,4*(1-)) – (4*(1-),4*(1-)) – (4, 4*(1-+*)) – (4,4); [fill=gray!20]
(4*(1-),4*(1-)) – (4,4*(1-)) – (4, 4*(1-)^2) – (4*(1-),4*(1-)); [fill=gray!20]
(0,4*(1-)) – (4,0) – (4*(1-),0) – (0, 4*(1-)^2) – (0,4*(1-)); (0,0) rectangle
(4,4); [scale=0.8] at (-0.4,4-4*) $1-\eta$; [scale=0.8] at (4.2, 4*(1-+*))
$y_{2}$; [scale=0.8] at (4.2, 4*(1-)^2) $y_{1}$; at (2-0.5*4*,2-0.5*4*)
$A_{2}$; at (4-0.3*4*, 4*(1-1.25*) ) $A_{2}$; at (4-0.5*4*,2-0.5*4*) $A_{1}$;
at (4-0.5*4*,4-0.5*4*) $A_{3}$; at (0.5*4*,4-0.3*4*) $A_{3}$; at (2*,2*)
$A_{1}$; ; [scale=1.2] at (5,-0.2) [fill = gray!40] (0,4-4*) – (4**(1-),4) –
(0,4) – (0,4-4*); [fill = gray!40] (0,0) – (4-4*,4) – (4-4*,4-4*) –
(4**(1-),0) – (0,0); [fill=gray!20] (4-4*,0) – (4-4*,4-4*) – (4-4**,0) –
(4-4*,0); [fill=gray!20] (4,0) – (4-4*,4) – (4-4**,4)– (4,4-4*) – (4,0); (0,0)
rectangle (4,4); [scale=0.8] at (4**(1-),-0.2) $x_{1}$; [scale=0.8] at
(4*(1-*),-0.2) $x_{2}$; at (0.6*4*(1-),2) $A_{3}^{\prime}$; at (0.7*4*(1-),1)
$A_{1}^{\prime}$; at (0.4*4*(1-),3) $A_{1}^{\prime}$; at (0.25*4*,4-0.3*4*)
$A_{3}^{\prime}$; at (4-0.3*4*, 4*(1-1.25*) ) $A_{2}^{\prime}$; at (4-2.5*,
2*) $A_{2}^{\prime}$; ;

Figure : A partition of the torus into four rectangles $R_{j}$, and their
images $A_{j}$ under $F_{\eta}^{-1}$, $A_{j}^{\prime}$ under $G_{\xi}$. The
smallest partition elements $A_{4}$ and $A_{4}^{\prime}$ are left unlabelled.
Case illustrated $\xi=\eta=0.2$.

Write $H=H_{(\xi,\eta)}$ and partition the torus into four rectangles $R_{j}$
using the lines $x=0$, $y=0$, $x=1-\xi$, and $y=1-\eta$, as shown in Figure 1.
Letting $A_{j}=F_{\eta}^{-1}(R_{j})$, the derivative $DH$ is defined
everywhere outside of the set $\mathcal{D}=\cup_{j}\partial A_{j}$ and is
constant on each of the $A_{j}$. The matrices $M_{j}=DH|_{A_{j}}$ are given by

$M_{1}=\pmatrix{1}&\frac{1}{1-\eta}\\\
\frac{1}{1-\xi}1+\frac{1}{(1-\xi)(1-\eta)},\quad
M_{2}=\pmatrix{1}&\frac{1}{1-\eta}\\\ -\frac{1}{\xi}1-\frac{1}{\xi(1-\eta)},$
$M_{3}=\pmatrix{1}&-\frac{1}{\eta}\\\
\frac{1}{1-\xi}1-\frac{1}{\eta(1-\xi)},\quad
M_{4}=\pmatrix{1}&-\frac{1}{\eta}\\\ -\frac{1}{\xi}1+\frac{1}{\eta\xi}.$

Letting $A_{j}^{\prime}=G(R_{j})=H(A_{j})$, the derivative of $H^{-1}$ is
defined everywhere outside of the set $\mathcal{D}^{\prime}=\cup_{j}\partial
A_{j}^{\prime}$ and is constant on each of the $A_{j}^{\prime}$. The labelled
intersections with the axes are $y_{1}=(1-\xi)(1-\eta)$,
$y_{2}=1-\eta+\xi\eta$, $x_{1}=\eta(1-\xi)$, and $x_{2}=1-\xi\eta$.

[][][width=0.33]images/exp.png[][][width=0.33]images/islands.png[][][width=0.33]images/otm.png

Figure : Poincaré sections for $H_{(\xi,\eta)}$ at parameter values (a):
$\xi=\eta=1/10$, (b): $\xi=3/10,\eta=6/10$, (c): $\xi=\eta=1/2$. Starting
point is $z_{0}=(1/\sqrt{2},1/\sqrt{3})$, 50,000 iterates shown. In (a) we see
fully ergodic behaviour, in (b) island structures are present, and in (c)
orbits can become trapped near invariant sets of line segments.

The aim of this paper is to establish mixing results on some subset of the
parameter space $0<\xi,\eta<1$. Figure 1 shows the orbit starting at
$z_{0}=(1/\sqrt{2},1/\sqrt{3})$ for three parameter choices, highlighting a
variety of ergodic and elliptic behaviour over the parameter space. A sketch
of the results in the present work is given in Figure 1. In section
LABEL:sec:outline we state two results which give an efficient scheme for
proving the mixing property. Sections LABEL:sec:lyp, LABEL:sec:man,
LABEL:sec:repMan give a simplified proof of the mixing property over a subset
of the one dimensional parameter space $\xi=\eta$. In particular we show: The
map $H=H_{(\eta,\eta)}$ has the Bernoulli property for $0<\eta<\eta_{1}\approx
0.2024$. In section LABEL:sec:rates we prove exponential decay of correlations
for $H$ over this subset. In section LABEL:sec:param we generalise our results
to the wider two dimensional parameter space; making a small technical
adjustment to Lemmas LABEL:lemma:unstableGrowth and LABEL:lemma:stableGrowth,
establishing parameter space symmetries, and proving our main result, Theorem
LABEL:thm:2dMixingResults. Section LABEL:sec:otm looks at the special case
$H_{(\xi,\eta)}$ with $\xi=\eta=\frac{1}{2}$, which exhibits some unique
dynamical features. We conclude with some final remarks in section
LABEL:sec:discussion.

[scale=1.8] at (2.487,2.542) whiteRGB255,255,255 [y=0.80pt, x=0.80pt,
yscale=-0.231*18/25, xscale=0.231*18/25, inner sep=0pt, outer sep=0pt]
[scale=2.000] [fill=white,rounded corners=0.0000cm] (0.0000,0.0000) rectangle
(1000.0000,1000.0000); [draw=black,draw opacity=0.900,line cap=round,line
join=round,line width=0.500pt,miter limit=10.00] (9.3039,989.6961) –
(9.3039,989.6961) – (9.3039,9.3039); [draw=black,draw opacity=0.900,line
cap=round,line join=round,line width=0.500pt,miter limit=10.00]
(9.3039,9.3039) – (9.3039,9.3039) – (989.6961,9.3039); [draw=black,draw
opacity=0.900,line cap=round,line join=round,line width=0.500pt,miter
limit=10.00] (9.3039,989.6961) – (9.3039,989.6961) – (989.6961,989.6961);
[draw=black,draw opacity=0.900,line cap=round,line join=round,line
width=0.500pt,miter limit=10.00] (989.6961,989.6961) – (989.6961,989.6961) –
(989.6961,9.3039); [fill=gray!10,draw=black,draw opacity=0.900,line
cap=round,line join=round,line width=0.500pt,miter limit=10.00]
(9.3039,989.6961) – (9.3039,744.5980) – (9.3039,744.5980) – (45.3984,735.2295)
– (79.8223,725.6021) – (112.5371,715.7524) – (143.5430,705.7141) –
(173.3281,695.3531) – (201.6484,684.7752) – (228.7480,673.9160) –
(254.6270,662.7987) – (279.5293,651.3358) – (303.4551,639.5389) –
(326.6484,627.2905) – (348.8652,614.7263) – (370.3496,601.7191) –
(391.1016,588.2688) – (411.3652,574.2035) – (430.8965,559.6811) –
(449.6953,544.7112) – (467.5176,529.5236) – (484.3633,514.1833) –
(500.2324,498.7665) – (515.3691,483.0999) – (529.7734,467.2339) –
(543.6895,450.9323) – (557.1172,434.2085) – (570.0566,417.0802) –
(582.5078,399.5705) – (594.4707,381.7083) – (605.9453,363.5287) –
(617.1758,344.6515) – (627.9180,325.4984) – (638.4160,305.6487) –
(648.6699,285.0809) – (658.6797,263.7737) – (668.4453,241.7066) –
(677.9668,218.8599) – (687.2441,195.2155) – (696.2773,170.7566) –
(705.0664,145.4686) – (713.6113,119.3394) – (722.1562,91.5412) –
(730.4570,62.7825) – (738.5137,33.0519) – (744.5980,9.3039) –
(989.6961,9.3039); [fill=gray!10,draw=black,draw opacity=0.900,line
cap=round,line join=round,line width=0.500pt,miter limit=10.00]
(9.3039,989.6961) – (254.4020,989.6961) – (254.4020,989.6961) –
(262.4615,958.4844) – (270.7092,928.5365) – (279.1173,899.8906) –
(287.7688,872.2214) – (296.5457,845.8542) – (305.5369,820.4635) –
(314.7295,796.0495) – (324.1076,772.6120) – (333.7944,749.8255) –
(343.6413,728.0156) – (353.7828,706.8568) – (364.2171,686.3490) –
(374.9407,666.4922) – (385.9478,647.2865) – (397.2300,628.7318) –
(408.7761,610.8281) – (420.8012,593.2500) – (433.3289,575.9974) –
(446.1258,559.3958) – (459.4445,543.1198) – (473.3091,527.1693) –
(487.7446,511.5443) – (502.7770,496.2448) – (518.7846,480.9453) –
(535.4903,465.9714) – (552.9281,451.3229) – (571.5612,436.6745) –
(591.0581,422.3516) – (611.4642,408.3542) – (632.3005,395.0078) –
(654.0668,381.9870) – (676.2019,369.6172) – (699.2654,357.5729) –
(723.2900,345.8542) – (748.3075,334.4609) – (774.3481,323.3932) –
(801.4400,312.6510) – (829.6085,302.2344) – (858.8754,292.1432) –
(890.3078,282.0521) – (923.0231,272.2865) – (957.0435,262.8464) –
(989.6961,254.4020) – (989.6961,9.3039); [fill=gray!40,draw=black,draw
opacity=0.900,line cap=round,line join=round,line width=0.500pt,miter
limit=10.00] (221.8336,989.4902) – (221.8336,989.4902) – (230.6817,908.9433) –
(235.8281,857.9290) – (239.5646,815.9062) – (241.9118,784.1282) –
(243.4491,756.5801) – (242.4199,755.5509) – (242.4199,755.5509) –
(214.8717,757.0882) – (183.0938,759.4355) – (141.0710,763.1719) –
(90.0567,768.3183) – (9.5098,777.1664) – (9.5098,989.4902);
[fill=gray!40,draw=black,draw opacity=0.900,line cap=round,line
join=round,line width=0.500pt,miter limit=10.00] (777.1664,9.5098) –
(777.1664,9.5098) – (769.6814,77.6250) – (764.2275,131.4150) –
(760.6707,170.5180) – (757.4084,212.4865) – (755.5602,242.4199) –
(756.5801,243.4491) – (756.5801,243.4491) – (784.1283,241.9118) –
(803.0731,240.5731) – (803.0731,240.5731) – (835.5474,237.8911) –
(886.2188,233.0369) – (956.5313,225.5233) – (989.4902,221.8336) –
(989.6961,9.3039); [fill=gray!80,draw=black,draw opacity=0.900,line
cap=round,line join=round,line width=0.500pt,miter
limit=10.00](9.5098,989.4902) – (221.8336,989.4902) – (221.8336,989.4902) –
(230.6817,908.9433) – (235.8281,857.9290) – (238.0254,831.8217) –
(238.0254,831.8217) – (235.4516,824.6953) – (232.8931,818.8413) –
(231.4307,815.9062) – (225.6891,806.4422) – (223.8117,803.8898) –
(216.4890,795.6946) – (211.6582,791.4400) – (209.2911,789.5762) –
(209.2911,789.5762) – (205.0664,784.4073) – (200.6719,779.9106) –
(191.3599,772.4838) – (190.5775,771.9609) – (179.0758,765.5994) –
(170.8489,762.2332) – (167.1783,760.9746) – (141.0710,763.1719) –
(90.0567,768.3183) – (9.5098,777.1664) – (9.5098,989.4902);
[fill=gray!80,draw=black,draw opacity=0.900,line cap=round,line
join=round,line width=0.500pt,miter limit=10.00](989.6961,9.3039) –
(777.1664,9.5098) – (777.1664,9.5098) – (769.6814,77.6250) –
(764.2275,131.4150) – (760.9746,167.1783) – (760.9746,167.1783) –
(761.1063,167.5812) – (761.1063,167.5812) – (763.5484,174.3047) –
(766.1069,180.1587) – (767.5693,183.0937) – (773.3109,192.5578) –
(775.1883,195.1102) – (782.5110,203.3054) – (787.3418,207.5600) –
(789.4949,209.2305) – (789.7089,209.4238) – (789.7089,209.4238) –
(790.3502,210.2721) – (790.3502,210.2721) – (795.2782,216.0527) –
(798.3281,219.0894) – (807.6401,226.5162) – (808.4225,227.0391) –
(819.9242,233.4006) – (828.1511,236.7668) – (831.8217,238.0254) –
(835.5474,237.8911) – (886.2188,233.0369) – (956.5313,225.5233) –
(989.4902,221.8336) – (989.6961,9.3039); [pattern= dots, draw=black,draw
opacity=0.900,line cap=round,line join=round,line width=0.500pt,miter
limit=10.00] (499.5000,989.6961) – (499.5000,499.5000) – (989.6961 ,499.5000)
– (989.6961 ,989.6961) – (499.5000,989.6961); [pattern= dots, draw=black,draw
opacity=0.900,line cap=round,line join=round,line width=0.500pt,miter
limit=10.00] (499.5000,499.5000) – (499.5000, 9.3039) – (9.3039 ,9.3039) –
(9.3039, 499.5000) – (499.5000,499.5000); [pattern= dots, draw=black,draw
opacity=0.900,line cap=round,line join=round,line width=0.500pt,miter
limit=10.00] (126.9711,618.3268) – (126.9711,618.3268) – (51.2362,600.4115) –
(9.3039,590.7454) – (9.3039,662.8987) – (9.3039,662.8987) – (41.2480,651.8920)
– (72.0098,640.5685) – (101.5508,628.9560) – (126.9711,618.3268); [pattern =
dots, draw=black,draw opacity=0.900,line cap=round,line join=round,line
width=0.500pt,miter limit=10.00] (380.6732,872.0289) – (380.6732,872.0289) –
(396.2285,937.6325) – (408.2546,989.6961) – (336.1013,989.6961) –
(336.1013,989.6961) – (345.6937,961.7396) – (355.3931,935.0469) –
(365.2818,909.3307) – (375.4795,884.2656) – (380.6732,872.0289); [pattern =
dots,draw=black,draw opacity=0.900,line cap=round,line join=round,line
width=0.500pt,miter limit=10.00] (989.6961,408.2546) – (989.6961,408.2546) –
(903.3901,388.1719) – (872.0289,380.6732) – (872.0289,380.6732) –
(872.0289,380.6732) – (900.1348,368.9565) – (929.4316,357.5051) –
(959.9492,346.3272) – (989.6961,336.1013); [pattern = dots, draw=black,draw
opacity=0.900,line cap=round,line join=round,line width=0.500pt,miter
limit=10.00] (590.7454,9.3039) – (590.7454,9.3039) – (608.1426,84.2567) –
(618.3268,126.9711) – (618.3268,126.9711) – (618.3268,126.9711) –
(628.4927,102.6901) – (638.5444,77.2995) – (648.2889,51.2578) –
(657.8431,24.2396) – (662.8987,9.3039); [pattern=dots, draw=black,draw
opacity=0.900,line cap=round,line join=round,line width=0.500pt,miter
limit=10.00] (224.2290,499.5000) – (224.2290,499.5000) – (231.6777,535.7173) –
(237.1550,563.9524) – (237.1550,563.9524) – (237.1550,563.9524) –
(260.7305,550.1834) – (283.4355,536.0539) – (305.6523,521.3202) –
(327.1367,506.1325) – (336.1013,499.5000); [pattern=dots, draw=black,draw
opacity=0.900,line cap=round,line join=round,line width=0.500pt,miter
limit=10.00] (499.5000,774.7710) – (499.5000,774.7710) – (457.1118,766.1016) –
(435.0476,761.8450) – (435.0476,761.8450) – (435.0476,761.8450) –
(446.7654,741.6875) – (458.9566,721.8307) – (471.4125,702.6250) –
(484.3487,683.7448) – (497.7874,665.1901) – (499.5000,662.8987); [pattern =
dots, draw=black,draw opacity=0.900,line cap=round,line join=round,line
width=0.500pt,miter limit=10.00] (761.8450,435.0476) – (761.8450,435.0476) –
(769.2754,473.0152) – (774.7710,499.5000) – (662.8987,499.5000) –
(662.8987,499.5000) – (684.3145,483.9478) – (706.2871,468.9844) –
(728.7480,454.6333) – (751.9414,440.7252) – (761.8450,435.0476); [pattern =
dots, draw=black,draw opacity=0.900,line cap=round,line join=round,line
width=0.500pt,miter limit=10.00] (563.9524,237.1550) – (563.9524,237.1550) –
(520.7843,228.6667) – (499.5000,224.2290) – (499.5000,336.1013) –
(499.5000,336.1013) – (512.8100,317.8594) – (525.8453,298.9792) –
(538.3940,279.7734) – (550.6736,259.9167) – (562.6735,239.4089) –
(563.9524,237.1550); ; [-¿] (0,0) – (0,2.3); [-¿] (0,0) – (2.3,0); at
(5.1,-0.1) 1; at (-0.075,5.1) 1; at (-0.075,-0.075) 0; [scale=1.2] at
(2.3,-0.15) $\xi$; [scale=1.2] at (-0.1,2.3) $\eta$; [fill=white] (0.2,1.89)
circle (3pt); at (0.2,1.89) $\mathcal{I}_{2}$; [fill=white] (1.3,2.4) circle
(3pt); at (1.3,2.4) $\mathcal{I}_{3}$; [fill=white] (1.25,3.75) circle
(4.5pt); [scale=1.5] at (1.25,3.75) $\mathcal{I}_{1}$; [pattern = dots]
(5.5,3.3) rectangle (5.7,3.5); at (6.43,3.38) Elliptic islands; [fill=gray!80]
(5.5,2.7) rectangle (5.7,2.9); at (6.34,2.78) Exp. mixing; at (5.6,2.5)
$\subset$; [fill = gray!40] (5.5,2.1) rectangle (5.7,2.3); at (6.12,2.18)
Mixing; at (5.6,1.9) $\subset$; [fill = gray!10] (5.5,1.5) rectangle
(5.7,1.7); at (6.28,1.58) Hyperbolic;

Figure : Proven behaviour of $H_{(\xi,\eta)}$ over the parameter space
$0<\xi,\eta<1$. Regions $\mathcal{I}_{n}$ and their reflections exhibit
persistent elliptic island structures, explored in section LABEL:sec:islands.