,

# Learning Linear Symmetries in Data Using Moment Matching

Colin Hagemeyer

###### Abstract

It is common in machine learning and statistics to use symmetries derived from
expert knowledge to simplify problems or improve performance, using methods
like data augmentation or penalties. In this paper we consider the
unsupervised and semi-supervised problems of learning such symmetries in a
distribution directly from data in a model-free fashion. We show that in the
worst case this problem is as difficult as the graph automorphism problem.
However, if we restrict to the case where the covariance matrix has unique
eigenvalues, then the eigenvectors will also be eigenvectors of the symmetry
transformation. If we further restrict to finding orthogonal symmetries, then
the eigenvalues will be either be $1$ or $-1$, and the problem reduces to
determining which eigenvectors are which. We develop and compare theoretically
and empirically the effectiveness of different methods of selecting which
eigenvectors should have eigenvalue $-1$ in the symmetry transformation, and
discuss how to extend this approach to non-orthogonal cases where we have
labels

##  1 Introduction

Symmetries exist throughout mathematics and the sciences. They sometimes
appear explicitly in objects or data points, but even more commonly in the
laws and distributions that produce those objects or data points. For example,
in classical mechanics, the laws of motion are invariant under rotations and
translations of space. It doesn’t matter which direction you call “up” and
“down:” the laws of physics remain the same. These symmetries can be exploited
in various ways to make solving problems easier, such as producing
conservation laws in physics.

Much research has been done in using these symmetries in machine learning.
Previous related research has tended to focus on one or more specific goals:

  1. 1.

Exploiting known symmetries in the data. Examples include Group-Equivariant
Convolutional Networks [CW16], data augmentation methods such as image
rotations, and even standard Convolution Neural Networks which implicitly use
translational symmetry.

  2. 2.

Implicitly learning (local) continuous symmetries by learning representations,
such as variational autoencoders, and disentangled representation learning.

  3. 3.

Using known symmetries in the underlying space which may not strictly hold for
the data. Using harmonic analysis one can utilize the prior that the desired
distributions will tend to be more symmetric. This approach is taken in
[Dia88][Kon08].

  4. 4.

Finding symmetries in the data points themselves, such as on abstract grids
like [SKH86].

  5. 5.

Learning identity preserving transformations using identity labels such as
with [She+18] or giving the model pairs equivalent data points such as in
[Deh+21].

  6. 6.

Learning symmetries using the associated differential equations as in the
recent preprint [LT21].

In our approach, we will attempt to learn symmetries in data, similar to the
goal 2, but we will look to solve the potentially more difficult problem of
learning discrete symmetries. By this we mean that there is no way of
continuously deforming the identity map into this transformation, in such a
way that each subsequent deformation is still a symmetry. For example, the
horizontal reflection of a face should still look like a face. This is a
discrete symmetry, which we can see because handedness is binary. On the other
hand, rotation of a face is still a face, and we can continuously deform a 90
degree rotation to a 0 degree rotation by just setting the degrees to be
$d(t)=90(1-t)$.

One of the key motivations of this question is that many symmetries used in
data augmentation are discrete. This is especially the case when moving
outside the the realm of images, such as in natural language processing where
there is no way to continuously deform a sentence. In particular, this leads
to the potential for a novel approach to semi-supervised learning, where a
symmetries could be extracted from a large unlabeled dataset, and then these
symmetries used to improve the performance of supervised learning on a smaller
labeled dataset using data augmentation or other methods.

##  2 Basic Definitions

First we start by defining formally what we mean by a symmetry in a
distribution. To start with, by a probability space we will always mean a
complete, Hausdorff space together with a probability measure on the Borel
$\sigma$-algebra. We’ll typically abuse notation and use $X$ to refer to both
the probability space and the underlying topological space. If the reader is
unfamiliar with this terminology, just consider the case of a probability
distribution on $\mathbb{R}^{n}$, and $\mu(Z)=\int_{Z}p(x)dx$ is the
probability of the event $Z$.

######  Definition 1.

Let $X$ be a probability space, and let $f:X\to X$ be a continuous map. We say
$X$ is $f$-invariant if $\mu(f^{-1}Z)=\mu(Z)$ for all measurable $Z\subseteq
X$. We will also say $f$ is a symmetry of $X$.

Typically we’ll assume the distribution has a continuous PDF $p_{\mu}$ and $f$
is one-to-one and volume preserving, which makes this equivalent to
$p_{\mu}(x)=p_{\mu}(f(x))$. We can also easily extended this definition to a
full groups of symmetries:

######  Definition 2.

Let $X$ be a probability space, and $G$ be a group acting on $X$ by continuous
maps. We say $X$ is $G$-invariant if $\mu(g^{-1}\cdot Z)=\mu(Z)$ for all $g\in
G$ and measurable $Z\subseteq X$. We will say $G$ is a symmetry group of $X$.

In this paper, we will focus on the case where $X$ is also a finite
dimensional vector space. In this case, we want to further impose that any
symmetry or symmetry-group action is linear, and in general we will use a more
restrictive definition of linear as defined below.

######  Definition 3.

Let $X$ be probability space which is also a finite dimensional vector space.
A linear symmetry is a symmetry $f$ such that $f$ is a linear transformation.
An affine symmetry is any symmetry of the form $f(x)=Ax+b$ where $A$ is
linear.

The distinction between linear and affine symmetries is important. In the
remainder of this paper, we will assume that any symmetries are linear in this
more specific sense unless otherwise specified.

##  3 Intractability in the Worst Case

Before we start discussing how to solve find symmetries in data, let’s start
by determining how hard it can be. Unfortunately, even if we know the
probability distribution explicitly, the answer is that it is very hard in
general.

######  Proposition 1.

The graph automorphism problem on $n$ vertices can be reduced to the decision
problem of determining whether a non-trivial linear symmetry exists for a
probability measure on $\mathbb{R}^{n}$.

###### Proof.

To prove this, we will encode the graph automorphism problem into the problem
of finding a non-trivial linear symmetry of a vector space. Let $\Gamma$ a
graphs with $n$ vertices. For simplicity, pick an arbitrary label of the
vertices as $\\{i\\}_{i=1}^{n}$. Let $X=\mathbb{R}^{n}$. We define a
probability measure $\mu$ on $X$ encoding the graph.

$\mu=\frac{\sum\limits_{i=1}^{n}\delta_{e_{i}}+\sum\limits_{(i,j)\in
edges(\Gamma)}2*\delta_{e_{i}+e_{j}}}{n+2|edges(\Gamma)|}$. Where if $v$ is a
vector, $\delta_{v}$ is the Dirac delta measure on $v$.

Since the span of points with non-zero measure is all of $X$, we know that any
linear symmetry must be surjective, and hence also injective. This also forces
points with non-zero measure to be sent to points with equal measure.
Therefore, any non-trival linear symmetry will send unit directions to unit
directions, thus corresponds to an permutation of the graph vertices.
Moreover, if $(i,j)\in edges(\Gamma)$ the symmetry must send $e_{i}+e_{j}$ to
$e_{k}+e_{\ell}$, where $(k,\ell)$ is some other edge in $\Gamma$. Therefore,
the above permutation of graph vertices also respects edges, and hence is an
automorphism. ∎

Note: this problem is not unique to discrete measures. One could find a
similar encoding by swapping out the Dirac delta measures with suitable
continuous measures.

##  4 Symmetries in Distributions

To start looking for symmetries in a set of data $\mathcal{D}$, we will want
to find statistics which will reflect the overall symmetries in the data. The
most direct approach would be to use $\mathcal{D}$ to construct an
approximation $\hat{q}$ of a distribution $q$, and then look for symmetries in
$\hat{q}$. Unfortunately, finding a sufficiently precise approximation of a
distribution in large dimensions is often infeasible unless you are in a very
restrictive setting and have a lot of data.

Instead, we will look for statistics which summarize the data to a lesser
extent, but will still possess the symmetry from the general data. The most
straightforward example would be the cumulants. Many linear transformations
can fix the mean, so let’s consider the first two cumulants: the mean and
covariance matrix. For the remainder of this section, we will consider the
problem of identifying symmetries in the full distribution, and then later we
will discuss how this translates to a sample. In particular, we will walk a
fine line of only relying on quantities that can be approximated with a
sample.

###  4.1 Identifying Symmetries in the Covariance Matrix

######  Proposition 2.

If $q$ is a distribution with linear symmetry $A$ then:

  1. 1.

$A(\mu)=\mu$ for $\mu=\mathbb{E}_{q}[X]$

  2. 2.

$A\Sigma A^{T}=\Sigma$ where $\Sigma$ is the covariance matrix of $q$

###### Proof.

Follows from the linearity of expected values, and a straightforward
calculation. ∎

Since the sample mean and covariance will approach the distribution mean and
covariance, these symmetries should hold approximately for the sample
statistics.

We now wish to make the second main restriction in this paper.

######  Assumption 1.

$A$ is a normal, linear symmetry of finite order.

This is actually more restrictive. Since we assume the symmetry has finite
order, the eigenvalues must be roots of unity, and hence the symmetry will be
orthogonal. Since common known symmetries like image flips and rotations are
orthogonal, this is still a large class of symmetries. Moreover, we still
haven’t escaped the setting where the intractibility of proposition 1 applies.

However, the results of this assumption is that the action on the covariance
matrix becomes an action by conjugation, ie $A\Sigma A^{T}=A\Sigma
A^{-1}=\Sigma$. This means that $A$ and $\Sigma$ must commute. Therefore, $A$
must send each eigenspaces of $\Sigma$ onto itself.

######  Proposition 3.

Let assumption 1 be true. If $\Sigma$ exists and has distinct eigenvalues,
then $A$ is an involution.

###### Proof.

By contradiction, assume that $\Sigma$ has unique eigenvalues but $A$ is not
an involution. Thus, $A$ must have complex non-real eigenvalues (in particular
roots of unity for $n>2$). Since $\Sigma$ has unique eigenvalues, the
eigenvectors of $\Sigma$ are also eigenvectors of $A$. Let $v_{\lambda}$ be an
eigenvector of $\Sigma$ which is in a non-real eigenspace $\lambda$ of $A$.
Since $\Sigma$ is symmetric, its eigenvectors are real, so $v_{\lambda}$ is a
real vector, but $Av_{\lambda}=\lambda v_{\lambda}$ is not. However, since $A$
is a real matrix, this is a contradiction. ∎

######  Corollary 1.

Under the premises of proposition 3, any finite symmetry group $G$ of the
distribution $q$ is isomorphic to $(\mathbb{Z}/2\mathbb{Z})^{n}$

###### Proof.

We know all elements of $G$ must be of order 2. But then $e=(ab)^{2}=abab$ for
all $a,b\in G$. Multiplying by $ba$ on both sides, we get $ab=ba$, so the
group is Abelian, and the result then follows from the fundamental theorem of
finitely generated Abelian groups. ∎

By proposition 3, we can relatively easily identify which sets of data may
higher order symmetries by looking at the eigenvalues of the covariance
matrix. Moreover, if the covariance matrix has degenerate eigenvalues, then
learning an involutional symmetry will be significantly harder. So in order to
escape proposition 1, we will make the following assumption:

######  Assumption 2.

The distribution generating $\mathcal{D}$ has finite covariance with distinct
eigenvalues

This now forces the eigenvectors of $\Sigma$ to be eigenvectors of $A$, and we
know the eigenvalues of $A$ must be $\pm 1$, so if we know $\Sigma$ or its
eigenvectors, we immediately get the following:

######  Proposition 4.

Let assumptions 1 and 2 be true. If $\Sigma=V\Lambda V^{T}$ is an orthogonal
eigenvalue decomposition of $\Sigma$, then $A=VDV^{T}$, where $D$ is a
diagonal matrix with diagonal entries $\pm 1$

###### Proof.

See above discussion. ∎

###  4.2 Finding Symmetric Dimensions

Unfortunately, Assumption 2 is insufficient to give us exact results. In fact,
every matrix of the form $VDV^{T}$ will commute with $\Sigma$, so $\Sigma$
doesn’t contain any more information for this purpose. Another way to see this
is a central Gaussian in $d$ dimensions has a linear symmetry group
$(\mathbb{Z}/2\mathbb{Z})^{d}$ which acts by reflecting each principal axis.

Most distributions don’t have as many innate symmetries as a Gaussian, so we
shouldn’t expect all of these symmetries to also be symmetries of the data we
encounter in practice. To see which of these symmetries are real, let’s
consider other statistics starting with the mean.

######  Proposition 5.

Under the assumptions of proposition 4, let $v_{i}$ be the eigenvectors of
$\Sigma$. These form a basis. Let $\mu=a_{1}v_{1}+...+a_{d}v_{d}$ be it’s the
unique decomposition. If $Av_{k}=-v_{k}$ for some $k$, then $a_{k}=0$.

###### Proof.

$A\mu=b_{1}v_{1}+...+b_{d}v_{d}$ where $b_{i}=\pm a_{i}$. Basis decompositions
are unique, so $a_{i}=b_{i}$. Thus, $A(v_{k})=v_{k}$ or $a_{i}=0$. ∎

This is where we needed $A$ to be linear. if $A$ is allowed to be affine with
fixed center, then we would need to first shift the mean before looking for
zeros. If we don’t know the center, then the zeros could be everything or
nothing depending on the offset. Therefore, this approach should only be used
when the $0$ point of our data is meaningful.

In particular, let’s consider a Bayesian approach perspective, and condition
on the covariance. As long as the prior distribution on the mean doesn’t have
any innate discrete concentrations at $0$, we should expect any zeros in the
decomposition of $\mu$ come from a symmetry.

We can take alternative approach related to the median instead of the mean,
and in practice, the corresponding sample approach is more robust.

######  Proposition 6.

Under the assumptions of proposition 4, let $v_{i}$ be the eigenvectors of
$\Sigma$. If $Av_{k}=-v_{k}$ for some $k$, then $median(Proj_{v_{k}}(X))=0$ or
equivalently

$Prob(Proj_{v_{k}}(X)<0)=Prob(Proj_{v_{k}}(X)>0)$.

###### Proof.

$Proj_{v_{k}}(A(a_{1}v_{1}+...+a_{n}v_{n}))=-a_{k}=-Proj_{v_{k}}(a_{1}v_{1}+...+a_{n}v_{n})$,
so the projected distribution is symmetric about $0$ ∎

The mean and median approaches have the advantage of simplicity, but they also
have two main weaknesses. First, it restricts us to finding linear
transformations, and second it may cause a false positive if the mean happens
to be zero in an eigenvector direction for other reasons.

The first disadvantage can be solved by instead looking at measures of
skewness. These will allow us to find affine symmetries. The two simplest
choices would be the 3rd cumulent, and the non-parametric skewness. If the
projected distribution has 0 skew, then that is strong evidence that there is
really a symmetry here, and we can find the offset by looking at either the
mean or median.

However, to be sure that at least the distribution on $Proj_{v_{k}}(X)$ is
symmetric, we can use the distance skewness. If the distance skewness is $0$,
then the projected distribution must be symmetric, thus avoiding both
weaknesses. However, it’s still possible that the symmetry in the projection
doesn’t come from a symmetry in the full distribution, although that would be
quite the coincidence.

######  Proposition 7.

Under the assumptions of proposition 4, let $v_{i}$ be the eigenvectors of
$\Sigma$. If $Av_{k}=-v_{k}$ for some $k$, then the skewness, non-parametric
skewness, and distance skewness (around $0$) of $Proj_{v_{k}}(X)$ are all $0$

###### Proof.

Again follows from the projected distribution being symmetric around $0$ ∎

To confirm that a symmetry is a true symmetry with complete certainty, we
could use a non-degenerate Maximum Mean Discrepency or a KS-test, but the key
issue here is the curse of dimensionality. In particular, let’s define an
unfixed vector $v$ to be an eigenvector of $\Sigma$ such that $\exists g\in G$
with the property that $gv=-v$. The above propositions give approaches to
determine which vectors may be unfixed vectors.

###  4.3 Distributions with Multiple Non-trivial Symmetries

Previously, we identified how symmetries affected a variety of statistics, and
how to use these to infer the unfixed vectors. If there is only a single non-
trivial symmetry in the distribution, then we can simply negate all the
unfixed vectors and combine them into a transformation via $VDV^{T}$. However,
if we don’t know how many symmetries the distribution has, you could just as
easily have a symmetry group which negates each unfixed vectors independently,
or anything in between. In particular, the negation of all unfixed vectors
need not be a symmetry. Therefore, we need a way to distinguish between these
possibilities. Below we provide one possible approach.

Let $G$ be the complete group of linear symmetries of the distribution. The
key is that we have an action on each principal axis $v_{\lambda}$,
$gv_{\lambda}=\pm v_{\lambda}$ for all $g\in G$. In particular, if
$v_{\lambda}$ is an unfixed vector, then $gv_{\lambda}=-v_{\lambda}$ for some
$g$, and so $\\{-v_{\lambda},v_{\lambda}\\}$ is an orbit under $G$. Then by
the orbit-stabilizer theorem, we know that $|stab_{G}(v_{\lambda})|=|G|/2$.
One approach to computing $stab_{G}(v_{\lambda})$ is to use
$stab_{G}(v_{\lambda})=Aut(\\{x|\langle v_{\lambda},x\rangle>0\\})$. In
particular, we can tell that $Aut(\\{x|\langle v_{\lambda},x\rangle<0\\})$ is
trivial if there are no unfixed vectors. Therefore, we can recursively compute
$G$.

To be explicit, let’s first define the set $X_{V}^{+}:=Aut(\\{x|\langle
v,x\rangle>0|\forall v\in V\\})$. The algorithm goes as follows: we have two
variables: $S=\\{v_{\lambda}\\}_{\lambda\in\Lambda}$ which is a set of unfixed
eigenvectors, and a sequence of eigenvectors $F=(v_{\lambda_{i}})_{i}$ keeping
track of which vectors we’ve fixed. First, compute the unfixed vectors of $X$,
and place them in $S$. Second, take the largest $\lambda$ corresponding to an
unfixed vector which we’ll denote by $\lambda_{0}$, and append
$v_{\lambda_{0}}$ to $F$. Then replace $S$ with the unfixed vectors of
$X_{F}^{+}$ which are also unfixed vectors of $X$. Repeat this process until
$S=\emptyset$.

From this, if $n$ is the length of $S$, we can conclude that $|G|=2^{n}$ (and
thus in particular $G\cong(\mathbb{Z}/2\mathbb{Z})^{n})$). Moreover, the last
non-empty set of unfixed vectors, $S_{n-1}$, gives us a true (minimal)
symmetry $A_{0}$ by negating all vectors in $S_{n-1}$ while fixing the rest of
the basis.

######  Lemma 1.

At each step in the above algorithm, $G(X^{+}_{F})=X$ and any symmetry $h$ of
$X^{+}_{F}$ which commutes with $G$ is an element of $G$.

###### Proof.

By induction, the basis is clear since $G(X)=X$, and a symmetry of $X$ is an
element of $G$ by definition. Next, assume that $G(X^{+}_{F_{k}})=X$ and each
symmetry of $X^{+}_{F_{k}}$ is in $G$. By construction $v_{\lambda_{k+1}}$ is
an unfixed vector, so by the inductive hypothesis there exists a $g\in G$ such
that $gv_{\lambda_{k+1}}=-v_{\lambda_{k+1}}$ but also by construction
$gv_{\lambda_{i}}=v_{\lambda_{i}}$ for all $i<k$. Therefore,
$gX^{+}_{F_{k+1}}=X^{+}_{F_{k}}\cap X^{-}_{\\{v_{\lambda_{k+1}}\\}}$, and
similarly for combinations of signs of $F_{k+1}$, and hence
$G(X^{+}_{F_{k+1}})=X$. Now let $h$ be a symmetry of $X^{+}_{F_{k+1}}$.
Consider the summing operator $sum_{G}:=\sum_{g\in G}g_{*}$. We get the
following:

$2^{m}\mu=sum_{G}\mu\downharpoonright\
X^{+}_{F_{k+1}}=sum_{G}h_{*}\mu\downharpoonright
X^{+}_{F_{k+1}}=h_{*}sum_{G}\mu\downharpoonright
X^{+}_{F_{k+1}}=2^{m}h_{*}\mu$ (1)

Where $\mu$ is the distribution on $X$, and $2^{m}$ is the order of the
stabilizer of $X^{+}_{F_{k+1}}$ in $G$. Thus $h_{*}$ respects $\mu$, and so by
definition is an element of $G$. ∎

In particular, this implies that $A_{0}$ is an element of $G$, and
$G(X_{F_{n}})=X$. In fact, we can get a full generating set as follows. Let
$F_{-j}=(v_{\lambda_{i}})_{i\neq j}$ be the subsequence of a $F$ where the
$j$th element is omitted. Since $F_{n}\subseteq F_{-j}$, we know that
$G(X_{F_{-j}})=X$.

######  Proposition 8.

The set of unfixed vectors of $X_{F_{-j}}^{+}$ which are unfixed vectors of
$X$ gives a non-trivial element $g_{j}\in G$, by negating each unfixed vector,
fixing each fixed eigenvector of $\Sigma$, and then extending by linearity.
Together these form a generating set of the group of symmetries $G$.

###### Proof.

First let’s prove that the $g_{j}$ are elements of $G$. By insisting that we
only consider unfixed vectors of $X$ we force the $g_{j}$ to commute with $G$.
Then the result follows by the same summing argument as Lemma 1.

Finally we will prove that $g_{j}(v_{\lambda_{j}})=-v_{\lambda_{j}}$, and
since we know $g_{i}(v_{\lambda_{j}})=v_{\lambda_{j}}$ for all $i\neq j$ by
construction, this action induces a surjection onto $(\mathbb{Z}/2/Z)^{n}$,
which must be an isomorphism by the equality of cardinality. By strong
induction, first we look at $g_{n}$. We already know that it negates
$v_{\lambda_{n}}$ by construction. Next, assume that we know
$g_{k+1}v_{\lambda_{i}}=-v_{\lambda_{i}}$ for all $i>k$. We also know by
construction that there exists an element $h_{k}\in Stab(X_{S_{k}}^{+})$ such
that $h_{k}v_{\lambda_{k}}=-v_{\lambda_{k}}$. Unfortunately, it may not fix
each $v_{\lambda_{i}}$ for $i>k$. However, for each $i>k$ such that
$h_{k}v_{\lambda_{i}}=-v_{\lambda_{i}}$, we can multiply $h_{k}$ by $g_{i}$ to
prevent this. Assume $\\{\ell_{k}\\}_{k}$ is the set of all $\ell>k$ such that
$h_{k}v_{\lambda_{\ell}}=-v_{\lambda_{\ell}}$. Letting
$g_{k}^{\prime}=g_{i_{1}}\cdots g_{i_{\ell}}h_{k}$ then produces an element of
$G$ which stabilizes $X_{F_{-j}}^{+}$ while negating $v_{\lambda_{j}}$, which
implies that $v_{\lambda_{j}}$ is unfixed vector of $X_{F_{-j}}^{+}$, and
hence $g_{k}$ negates it by construction. ∎

In practice the main issue with this algorithm is that we have reduced the
size of our space by a factor of $|G|$. Therefore, for large $|G|$, an
alternative way of determining the elements of $G$ which fix $F_{k}$ may be
needed.

###  4.4 Non-Orthogonal Symmetries

The case of learning a symmetry that isn’t orthogonal is a bit more tricky.
The fundamental problem is that the covariance matrix lives in $X\otimes X$,
and not in $End(X)\cong X\otimes X^{*}$, so in reality the action on the
covariance is $A\cdot\Sigma=A\Sigma A^{T}$, and not $\Sigma=A\Sigma A^{-1}$.
The former action doesn’t necessarily respect eigenvectors. On the other hand,
if there is an invariant inner product, we get an isomorphism
$\phi:X\xrightarrow{\sim}X^{*}$ which is respected by the symmetry $A$. From
this we get $(id\otimes\phi)(\Sigma)\in End(X)$, where the action by $A$ is
$A\cdot(id\otimes\phi)(\Sigma)=A(id\otimes\phi)(\Sigma)A^{-1}$, and therefore
the eigenvectors of $(id\otimes\phi)(\Sigma)$ are eigenvectors of $A$. Note
that we could just have easily asserted that $A$ was symplectic or really that
it respects any known non-degenerate form.

If we don’t know of such a form, then $\Sigma$ can most naturally be thought
of as a (positive semi-definite) symmetric form. Without any additional
information, the most we can then say about $A$ is that it is orthogonal with
respect to this form, which means we still have a continuous space of
possibilities.

Luckily, there is a common setting where this issue can be resolved. Assume we
have supervised data, with a non-trivial discrete labeling, for example
labeled images. Assume moreover the labeling is invariant under the action of
$A$. Then we know that $A$ must respect the covariance matrix of each of the
labeled subdistributions $p(X|y=c)$ separately. Thus, if we assume one of the
subdistributions has a non-degenerate covariance matrix, we can rewrite the
space $X$ in terms of its eigenvector basis (or the eigenvectors for the
combined dataset), and then use the previous methods on the other
subdistributions. Of course, all of this assumes the existence of a symmetry
fixing multiple subdistributions simultaneously, which is going to be rarer
than a symmetry fixing a single subdistribution.

###  4.5 Labeled Distributions and Subrepresentations

If we do have access to some invariant labels in a dataset but also know that
the symmetry is orthogonal, then the problem becomes somewhat easier. Now, we
have a set of covariance matrices $\Sigma_{c}$ for each label $c$, and the
symmetry transformation $A$ must commute with all of these matrices. In
mathematical terms, $A$ must be an intertwiner for the defining representation
of the algebra $\langle\Sigma_{c}\rangle_{c}$. If $A$ is an involution, such
as is the case if we take assumption 2 that $\Sigma$ have distinct
eigenvalues, this implies that the $+1$ and $-1$ eigenspaces of $A$ must be
subrepresentations, and in particular that $\langle\Sigma_{c}\rangle_{c}$ is
decomposable.

In more concrete terms, there is a basis $\\{v_{i}\\}_{i}$ such that the set
of matrices $\\{\Sigma_{c}\\}_{c}$ is simultaneously block diagonal. These
subspaces then take the place of the eigenvectors in the non-labeled case. We
just need to select some number of these subspaces to be negated in order to
get a matrix which respects the covariances. Generically, we should expect
that as long as the number of labels is sufficiently large (maybe equal to or
greater than $4$ cf. [AKM13]), then any decomposition $\\{\Sigma_{c}\\}_{c}$
should come from the symmetry. In particular, if there is only a single non-
trivial symmetry, then there should be only $2$ subrepresentations to choose
from.

Unfortunately, the noise coming from the sampling requires that our method of
detecting subrepresentations needs to be robust. Without distinct eigenvalues,
we might need to use a more complicated approach as in [MM11], but this is
problematic because it requires finding the eigenvectors of a $d^{4}$
dimensional matrix. As a more direct approach, we can diagonalize the total
covariance matrix $\Sigma$ (including the non-labeled data), then rewrite each
$\Sigma_{c}$ in the corresponding basis, and then in this basis each
$\Sigma_{c}$ should be almost block diagonal. In order to control for
different variances, we’ll use the correlation instead. We can write
$\tilde{P}=\frac{1}{C}\sum_{c}abs(P_{c})$ where $abs(P_{c})$ is the pointwise
absolute value of the correlation matrix associated to $c$ in the basis which
diagonalizes $\Sigma$ (alternatively, we can make this distribution closer to
a normal distribution by taking the average of
$\tilde{\theta}_{c}=\frac{1}{2}(\log(1+abs(P_{c}))-\log(1-abs(P_{c}))$
instead). To find the block diagonal structure, we can interpret this $d\times
d$ matrix as a similarity matrix, and then do clustering in order to determine
the blocks. Moreover, we can then improve the clustering and determine which
eigenvectors are fixed by adding in an auxiliary point $p_{aux}$ whose
similarity to each eigenvector is equal to the average of normalized absolute
means:
$similarity(p_{aux},p_{\lambda})=\frac{1}{C}\sum_{c}\frac{\lvert\langle\mu_{c},v_{\lambda}\rangle\rvert}{\sqrt{\lambda}}$.

Note: the best (similarity-based) clustering algorithm will depend on the the
size of $d$. For small $d$, we can directly find an optimal normalized cut. In
the case where this becomes intractable, we need to be careful about the
potential for eigenvectors with large error. These eigenvectors will tend to
have a large affinity to $2$ or more clusters, making it hard for certain
algorithms to split those clusters apart. Methods such as $k$-medoids or
affinity propagation may select these erroneous eigenvectors as exemplars.

As a simple example of a fast algorithm which is robust to the above issue, we
can take a modified single-step version of dissimilarity analysis [MS+64] by
first making two clusters with the auxiliary point as one cluster, and all
other points in the other. Then we iteratively add points to the auxiliary-
point cluster such $d_{\lambda}:=d_{avg}(p_{\lambda},\text{auxiliary
cluster})-d_{avg}(p_{\lambda},\text{swap cluster})$ is maximal, until
$d_{\lambda}$ is negative for all $\lambda$, or a desired number of vectors
are selected. If we included the auxiliary node, then this algorithm can be
interpreted as an extension of the unsupervised mean-based ranking method by
incorporating label-specific covariance data.

##  5 Symmetries in Samples

The previous section looked at the case of finding symmetries assuming we know
the true distribution and can compute expected statistics. Now we move to the
more practical setting where we don’t know the true distribution, and only
have a sample which we’ll denote $\mathcal{D}$. However, everything works
more-or-less the same asymptotically except using sample statistics. Assuming
the sample size is large, we know the sample mean and sample covariance will
become arbitrarily close to the distribution mean, median and covariance.
Therefore, assuming distinct eigenvalues of the covariance, the sample
eigenvectors should converge to true eigenvectors. Therefore, as before we can
find a (now approximate) symmetry by negating the eigenvectors such that the
corresponding direction has mean, median, or skew that is close enough to $0$.
If the symmetry is actually a symmetry of the sample, then these statistics
will be $0$ on the nose, but otherwise we need to develop a test for when they
are close enough to $0$.

The key challenge is that we have two points of approximation. First, we’re
approximating the covariance eigenvectors, and second we are approximating
some other statistic which dependent on that eigenvector to determine if it is
likely an unfixed vectors. One approach would be to ignore the uncertainty in
the eigenvectors, and then look for the dimensions where the confidence
interval of the desired statistic overlaps with $0$. Unfortunately, this works
poorly in tests with synthetic data.

###  5.1 Ranking Approaches

A more effective approach is to rank the dimensions by how appropriate it
would be to negate them. When using the mean, it is possible to just use the
magnitude of a eigenvectors coefficient in the decomposition
$\mathbb{E}[X]=a_{1}v_{\lambda_{1}}+...+a_{d}v_{\lambda_{d}}$. The advantage
of this is that quadratic mean of the $a_{i}$’s for negated $v_{\lambda_{i}}$
is exactly the distance that the mean is moved by the linear transformation.
Moreover, as the sample size grows, these $a_{i}$ will converge to $0$.
Unfortunately, this doesn’t account for the fact that some dimensions have
larger covariance, so using the raw $a_{i}$ will tend to select
transformations that don’t change datapoints much. Therefore, to correct for
this, we can divide out the standard deviation, giving
$\\{a_{i}/\sqrt{\lambda_{i}}\\}_{i=1}^{d}$.

Similarly, we can use other statistics like the the (normalized) median, a
sign test statistic, or the skew. In experiments with synthetic (Gumbel
mixture model) data, the mean and median seem to be the most effective. In
fact, a mixture of the two seems to work best. This makes sense because the
skew requires more data to converge, while the sign statistic will tend to be
misled by mixture data.

Once we have an effective ranking, we have reduced the number of models to
check from $2^{d}$ to $d$, which allows us to take a more global approach to
model selection. There are a number of non-parametric approaches to comparing
samples. In this case, we will use maximum mean discrepancy (MMD)[SS02] with
an Squared Exponential kernel $k(x,x^{\prime})=e^{\lVert
x-x^{\prime}\rVert^{2}}$. This has the two main advantages that it is zero
only when the distributions are equal, and it is differentiable. We will use
the differentiability later in order to fine-tune the transformation. As a
variant, we can use a weighted $L^{2}$ norm to greater emphasize errors coming
from lower variance dimensions. In particular, when doing model selection with
MNIST data we will use a regularized matrix of the form
$\Omega^{\prime}=((1-\alpha)\Sigma+\alpha I)^{-1}/h$, giving us
$k_{\Omega^{\prime}}(x,x^{\prime})=e^{-(x-x^{\prime})^{T}\Omega^{\prime}(x-x^{\prime})/2}$.

After fixing a scoring statistic, we then have a couple approaches to applying
it. The most direct approach is to fit a separate transformation $T_{i}$ on
all the data for each choice $i$ of the number of negative eigenvalues. Then
we can pick the $i$ which minimizes $error(X,X*T)$.

The second approach is to use $k$-fold cross-validated scoring. We split the
data into k disjoint folds $\\{X_{i}\\}$, then fit the linear transformation
$T_{i}$ on the union of $k-1$ folds. Finally we can compute
$error(X_{i},X_{i}*T_{i})$, and average these over $i=1,2,...,k$. This has the
advantage of reducing the bias in the scoring due to fitting process, and so
the results should be more generalizable. It has the disadvantage of being
considerably slower, and we don’t necessarily have a guarantee that the model
we finally fit works well for the full dataset.

There are a few important practical details to note:

  1. 1.

The iterative approach to approximating the MMD tends to have a significant
variance. Because of this, in practice it is important to use repeated
k-folds.

  2. 2.

If the ordering of eigenvectors is incorrect, this approach has a tendency to
select only a single negative eigenvector. In order to select a decent (but
necessarily imperfect) transformation in this case, we can use the one
standard error rule to select the most negative eigenvectors whose error is
less than a standard error above the minimum error.

  3. 3.

If there are more than one non-trivial symmetry, some certain combinations of
unfixed vectors may be incompatible. The simplest example is

$G=\\{I_{3},\>diag(-1,-1,1),\>diag(-1,1,-1),\>diag(1,-1,-1)\\}$

where $diag(-1,-1,-1)$ is not a symmetry. In this case, we may be unable to
get a complete symmetry using this method, so alternative methods would be
necessary.

####  5.1.1 CLT-Based Bayesian Model Selection

A much faster, but biased approach is to use (partial) Bayesian Model
selection on the mean of each dimension. We’ll see that we obtain similar
results. First, as before, we will ignore the uncertainty of the covariance
and assume that the covariance matrix is equal to the sample covariance
matrix, and hence so are the eigenvectors and eigenvalues. For each
eigenvector $v_{\lambda_{i}}$ we define $\mathcal{D}_{i}$ to be the data
projected onto $v_{\lambda_{i}}$, ie
$\\{Proj_{v_{\lambda_{i}}}(x)|x\in\mathcal{D}\\}$, or if $X$ is the design
matrix then $\mathcal{D}_{i}$ has a $1D$ design matrix
$X_{i}=Xv_{\lambda_{i}}$. For each $i$, let $m_{0}$ be a model where
$\mathcal{D}_{i}$ has an unknown distribution with variance
$\lambda_{i}=\hat{\lambda}_{i}$, which is our model if $v_{\lambda_{i}}$ is an
unfixed vector. Let $m_{1}$ be the model that $\mathcal{D}_{i}$ be an unknown
distribution with mean $\mu_{i}$ and variance $\lambda_{i}=\hat{\lambda}_{i}$,
where $\mu_{i}\sim\mathcal{N}(0,\phi(\lambda_{i}))$, for some function $\phi$.
This is our model for a fixed vector.

Next we consider the case where we only know the first two moments of the
data. By the Central Limit Theorem we can approximate the two top-level
distributions of $\hat{\mu}_{i}$ as a Gaussian with variance
$\hat{\lambda}_{i}/N$. Thus we get the following approximations:

$\displaystyle\hat{\mu}|\Sigma,m_{0}$
$\displaystyle\sim\mathcal{N}(0,\hat{\lambda}_{i}/N)$ (2)
$\displaystyle\hat{\mu}|\Sigma,m_{1}$
$\displaystyle\sim\mathcal{N}(\mu_{i}|0,\varphi(\lambda_{i}))+\mathcal{N}(\varepsilon|0,\hat{\lambda}_{i}/N)=\mathcal{N}(0,\varphi(\lambda_{i})+\hat{\lambda}_{i}/N)$
(3)

Which we’ll denote as $P_{CLT}$. We want to select the model which makes the
observed sample mean $\overline{x}$ most likely.

######  Proposition 9.

If $\phi(\lambda)=\lambda$ and $\lambda_{i}\neq 0$, then
$P_{CLT}(\hat{\mu}|\Sigma,m_{0})>P_{CLT}(\hat{\mu}|\Sigma,m_{1})$ if and only
if

$\frac{\overline{x}_{i}}{\sqrt{\hat{\lambda}_{i}}}>\frac{\sqrt{(N+1)\ln(N+1)}}{N}$
(4)

###### Proof.

$\displaystyle\mathcal{N}(x|0,\hat{\lambda}_{i}/N)=\frac{1}{\sqrt{2\pi\hat{\lambda}_{i}/N}}e^{-\frac{x^{2}}{2\hat{\lambda}_{i}/N}}$
$\displaystyle<\frac{1}{\sqrt{2\pi\hat{\lambda}_{i}\frac{N+1}{N}}}e^{-\frac{x^{2}}{2\hat{\lambda}_{i}\frac{N+1}{N}}}=\mathcal{N}(x|0,\hat{\lambda}_{i}+\hat{\lambda}_{i}/N)$
$\displaystyle\Longleftrightarrow$
$\displaystyle\sqrt{N+1}e^{-\frac{x^{2}}{2\hat{\lambda}_{i}/N}}$
$\displaystyle<e^{-\frac{x^{2}}{2\hat{\lambda}_{i}\frac{N+1}{N}}}$
$\displaystyle\Longleftrightarrow$
$\displaystyle\frac{1}{2}\ln(N+1)-\frac{x^{2}}{2\hat{\lambda}_{i}/N}$
$\displaystyle<-\frac{x^{2}}{2\hat{\lambda}_{i}\frac{N+1}{N}}$
$\displaystyle\Longleftrightarrow$ $\displaystyle\ln(N+1)$
$\displaystyle<\frac{x^{2}}{\hat{\lambda}_{i}}(N-\frac{N}{N+1})$
$\displaystyle\Longleftrightarrow$
$\displaystyle\frac{x}{\sqrt{\hat{\lambda}_{i}}}$
$\displaystyle>\frac{\sqrt{(N+1)\ln(N+1)}}{N}$

∎

Therefore, if we use the prior that the two models are equally common, then we
obtain a criterion for determining the cut-off point for the mean-based
ranking approach. Notice that we’ve made a significant assumption that the
mean’s prior distribution has variance equal to the observed variance. This
was chosen primarily to make the math nicer, but asymptotically the results
should converge to the same result.

###  5.2 Combined Bootstrapping Approaches

Another principled, non-parametric approach would be to use end-to-end
bootstrapping. Take the null hypothesis to be that the dimension is symmetric,
and fix a particular univariate statistic $T$, distance function, and
significance $\alpha$. We’ll record data in sets, $S_{v}$ and $N_{v}$, for
each eigenvector of the original covariance matrix $\hat{\Sigma}$. First,
resample the data $m$ times, giving new samples
$\\{\mathcal{D}_{i}\\}_{i=1}^{m}$. Second, compute the covariance matrices
$\hat{\Sigma}_{i}$ for each sample. Third, for each eigenvector $v$ of the
original covariance $\hat{\Sigma}$, determine the closest eigenvector
$v^{\prime}$ of $\hat{\Sigma}_{i}$. Record that distance in $N_{v}$. Finally,
compute $T$ on the projection of $\mathcal{D}$ onto $\mathbb{R}v^{\prime}$,
and record that value in $S_{v}$.

We can use the distances $N_{v}$ to determine if the eigenvalue for $v$ is
likely distinct by comparing them to a uniform distribution on a circle, but
for simplicity, let’s assume that we know that the true eigenvalues of
$\Sigma$ are distinct, and the distances are consistently small, so we know
that we have a consistent correspondence. Then we can order the values in
$S_{v}$, and use the bootstrap percentile interval (under the null hypothesis
this is equivalent to the bootstrap pivotal interval). If the confidence
interval contains $0$, then we fail to reject the null hypothesis, and we add
$v$ to the potential unfixed vector. Otherwise, we conclude that $v$ is a
fixed vector under the action of $G$.

Unfortunately there is significant variance in the eigenvectors even for
relatively large datasets, which can make the closest eigenvector $v^{\prime}$
a bad approximation, and leading to an inability to reject the null
hypothesis, at least without a very large dataset.

###  5.3 Theoretical Bounds

With some fairly direct calculations, we can bound the error of these methods
in terms of the error of the covariance matrix and mean, which will give us
consistency. For now $\hat{\mu}$ and $\hat{\Sigma}$ can be any estimators of
$\mu$ and $\Sigma$, but we’ll be primarily interested in the sample mean and
covariance matrix. Let $\\{\lambda_{i}\\}_{i}$ be the eigenvalues of the
population covariance $\Sigma$, with (some choice of) corresponding unit
eigenvectors $\\{v_{\lambda_{i}}\\}$,
$\mu_{i}=\langle\mu,v_{\lambda_{i}}\rangle$, and
$\tilde{\mu}_{i}=\langle\hat{\mu},v_{\lambda_{i}}\rangle$. Note:
$\tilde{\mu}_{i}$ is not a statistic since it depends on $v_{\lambda_{i}}$,
but importantly its difference with $\mu_{i}$ can be controlled using the
Central Limit Theorem. First, we assume the following bounds:

$\displaystyle\lvert\tilde{\mu}_{i}-\mu_{i}\rvert$
$\displaystyle<\varepsilon_{1i}\;\forall i$ (5)
$\displaystyle\lVert\hat{\Sigma}-\Sigma\rVert_{op}$
$\displaystyle<\varepsilon_{2}$ (6)
$\displaystyle\lvert\lambda_{i}-\lambda_{j}\rvert$ $\displaystyle\geq
2\delta\;\forall i\neq j$ (7)

For some $\varepsilon_{i}>0$ and $\delta>0$. First, we want to bound the error
for the eigenvalues and eigenvectors. Let $\\{\hat{\lambda}_{k}\\}_{k}$,
$\hat{v}_{k}$ be the set of eigenvalues and corresponding eigenvectors for
$\hat{\Sigma}$. However, a priori these may not be in natural correspondence
with the eigenvalues of the population covariance $\Sigma$. Fix a particular
estimated eigenvalue $\hat{\lambda}_{k}$. Since the eigenvalues are distinct
by assumption, we can decompose $\hat{v}_{k}$ in terms of the population
eigenvectors:

$\hat{v}_{k}=\sum\limits_{i=1}^{d}a_{i}v_{\lambda_{i}}$ (8)

Where we know $\sum\limits_{i=1}^{d}a_{i}^{2}=1$, since by assumption
$\hat{v}_{k}$ is normalized, and $v_{\lambda_{i}}$ form an orthonormal basis.
Then can use the bound on the covariance to bound the difference of
eigenvalues:

$\displaystyle\varepsilon^{2}_{2}$
$\displaystyle>\lVert(\hat{\Sigma}-\Sigma)\hat{v}_{k}\rVert_{2}^{2}$ (9)
$\displaystyle=\lVert\hat{\Sigma}v_{\lambda_{k}}-\Sigma\sum\limits_{i=1}^{d}a_{i}v_{\lambda_{i}}\rVert_{2}^{2}$
(10)
$\displaystyle=\lVert\hat{\lambda}_{k}\hat{v}_{k}-\sum\limits_{i=1}^{d}a_{i}\lambda_{i}v_{\lambda_{i}}\rVert_{2}^{2}$
(11)
$\displaystyle=\lVert\sum\limits_{i=1}^{d}(a_{i}\hat{\lambda}_{k}-a_{i}\lambda_{i})v_{\lambda_{i}}\rVert_{2}^{2}$
(12)
$\displaystyle=\sum\limits_{i=1}^{d}a_{i}^{2}(\hat{\lambda}_{k}-\lambda_{i})^{2}$
(13)

Since the $\lambda_{i}$ have minimum pairwise distance $2\delta$, we know that
at most one of the $\lambda_{i}$ satisfies
$\lvert\hat{\lambda}_{k}-\lambda_{i}\rvert<\delta$. Let’s consider the case
where none of them satisfies this inequality. Then:

$\displaystyle\varepsilon^{2}_{2}>\sum\limits_{i=1}^{d}a_{i}^{2}(\hat{\lambda}_{k}-\lambda_{i})^{2}\geq\delta^{2}\sum\limits_{i=1}^{d}a_{i}^{2}=\delta^{2}$
(14)

And therefore, by contrapositive, as long as the covariance error is small
enough, (ie $\varepsilon_{2}<\delta$), there is exactly one such
$v_{\lambda_{i}}$ near $\hat{v}_{k}$, which produces a correspondence. So from
now on, we’ll assume that this inequality holds, and therefore after
relabeling the eigenvectors we get that
$\lvert\hat{\lambda}_{k}-\lambda_{i}\rvert<\delta$ if and only if $i=k$.

Next, we want to bound the distance between $\hat{v}_{k}$ and
$v_{\lambda_{k}}$. From the above assumption, we get the following breakdown
of (13):

$\displaystyle\varepsilon_{2}^{2}>\sum\limits_{i=1}^{d}a_{i}^{2}(\hat{\lambda}_{k}-\lambda_{i})^{2}$
$\displaystyle=a_{k}^{2}(\hat{\lambda}_{k}-\lambda_{k})^{2}+\sum\limits_{i\neq
k}a_{i}^{2}(\hat{\lambda}_{k}-\lambda_{i})^{2}$ (15) $\displaystyle\leq
a_{k}^{2}(\hat{\lambda}_{k}-\lambda_{k})^{2}+\delta^{2}\sum\limits_{i\neq
k}a_{i}^{2}$ (16)

Which in particular gives us:

$\sum\limits_{i\neq k}a_{i}^{2}<\frac{\varepsilon_{2}^{2}}{\delta^{2}}$ (17)

But using $\sum\limits_{i=1}^{d}a_{i}^{2}=1$, we also get a bound on $a_{k}$:

$\displaystyle a_{k}^{2}$ $\displaystyle=1-\sum\limits_{i\neq k}a_{i}^{2}$
(18) $\displaystyle>1-\frac{\varepsilon_{2}^{2}}{\delta^{2}}$ (19)

So as $\varepsilon_{2}\to 0$, $a_{i}\to\delta_{ik}$ where $\delta_{ik}$ is the
Kronecker delta function, and so $\hat{v}_{k}\to v_{\lambda_{k}}$. In
particular, we have:

$\displaystyle\lVert\hat{v}_{k}-v_{\lambda_{k}}\rVert^{2}_{2}$
$\displaystyle=(1-a_{k})^{2}+\sum\limits_{i\neq k}a_{i}^{2}$ (20)
$\displaystyle<\left(1-\sqrt{1-\frac{\varepsilon_{2}^{2}}{\delta^{2}}}\right)^{2}+\frac{\varepsilon_{2}^{2}}{\delta^{2}}$
(21)
$\displaystyle=2\left(1-\sqrt{1-\frac{\varepsilon_{2}^{2}}{\delta^{2}}}\right)=:g(\varepsilon_{2},\delta)^{2}$
(22)

Using the other half of (16), we are also able to bound the difference of
eigenvalues:

$(\hat{\lambda}_{k}-\lambda_{k})^{2}<\frac{\varepsilon_{2}^{2}}{1-\frac{\varepsilon_{2}^{2}}{\delta^{2}}}\\\
$ (23)

So as $\varepsilon_{2}\to 0$ we also get $\hat{\lambda}_{k}\to\lambda_{k}$.
Now let’s consider more specifically our topic, and use the mean inequality
(5). What we really care about is the estimator
$\hat{\mu}_{k}=\langle\hat{\mu},\hat{v}_{k}\rangle$ which we use to determine
which eigenvectors to negate.

$\displaystyle\lvert\hat{\mu}_{k}-\mu_{k}\rvert$
$\displaystyle=\lvert\sum\limits_{i=1}^{d}a_{i}\langle\hat{\mu},\hat{v}_{i}\rangle-\mu_{k}\rvert$
(24)
$\displaystyle=\lvert\sum\limits_{i=1}^{d}a_{i}\langle\hat{\mu},\hat{v}_{i}-v_{\lambda_{i}}\rangle+\sum\limits_{i=1}^{d}a_{i}\langle\hat{\mu},v_{\lambda_{i}}\rangle-\mu_{k}\rvert$
(25)
$\displaystyle\leq\lVert\vec{a}\rVert_{2}\sum\limits_{i=1}^{d}\lVert\hat{\mu}\rVert_{2}\lVert\hat{v}_{i}-v_{\lambda_{i}}\rVert_{2}+\lvert\sum\limits_{i=1}^{d}a_{i}\tilde{\mu}_{i}-\mu_{k}\rvert$
(26)
$\displaystyle<d\cdot\lVert\hat{\mu}\rVert_{2}g(\varepsilon_{2},\delta)+\lvert\tilde{\mu}_{k}-\mu_{k}\rvert+\lvert(a_{k}-1)\tilde{\mu}_{k}+\sum\limits_{i\neq
k}a_{i}\tilde{\mu}_{i}\rvert$ (27)
$\displaystyle<d\lVert\hat{\mu}\rVert_{2}g(\varepsilon_{2},\delta)+\varepsilon_{1k}+g(\varepsilon_{2},\delta)\sqrt{\sum\limits_{i=1}^{d}\tilde{\mu}_{i}^{2}}$
(28)
$\displaystyle<(d+1)g(\varepsilon_{2},\delta)\lVert\hat{\mu}\rVert_{2}+\varepsilon_{1k}$
(29)

Putting it all together, we get the following theorem:

######  Theorem 1.

Assume inequalities (5) - (7) hold for some estimators $\hat{\mu}$ and
$\hat{\Sigma}$, and further $\varepsilon_{2}<\delta$, then we get the
following error bound:

$\lvert\hat{\mu}_{k}-\mu_{k}\rvert<\varepsilon_{1k}+\sqrt{2}(d+1)\lVert\hat{\mu}\rVert_{2}\sqrt{1-\sqrt{1-\frac{\varepsilon_{2}^{2}}{\delta^{2}}}}$
(30)

In particular, the operator norm is bounded by the Frobenius norm which is
$\sqrt{n}$ convergent, and
$\sqrt{2}\sqrt{1-\sqrt{1-\frac{\varepsilon_{2}^{2}}{\delta^{2}}}}\to\frac{\varepsilon_{2}}{\delta}$
as $\varepsilon_{2}\to 0$.

######  Corollary 2.

Let $\hat{\Sigma}$ and $\hat{\mu}$ be the sample covariance matrix and sample
mean respectively, and assume inequality (7) holds, then $\hat{\mu}_{k}$ is
strongly consistent for $\mu_{k}$ and if further and the fourth moment is
finite, then it is $\sqrt{n}$-consistent.

###### Proof.

Since $\hat{\mu}$ converges a.s. by the (strong) law of large numbers, we know
$\lVert\hat{\mu}\rVert_{2}$ is bounded. The former statement then follows by
the (strong) law of large numbers applied to $\mu$ and $\Sigma$, while the
latter follows from the CLT. ∎

Finally we can use use this to analyze our previous methods:

######  Corollary 3.

Make Assumptions 1 and 2 and let $\hat{\Sigma}$ and $\hat{\mu}$ be the sample
covariance matrix and sample mean respectively. Consider the population
property $T_{i}(P):=\begin{cases}1&\mu_{i}=0\\\ 0&\mu_{i}\neq 0\end{cases}$,
then any estimator of the form
$\hat{T}_{i}(\mathcal{D})=\mathbb{I}(\hat{\mu}_{i}<a_{n})$ where
$\frac{\sqrt{n}}{a_{n}}\to 0$ and $a_{n}\to 0$ as $n\to\infty$ is consistent.

In particular, this shows that method in subsection 5.1.1 is consistent
assuming that the only $\mu_{i}$ which are $0$ come from unfixed vectors.

###  5.4 Covariance-Adjusted Rankings

The theoretical results of subsection 5.3 suggest approach to incorporating
the covariance error via the approximate error bound
$\varepsilon_{1k}+(d+1)\lVert\hat{\mu}\rVert_{2}\frac{\varepsilon_{2}}{\delta}$.
Unfortunately, in the current form, the adjustment is a constant with respect
to $k$.

To obtain a more useful adjustment, we need to make things more local. To do
this, we consider the additional error that comes from each $a_{i}\neq 0$ for
$i\neq k$. Asymptotically, this contributes
$\lvert\hat{\mu}_{i}-\hat{\mu}_{k}\rvert\lvert a_{i}\rvert$ worth of error.
However, considering the actual problem of finding the symmetry, we are not
worried about the case when $\lvert\mu_{i}\rvert\leq\lvert\mu_{k}\rvert=0$,
because then $\lvert\mu_{i}\rvert$ wouldn’t be contributing to an incorrect
decision. Therefore, we replace it with
$(\lvert\hat{\mu}_{i}\rvert-\lvert\hat{\mu}_{k}\rvert)_{+}$, which
asymptotically only adds error when $v_{\lambda_{i}}$ is fixed and
$v_{\lambda_{k}}$ is unfixed. Next, we estimate the sample variance of $a_{i}$
by
$\frac{1}{\lvert\lambda_{k}-\lambda_{i}\rvert}\frac{Corr(\langle\hat{v}_{i},(X-\mu)\rangle^{2},\langle\hat{v}_{k},(X-\mu)\rangle^{2})}{\sqrt{n}}$,
giving a ranking statistic:

$T_{k}(\mathcal{D})=\hat{\mu}_{k}/\left(\sqrt{\lambda_{k}}+\sum\limits_{i\neq
k}\frac{(\lvert\hat{\mu}_{i}\rvert-\lvert\hat{\mu}_{k}\rvert)_{+}}{\lvert\lambda_{i}-\lambda_{k}\rvert}Corr(\langle\hat{v}_{i},(X-\hat{\mu})\rangle^{2},\langle\hat{v}_{k},(X-\hat{\mu})\rangle^{2})\right)$
(31)

In the experiment section we will see that this extra complication can
significantly improve performance.

##  6 Fine-tuning using a Linear MMD Network

We started by forcing the symmetry to respect the second cumulant (ie the
covariance) rather than a more holistic view including the first cumulant (the
mean). This is because the mean has insufficient information to narrow down
our search to finitely many transformations. However, this approach is
somewhat myopic and the eigenvectors found may have some error that could be
reduced by including more information. To address this, we can instead try to
learn a transformation that minimizes the MMD (Maximum Mean Discrepancy) for a
non-degenerate kernel (we’ll use an Squared Exponential Kernel). This is in
some sense a natural generalization of the previous approach since training
using the MMD corresponds to moment matching after embedding into a Hilbert
space. Moreover, unlike the previous moment matching

Assume we have a decomposed approximate transformation $A=VDV^{T}$ using the
previous approach, so that $V$ is the space of eigenvectors of the covariance,
and $D$ is a diagonal matrix with entries $\pm 1$. We consider a
transformation $A_{W}=WDW^{T}$ with matrix parameter $W$, where $D$ is fixed.
Initialize $W$ at $V$, and minimize the MMD between $A_{W}X$ and $X$. Since
the MMD is differentiable, this can be approached using a standard gradient
descent-based methods. We used Stochastic Gradient Descent with a momentum of
$0.5$. We also make a slight modification of the sample MMD where we omit the
contribution of point $x$ with its transformed self $A_{W}x$ in order avoid
biasing the algorithm towards learning an identity map.

######  Remark.

One might wonder why we didn’t just use the MDD network approach to begin with
if it is less myopic. The key issue is that this model’s objective function is
highly multimodal. This seems to correspond to the fact that if you have a $W$
which respects the covariance, then there are $2^{n}$ discrete local optima
corresponding to choosing which eigenvectors to negate. To travel between this
modes, the transformation would need to stop respecting some of the
eigenvectors, which makes these modes rather deep. Likely because of this, in
practice with a random initialization, the algorithm tends to get stuck in
poor local optima.

##  7 Experimental Results

In order to determine which approaches are optimal, and determine the
applicability of these methods we’ll use two types of experimentation. First,
we’ll consider synthetic data where we can compute the ground truth error
under certain conditions. Then we’ll consider the more realistic MNIST
dataset, and see how well these techniques might work in practice.

###  7.1 Synthetic Data

In order to make sure the assumptions of the paper hold, we can repeatedly
build synethetic datasets and use these approaches to test them. In order to
avoid models with additional symmetries, we used a mixture of skew base
distributions, in this case Gumbel distributions. We differentiate each Gumbel
distribution in the mixture by multiplying it by a random invertible matrix,
and shifting its mean by vector produced component-wise by a normal
distribution, truncated to stay within $\pm 2\sigma$. The invertible matrix is
produced as the product of a random strictly lower triangular matrix and upper
triangular matrix, where the diagonal is produced in a more involved way in
order to keep things more stable. We first take $d$ samples of a standard
normal distribution truncated between $0.4$ and $2$, and divided by $4$,
making all values $>0.1$. Then we produce a new array so that the $i$th entry
is the sum of the first $i$ random samples. Finally, we permute these values
randomly, and set the diagonal of the upper triangular matrix to be equal to
this array. This ensures that none of the diagonals are within $0.1$ of each
other.

To produce a sample with symmetry, we first fix a particular linear
transformation so that the results can be compared. In this case, we used the
transformation which swaps pairs of coordinate vectors $e_{0}\leftrightarrow
e_{1}$, $e_{2}\leftrightarrow e_{3}$, etc… Then, we take
$N/(2*num\\_clusters)$ samples from one of the Gumbel distributions, add them
to our dataset, sample another $N/(2*num\\_clusters)$ samples from that same
distribution, but then also multiply them by the fixed symmetry. Then we
repeat this for each cluster. Finally, we divide the whole dataset by the
scalar standard deviation for stability and hyperparamter consistency. In the
following, we will use only $2$ clusters.

####  7.1.1 Dataset Parameters and Ranking Methods

First let’s compare the different ranking methods. We’ll keep the model
selection method fixed by using the true number of swapped eigenvectors, and
also keep the computational cost down by not fine tuning. We measure the error
as the $L^{2}$ norm of the difference of the predicted and ground truth
matrices. Fixing the number of dimensions to be $10$, and using $1,000$
datasets, we get the following table containing the MMD and standard errors.

Samples | Mean | Median | MM Mix | Sign | Skew
---|---|---|---|---|---
2,000 | 0.239 $\pm$ 0.003 | 0.236 $\pm$ 0.003 | 0.227 $\pm$ 0.003 | 0.270 $\pm$ 0.003 | 0.326 $\pm$ 0.003
10,000 | 0.168 $\pm$ 0.004 | 0.167 $\pm$ 0.004 | 0.151 $\pm$ 0.004 | 0.211 $\pm$ 0.004 | 0.247 $\pm$ 0.004
50,000 | 0.111 $\pm$ 0.004 | 0.113 $\pm$ 0.004 | 0.097 $\pm$ 0.004 | 0.157 $\pm$ 0.004 | 0.173 $\pm$ 0.004
250,000 | 0.066 $\pm$ 0.003 | 0.064 $\pm$ 0.003 | 0.052 $\pm$ 0.003 | 0.111 $\pm$ 0.004 | 0.104 $\pm$ 0.004

As we can see, an even mixture of Mean and Median seems to be the most
effective approach of this set, although the difference with the mean and
median isn’t particularly large. Moreover, increasing the sample size by a
factor of $5$ seems to correspond to roughly a linear decrease in error on
order of $0.05$. On the other hand, the sign and skew approaches tend to lag
in error, with the skew starting with larger error and but eventually
overtaking the sign statistic.

However, all of these methods only look at the raw statistics without
considering the error in the covariance eigenvectors. Using the mean, but
adjusting for the error in the covariance eigenvalues we get a significantly
improved results which are unfortunately a bit slower to calculate. We can
compare them to the unadjusted best results below:

Samples | MM Mix | Cov Adjusted
---|---|---
2,000 | 0.227 $\pm$ 0.003 | 0.193 $\pm$ 0.003
10,000 | 0.151 $\pm$ 0.004 | 0.109 $\pm$ 0.003
50,000 | 0.097 $\pm$ 0.004 | 0.066 $\pm$ 0.003
250,000 | 0.052 $\pm$ 0.003 | 0.032 $\pm$ 0.002

Looking just at the case of 50,000 samples and a Mean-Median mix, we can plot
the histogram of errors giving Figure 1. There are two clear modes, one
corresponding to being near the global minimum, while the other seems to
correspond to selecting two of the eigenvectors incorrectly. Around 69% of the
tests avoid an incorrect selection, and therefore should be be near the global
optima.

Figure 1: Error histogram for $N=50,000$ and $d=10$

The effect of increasing the dimension is a bit more subtle. If the dimensions
is low enough, this approach will usually find a transformation near the
global optima, and so has a small error. As the dimension increases, the
probability of choosing all the eigenvectors correctly significantly
decreases, as can be seen in Figure 2.

Figure 2: Error histogram for $N=50,000$ and $d=22$. Compare to Figure 1.

However, even though the chance of being near the global optima decreases, the
error eventually stabilizes as can be seen in the first two rows of the
following table, which records the results of different dimensions and sample
sizes using the covariant adjusted statistic (even clearer results happen for
other statistics). This suggests that the bound in Theorem 5.3 which has $d$
dependence isn’t sharp with respect to $d$ complexity, at least at low
dimensions. For each pair of parameters, we produced $1,000$ synthetic
datasets, and tested the algorithm on each, recording the mean of ground truth
errors $\pm$ the standard errors.

N\d | 6 | 10 | 14 | 18 | 22
---|---|---|---|---|---
2,000 | 0.127 $\pm$ 0.005 | 0.190 $\pm$ 0.003 | 0.212 $\pm$ 0.002 | 0.211 $\pm$ 0.001 | 0.208 $\pm$ 0.001
10,000 | 0.064 $\pm$ 0.003 | 0.113 $\pm$ 0.003 | 0.146 $\pm$ 0.003 | 0.153 $\pm$ 0.002 | 0.157 $\pm$ 0.001
50,000 | 0.033 $\pm$ 0.002 | 0.065 $\pm$ 0.003 | 0.087 $\pm$ 0.002 | 0.097 $\pm$ 0.002 | 0.110 $\pm$ 0.002
250,000 | 0.017 $\pm$ 0.002 | 0.031 $\pm$ 0.002 | 0.045 $\pm$ 0.002 | 0.058 $\pm$ 0.002 | 0.069 $\pm$ 0.002

From this we can see that each $5$-fold increase in the sample size leads to a
reduction of mean error of around $0.04-0.05$ for larger dimensions in this
synthetic dataset model. We can also see that a large dataset in the hundreds
of thousands of sample points may be necessary to produce very low error.
Luckily, this is feasible for unsupervised image datasets, as long as standard
incremental methods are used to compute the desired statistics. In fact, we
can see that as the dimension gets very large, the error continues to
decrease. For each entry in the following table we only did $10$ tests but
otherwise using the same methodology as the previous table, but the standard
errors remained small due to a significant lower variance:

N\d | 100 | 200 | 400
---|---|---|---
10,000 | 0.116 $\pm$ 0.002 | 0.092 $\pm$ 0.001 | 0.068 $\pm$ 0.000
50,000 | 0.095 $\pm$ 0.002 | 0.079 $\pm$ 0.001 | 0.063 $\pm$ 0.000

####  7.1.2 Model Selection

So far, we have only considered the case where we actually know the dimension
of the $\lambda=-1$ eigenspace in the symmetry matrix. This is unlikely to be
the case in practice. Here we experimented with two of the approaches of model
selection laid out in Section 5.1, in particular the (5 times repeated) 5-fold
Cross Validation method and the full-dataset MMD method (we approximate it 5
times using batches of size $1,024$ and taking the mean of the results). For
simplicity, we’ll consider $d=10$ and $N=50,000$. We’ll also use the Mean-
Median mix statistic since it seems to perform the best.

Both of these approaches are quite slow, so we only use $50$ tests. The
following tables gives the average results of these tests for choices of the
bandwidth $h$.

| Mean Error $\pm$ SE | Correct Number of Swaps
---|---|---
Full Dataset ($h=1$) | 0.180 $\pm$ 0.020 | 38%
Full Dataset ($h=3$) | 0.132 $\pm$ 0.021 | 62%
Full Dataset ($h=6$) | 0.189 $\pm$ 0.024 | 46%
Cross Validation ($h=1$) | 0.240 $\pm$ 0.019 | 25%
Cross Validation ($h=3$) | 0.147 $\pm$ 0.022 | 54%
Cross Validation ($h=6$) | 0.153 $\pm$ 0.021 | 50%

The increase in error is relatively small, and the majority of the errors in
the methods occur when the ordering is already incorrect, so it seems that
most of the difficulty comes from picking a good order. The difference in
error between the full dataset and cross validation approaches seems
relatively small, although selected an optimal bandwidth may require more work
since the results seem more sensitive to this.

####  7.1.3 Fine-Tuning

The next important factor to consider is fine-tuning as described in Section
6. In particular, for this experiment we will again assume we know the correct
number of swapped eigenvectors (ie $5$). This process is significantly slower
and more finicky than the previous methods, so we only considered $50$
synthetic datasets instead of $1000$, and used $100$ epochs. In fact, addition
epochs seemed to often improve performance, but we kept it at $100$ to balance
performance and accuracy. For Stochastic Gradient Descent we used a learning
rate of $0.1$, and a momentum of $0.5$, with an additional penalty to push the
change-of-basis matrix to remain orthogonal (although other experiments have
indicated this may not be necessary).

If there was numeric instability in the process and NaN values appeared, we
would reduce the learning rate by a factor of $0.3$ and try again until it was
able to finish a full $100$ epochs without issues. The results are plotted in
Figure 3.

Figure 3: Error for $N=50,000$ and $d=10$ over $50$ tests after $100$ epochs
of fine-tuning

This experiment seems to indicate that if the initial results are close enough
to the global optima, then fine-tuning will have a large effect. In this case,
fine-tuning almost always results in a mean error of roughly $0.01$ even if
the error before fine-tuning was an order of magnitude larger. On the other
hand, if the initial results are near one of the non-global optima, then
little improvement should be expected from fine-tuning. Since image data is
high dimensional, and we’ve seen that some incorrect selections are almost
certain for high dimensions, this suggests that fine-tuning will have limited
utility without additional techniques. In particular, the average error
increased slightly due to a few outlier errors.

###  7.2 MNIST and EMNIST Datasets

As a simple example of a more realistic dataset we will consider the MNIST and
EMNIST datasets of handwritten numbers. The two primary reasons for this
decision were that the dimension is smaller than most image datasets, and the
dataset has some natural symmetries, for example the horizontal and vertical
flips for labels like $0$ or $8$. However, after experimentation it turned out
that these specific symmetries were not strong because handwritten letters
tend to be tilted.

In particular, we will consider a couple variants. First, we’ll consider the
semi-synethetic dataset where we include random horizontal flips. This allows
us to see how well our techniques work with realistic data but with a known
perfect symmetry. Second, we will consider subsets of the dataset with a fixed
label, since these more restrictive sets are more likely to have strong
symmetries. Finally, we will consider the full dataset, and attempt to apply
our approach to produce data augmentation to improve supervised learning.

####  7.2.1 Semi-Synthetic (Horizontal Flips)

To analyze the effectiveness of our techniques we will consider a different
metric from the synthetic datasets. Partly because the error in the real set
appears to more often come from error in the covariance, and partly because
the dataset is degenerate so it has no way of learning the full
transformation. So we will directly compute the fraction of covariant
eigenvectors are close to being eigenvectors of the horizontal flip
transformation, and for those that are, we will determine which fraction are
correctly swapped by our estimator transformation.

For the former, for each eigenvector we will compute the angle between
$v_{\lambda}$ and $T_{flip}v_{\lambda}$, and sort the results into three
buckets. If the angle is less than 60 degree, or greater than 120 degrees,
then it’s close to an eigenvector of $T_{flip}$ with eigenvalue $+1$ or $-1$
respectively. Otherwise, we consider it to be in error.

For the latter, we take the half of (non-trivial) eigenvectors with the
smallest statistic, then take those to be predicted unfixed vectors. We then
compare this choice with where they were placed in the $+1$ and $-1$ buckets
above, and calculate what fraction of vectors which are correctly swapped or
fixed.

To get a better sense of how dimension affects this, we resized the data to
$4\times 4$, $10\times 10$, $16\times 16$, $22\times 22$ and $28\times 28$
pixels. For half of each dataset, we did a horizontal flip, and left the other
half alone. For the covariance eigenvectors, we got the following results:

Table 1: Semi-Synthetic Covariance Eigenvector Accuracy dataset \ dim | 4 | 10 | 16 | 22 | 28
---|---|---|---|---|---
MNIST | 100% | 80% | 67% | 57% | 36%
EMNIST | 100% | 98% | 89% | 79% | 62%

There is a clear reduction in accuracy as the dimension increased. Moreover,
there is a very clear improvement in accuracy from using the larger dataset.
Nevertheless, even for the EMNIST dataset it appears that the accuracy for the
full sized images isn’t sufficient.

We get a similar result for selecting the correct eigenvectors. In low
dimensions, we get very high accuracy, which quickly deteriorates as the
dimension increases.

Table 2: MNIST Semi-Synthetic Eigenvector Selection Accuracy dim \ stat | mean | median | mm mix | sign | corr adj | label-based
---|---|---|---|---|---|---
4 | 88% | 88% | 88% | 88% | 88% | 100%
10 | 71% | 61% | 69% | 84% | 74% | 90%
16 | 68% | 57% | 64% | 76% | 68% | 81%
22 | 64% | 58% | 62% | 66% | 64% | 68%
28 | 63% | 57% | 58% | 62% | 61% | 65%

First, it should be noted that the label-based clustering approach makes use
of the labels, so while it is more accurate, it is also less useful for semi-
supervised learning. Putting that column aside, we get a surprising result. In
the semi-synthetic data the sign statistic generally outperforms the other
statistics, including the correlation adjusted mean. This is nearly the
opposite of the synthetic data where the sign statistic did very poorly, while
the correlation adjusted mean dominated. The difference in the performance of
the sign statistic seems likely to be due to the disconnected nature of the
synthetic data. In that case, the sign statistic may detect that the two sides
are balanced, but not notice that one cluster is farther from $0$ than the
other.

Moreover, the ground truth error keeps increasing, unlike the corresponding
synthetic datasets, suggesting that the synthetic datasets are missing some
important properties of more realistic large dimensional datasets.

For the EMNIST dataset, we get similar results. However, unlike with the
covariance eigenvectors, the improvement due to the larger dataset is
significantly smaller in the larger dimensions.

Table 3: EMNIST Semi-Synthetic Eigenvector Selection Accuracy dim \ stat | mean | median | mm mix | sign | corr adj | label-based
---|---|---|---|---|---|---
4 | 100% | 100% | 100% | 100% | 100% | 100%
10 | 86% | 77% | 83% | 89% | 86% | 92%
16 | 78% | 65% | 74% | 78% | 76% | 83%
22 | 68% | 59% | 66% | 69% | 68% | 69%
28 | 67% | 62% | 65% | 69% | 68% | 73%

Model selection seems to be tricky in this case. Even in the $10\times 10$
case, the dimensions that are nearly trivial can produce serious problems, as
seen in Figure 4. There is a noticeable dip at $50$ as there should be (with
the actual local minimum at $51$), but it is relatively shallow, and the first
couple dimensions have lower error because their variance is so small.

Figure 4: EMNIST Data-based MMD Errors using Sign-Based Ranking.
$\sigma^{2}=5$, $\lambda=0.03$, and using 25% of the data as a validation set
to do model selection

The dip seems to be more clear in lower dimensions, and can disappear in
higher dimensions, suggesting we may need to use alternative methods in order
to perform model selection when the dimension is large.

####  7.2.2 Fixed Label

Next we consider the unmodified (E)MNIST dataset, but restricted to a specific
label where symmetries should be easier to find. The results of the semi-
synthetic test suggest that the dimension of the full $28\times 28=784$
dimensional space is probably too large for accurate selections without an
even larger dataset, or incorporating other priors like spacial continuity.
Therefore, we consider the $10\times 10$ case where it is possible to see
what’s happening visually while still keeping a relatively small dimension.

Secondly, as with most unsupervised learning methods, the model is rarely
perfect, so we may need to add more flexibility to model selection methods. In
particular, it is sometimes the case that very low covariance dimensions
contribute very little to the scoring, and so swapping them is roughly the
same as doing nothing, causing the model selection to erroneously just pick
one of those vectors to negate.

If we use the spherical MMD, and just look at the EMNIST dataset of points
labeled $0$ we get the following graph of MMD versus the number of
eigenvectors negated:

Figure 5: EMNIST MMD Errors using Sign-Based Ranking

There are no large dips, and instead the error is dominated by places where
eigenvectors with large eigenvalues are added. We can address this somewhat by
rescaling the MMD so that dimensions with larger variances are adjusted down
as explained in subsection 5.1.

Figure 6: EMNIST MMD Errors using Sign-Based Ranking. $\sigma^{2}=5$ and
$\lambda=0.03$

This greatly reduces the size of the large jumps, but we still don’t get any
significant dips. Part of this is likely due to the innate bias in the ranking
method. Each time we add a dimension to swap, the transformation necessarily
respects the mean and standard deviation less, so reductions in MMD need to
come from reductions in larger moments. However, especially as the bandwidth
increases, the MMD focuses more on the lower order moments, making it very
difficult for a drop to appear in the graph. Unfortunately, setting the
bandwidth to be small is often infeasible because the scoring becomes too
unstable.

We can graph the effect of the learned transformation on a sample of the
dataset, with the top row being the original images, and each subsequent row
involving more and more change.

Figure 7: Potential transformations of $0$. From top to bottom $0$, $13$,
$23$, $45$, $65$, and $89$ swaps

We see that the first $3$ rows there is very little difference, but that for
the $4$th and $5$th we start to see some changes. In particular, some of
columns (especially the $4$th, $9$th and $10$th) seem to switch between
roughly vertical, and slanted top-left to bottom-right. We also see some
errors begin to appear, in particular in the $8$th column. We see a similar
result for images labeled $1$

Figure 8: Potential transformations of $1$. From top to bottom $0$, $20$,
$40$, $50$, $60$, and $80$ swaps

On the other hand, some labeled classes like $4$ seem to be too complicated
for this method, and produce incoherent transformations. It seems to be again
try to switch from straight to slanted, but in the process picks up a lot of
blurring.

Figure 9: Potential transformations of $4$. From top to bottom $0$, $20$,
$40$, $50$, $60$, and $80$ swaps

####  7.2.3 Full Dataset

Ideally we hope to be able to find a useful symmetry that works for all labels
simultaneously, because this could allow us to use it to do semi-supervised
data augmentation. Unfortunately, in this case even only a couple of swapped
dimensions quickly leads to large errors. This is somewhat unsurprising given
that this process didn’t work with some of the individual labels.

Figure 10: Potential transformations of full digits EMNIST dataset. From top
to bottom $0$, $5$, $10$, $15$, $20$, $40$, $60$ and $80$ swaps

Clearly, additional work needs to be done to make this functional for more
complicated datasets. Since the transformation seems to make some of the
images change the label (for example in column $1$), it seems to be necessary
to force the symmetry to respect each single-label dataset individually.
Moreover, since the algorithm has trouble with some labels, it’s probably also
necessary to allow more flexible transformations to be learned.

##  8 Further Directions

###  8.1 More-Global Scoring

One of the key limitations is that scoring of dimensions is primarily local.
We obtain a ranking of dimensions, then just find a cut-off point. However, if
there are multiple linear symmetries, or if we accidentally incorrectly rank
some of the dimensions, there may be no good symmetry that can be obtained by
selecting the top ranked dimensions. Alternative approaches such as greedily
selected dimensions has the potential to be more robust to changes in the
model assumptions.

###  8.2 Non-Linear Symmetries

The current paper focused on the case of linear symmetries because they are
relatively easy to work with, and because many symmetries in images should be
locally linear. However, in order to achieve the potential of this goal, it
will be necessary to learn symmetries which are not linear or affine, since
few symmetries are likely to be globally linear.

There are a number of ways we could try to leverage this work to the non-
linear setting. If the desired symmetry is close to linear, we could
initialize a more flexible model at the linear transformation and then fine
tune it. Alternatively, we could learn a symmetry within the encoding of a
non-linear generative model.

###  8.3 Higher Dimensional Symmetries

The methods in the current paper had difficulty when the dimension $d$ grew
much larger than $100$. Ideally, we would like to be able work with larger
images of order $d\approx 300\times 300\approx\num{100000}$. Beyond the
accuracy issues, this will start to run into computational issues since we’d
need to compute the covariance which would be
$\num{100000}\times\num{100000}=\num{10000000000}$ dimensional, which starts
to become impractical.

The most direct solution is to use dimension reduction techniques. In
particular, we never used the extra structure arising from translation and
dilation symmetry. Enforcing these strictly would restrict the allowable
linear symmetries to reflections and rotations of the image, but weaker
constraints like penalties on discontinuity or failure to respect rescaling
may greatly improve performance in higher dimensions.

###  8.4 Higher Order Symmetries

In this paper, we restricted our attention to the case where our models had
only order $2$ symmetries, so applying the transformation twice gave back the
identity transformation. This restriction was justified by the fact that we
can identify datasets which might have higher order symmetries by looking for
eigenvalues that are approximately equal. An important extension would be to
develop techniques to better handle cases where some eigenvalues are equal. If
there are too many of these, we will start to fall into the intractable
setting, but dealing with a smaller number of these equalities should be
feasible. One would need to identify which eigenvalues are likely equal,
perhaps determine a good basis for the corresponding eigenspace, then use this
to determine if the eigenspace has higher order or even continuous symmetries.

## References

  * [AKM13] Daishi Aiura, Naonori Kakimura and Kazuo Murota  “On the number of matrices to generate a matrix *-algebra over the real field”  In _Linear Algebra and its Applications_ , 2013
  * [CW16] Taco Cohen and Max Welling  “Group Equivariant Convolutional Networks”  In _Proceedings of The 33rd International Conference on Machine Learning_ , 2016
  * [Deh+21] Nima Dehmamy et al.  “Automatic Symmetry Discovery with Lie Algebra Convolutional Network” (https://proceedings.neurips.cc/paper/2021/file/
148148d62be67e0916a833931bd32b26-Paper.pdf)  In _35th Conference on Neural
Information Processing System_ , 2021

  * [Dia88] Persi Diaconis  “Group Representations in Probability and Statistics”  Institute of Mathematical Statistics, 1988
  * [Kon08] Risi Kondor  “Group Theoretical Methods in Machine Learning”, 2008
  * [LT21] Ziming Liu and Max Tegmark  “Machine-learning hidden symmetries” (arXiv: 1108.4616), 2021
  * [MM11] Takanori Maehara and Kazuo Murota  “Algorithm for Error-Controlled Simultaneous Block-Diagonalization of Matrices”  In _SIAM Journal on Matrix Analysis and Applications_ , 2011
  * [MS+64] P. Macnaughton-Smith, W. T. Williams, M. B. Dale and L. G. Mockett  “Dissimilarity Analysis: a new Technique of Hierarchical Sub-division”  In _Nature_ , 1964
  * [She+18] Yujun Shen et al.  “FaceID-GAN: Learning a Symmetry Three-Player GAN for Identity-Preserving Face Synthesis”  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018
  * [SKH86] Terrence J Sejnowski, Paul K Kienker and Geoffrey E Hinton  “Learning symmetry groups with hidden units: Beyond the perceptron”  In _Physica D: Nonlinear Phenomena_ , 1986
  * [SS02] B. Scho lkopf and A.J. Smola  “Learning with kernels: Support vector machines, regularization, optimization, and beyond”  MIT Press, 2002