# Learning Symmetric Rules with SATNet

Sangho Lim∗
School of Computing
KAIST
Daejeon, South Korea
<EMAIL_ADDRESS>
Eun-Gyeol Oh∗
Graduate School of Information Security
KAIST
Daejeon, South Korea
<EMAIL_ADDRESS>
Hongseok Yang
School of Computing and Kim Jaechul Graduate School of AI, KAIST
Discrete Mathematics Group, Institute for Basic Science (IBS)
Daejeon, South Korea
<EMAIL_ADDRESS>

###### Abstract

SATNet is a differentiable constraint solver with a custom backpropagation
algorithm, which can be used as a layer in a deep-learning system. It is a
promising proposal for bridging deep learning and logical reasoning. In fact,
SATNet has been successfully applied to learn, among others, the rules of a
complex logical puzzle, such as Sudoku, just from input and output pairs where
inputs are given as images. In this paper, we show how to improve the learning
of SATNet by exploiting symmetries in the target rules of a given but unknown
logical puzzle or more generally a logical formula. We present SymSATNet, a
variant of SATNet that translates the given symmetries of the target rules to
a condition on the parameters of SATNet and requires that the parameters
should have a particular parametric form that guarantees the condition. The
requirement dramatically reduces the number of parameters to learn for the
rules with enough symmetries, and makes the parameter learning of SymSATNet
much easier than that of SATNet. We also describe a technique for
automatically discovering symmetries of the target rules from examples. Our
experiments with Sudoku and Rubik’s cube show the substantial improvement of
SymSATNet over the baseline SATNet.

**footnotetext: These authors contributed equally to this work.

##  1 Introduction

Bringing the ability of reasoning to the deep-learning systems has been the
aim of a large amount of recent research efforts [Yang et al., 2017, Evans and
Grefenstette, 2018, Cingillioglu and Russo, 2019, Wang et al., 2019, Topan et
al., 2021]. One notable outcome of these endeavours is SATNet [Wang et al.,
2019], a differentiable constraint solver with an efficient custom
backpropagation algorithm. SATNet can be used as a component of a deep-
learning system and make the system capable of learning and reasoning about
sophisticated logical rules. Its potential has been demonstrated successfully
with the tasks of learning the rules of complex logical puzzles, such as
Sudoku, just from input-output examples where the inputs are given as images.

We show how to improve the rule (or constraint) learning of SATNet, when the
target rules have permutation symmetries. By having symmetries, we mean the
solutions of the rules are closed under the permutations of those symmetries.
For example, in Sudoku, if a completed $9\times 9$ Sudoku board is a solution,
permuting the numbers $1$ to $9$ in the board, the first three rows, or the
last three columns always gives rise to another solution. Thus, these
permutations are symmetries of Sudoku.

Our improvement is a variant of SATNet, called SymSATNet, which abbreviates
symmetry-aware SATNet. SymSATNet assumes that some symmetries of the target
rules are given a priori although the rules themselves are unknown. It then
translates these symmetries into a condition on the parameter matrix
$C\in\mathbb{R}^{n\times n}$ of SATNet (or our minor generalisation), and
requires that the parameters have a particular parametric form that guarantees
the condition. Concretely, the translated condition says that the matrix $C$
regarded as a linear map should be equivariant with respect to the group $G$
determined by the given symmetries, and the requirement is that $C$ should be
a linear combination of elements in a basis for the space of $G$-equivariant
symmetric matrices. The coefficients of this linear combination are the
parameters of our SymSATNet, and their number is often substantially smaller
than that of the parameters of SATNet.111SATNet assumes that $C$ is of the
form $S^{T}S$ for some $S\in\mathbb{R}^{m\times n}$ for $m<n$, so that the
number of parameters is $mn<n^{2}$. But often it is still substantially larger
than the number of parameters of SymSATNet. For Sudoku, the former is $18$,
while the latter is $729^{2}$ or $k\cdot 729$ for some $k\in\mathbb{N}$ at
best. The reduced number of parameters implies that SymSATNet has to tackle an
easier learning problem than SATNet, and has a potential to learn faster and
generalise better than SATNet.

Who provides symmetries for SymSATNet? The default answer is domain experts,
but a better alternative is possible. We present an automatic algorithm to
discover symmetries. Our algorithm is based on empirical observation that
symmetries emerge in the parameter matrix $C$ of SATNet in the early phase of
training, as clusters of similar entries. Our algorithm takes a snapshot of
$C$ at some training epoch of SATNet, and finds a group $G$ such that (i)
specific entries of $C$ share similar values by $G$-equivariance condition,
and (ii) the number of SymSATNet parameters under $G$ is minimised.

We empirically evaluate SymSATNet and our symmetry-discovering algorithm with
Sudoku and a problem related to Rubik’s cube. For both problems, our algorithm
discovered nontrivial symmetries, and SymSATNet with manually specified or
automatically found symmetries outperformed the baseline SATNet in learning
the rules, in terms of both efficiency and generalisation.

#### Related work

There have been multiple studies on discovering symmetries present in
conjunctive normal form (CNF) formulas in order to reduce the search space of
satisfiability (SAT) solvers. Crawford [Crawford, 1992] proved that the
symmetry-detection problem is equivalent to the graph isomorphism problem, and
showed how to reduce the complexity of pigeonhole problems using symmetries.
Crawford et al. [Crawford et al., 1996] proposed symmetry-breaking predicates
(SBPs), and Aloul et al. [Aloul et al., 2003, 2006] developed SBPs with more
efficient constructions. For automatic symmetry detection, Darga et al. [Darga
et al., 2004] presented a method that improves the partition refinement
procedure introduced by McKay [McKay, 1976, 1981], and Darga et al. [Darga et
al., 2008] proposed an algorithm that achieve efficiency by exploiting the
sparsity of symmetries. In contrast to these global and static methods,
Benhamou et al. [Benhamou et al., 2010] and Devriendt et al. [Devriendt et
al., 2017] handled local symmetries that dynamically arise during search. The
use of symmetries also appears in NeuroSAT [Selsam et al., 2019], which learns
how to solve SAT problems from examples. NeuroSAT solves a given SAT formula
by message passing over a graph constructed from the formula, and in so doing,
its learnable solver can exploit symmetries in the formulas. All of these
techniques use symmetries to help solve given formulas, whereas our approach
uses symmetries to help learn such formulas. Another difference is that those
techniques find hard symmetries of given formulas, whereas our approach
discovers soft or approximate symmetries in a given SATNet parameter matrix.
In the context of deep learning, Basu et al. [Basu et al., 2021] and Dehmamy
et al. Dehmamy et al. [2021] described algorithms that find and exploit
symmetries via group decompositions and Lie algebra convolutions. But these
techniques are not designed to find symmetries in logical formulas.

Our work is related to the studies on learning logical rules from examples
using gradients. Yang et al. [Yang et al., 2017] proposed neural logic
programming, an end-to-end differentiable system which learns first-order
logical rules, Evans and Grefenstette [Evans and Grefenstette, 2018] proposed
a differentiable inductive logic programming system which is robust to noise
of training data, and Cingillioglu and Russo [Cingillioglu and Russo, 2019]
introduced an RNN-based model to learn logical reasoning tasks end-to-end.
Want et al. [Wang et al., 2018, 2019] presented SATNet using the mixing
method, and Topan et al. Topan et al. [2021] further improved SATNet by
solving the symbol grounding problem, a key challenge of SATNet. Our work
extends these lines of work by proposing how to discover and exploit
symmetries from examples when learning logical rules with SATNet.

##  2 Background

We review SATNet, the formalisation of symmetries using groups, and
equivariant maps. For a natural number $n$, let $[n]=\\{1,\ldots,n\\}$, and
for a matrix $M$, let $M_{i,j}$ be the $(i,j)$-th entry of $M$.

###  2.1 SATNet

A good starting point for learning about SATNet is to look at its origin, the
mixing method [Wang et al., 2018], which is an efficient algorithm for solving
semidefinite programming problems with diagonal constraints. Let
$n,k\in\mathbb{N}$ and $C$ be a real-valued symmetric matrix in
$\mathbb{R}^{n\times n}$. The mixing method aims at solving the following
optimisation problem:

$\displaystyle\operatorname*{argmin}_{V\in\mathbb{R}^{k\times n}}\,\langle
C,V^{T}V\rangle\quad\text{subject to}\ \|v_{i}\|=1\ \text{for}\ i\in[n]$ (1)

where $v_{i}$ is the $i$-th column of the matrix $V$, and $\|v_{i}\|$ is
$L_{2}$ norm of $v_{i}$. The mixing method solves (1) by coordinate descent,
where each column $v_{i}$ of $V$ is repeatedly updated as follows:
$g_{i}\leftarrow\sum_{\begin{subarray}{c}j\in[n],j\neq
i\end{subarray}}C_{i,j}v_{j}$ and $v_{i}\leftarrow-\frac{g_{i}}{\|g_{i}\|}$.
This always finds a fixed point of the equations. In fact, it is shown that
almost surely this fixed point attains a global optimum of the optimisation
problem.

An example of the above optimisation problem most relevant to us is a
continuous relaxation of MAXSAT. MAXSAT is a problem of finding truth
assignments to $n$ boolean variables $b_{1},\ldots,b_{n}$. It assumes $m$
clauses of those variables, $F_{1},\ldots,F_{m}$, where $F_{\ell}$ is the
disjunction of some variables with or without negation:
$F_{\ell}=b_{i_{1}}\vee\ldots\vee b_{i_{p}}\vee\neg b_{i_{p+1}}\vee\ldots\neg
b_{i_{p+q}}$. Then, MAXSAT asks for a truth assignment on the variables that
maximises the number of true clauses $F_{i}$ under the assignment. The rules
of many problems, including Sudoku, can be expressed as an instance of MAXSAT.

To apply the mixing method to MAXSAT, we introduce relaxed vectors
$v_{1},\ldots,v_{n}\in\mathbb{R}^{k}$ that encode the boolean variables, and
construct the matrix $S\in\mathbb{R}^{m\times n}$ that encodes the $m$ clauses
of MAXSAT: the $(\ell,j)$-th entry of $S$ has $1$ if $F_{\ell}$ contains
$b_{j}$; and $-1$ if $F_{\ell}$ includes $\neg b_{j}$; and $0$ if neither of
these cases holds. Then, the problem in (1) is formed with $C=-S^{T}S$, and
solved by the mixing method.

SATNet is a variant of the mixing method where some of the columns of $V$ are
fixed and the optimisation is over the rest of the columns.222The original
SATNet assumes that $C$ has the form $S^{T}S$ for some $m\times n$ matrix $S$.
We drop this assumption and adjust the forward and backward computations of
SATNet accordingly. The main steps of derivations of the formulas for the
forward and backward computations are from the work on SATNet [Wang et al.,
2019]. Concretely, it assumes that the column indices in $[n]$ are split into
two disjoint sets, $\mathcal{I}$ and $\mathcal{O}$ (i.e.,
$\mathcal{I}\cup\mathcal{O}=[n]$ and $\mathcal{I}\cap\mathcal{O}=\emptyset$).
The inputs of SATNet are the columns $v_{i}$ of $V$ with $i\in\mathcal{I}$,
and the outputs are the rest of the columns (i.e., the $v_{o}$’s with
$o\in\mathcal{O}$). The symmetric matrix $C$ is the parameter of SATNet. Given
the input vectors, SATNet repeatedly executes the coordinate descent updates
on each output column, until it converges.

One important feature of SATNet is that it has a custom algorithm for
backpropagation. Let $V_{\mathcal{I}}$ be the matrix of the input columns to
SATNet, and $V_{\mathcal{O}}$ be that of the output columns computed by SATNet
on the input $V_{\mathcal{I}}$ under the parameter $C$. Assume that $l$ is a
loss of the output $V_{\mathcal{O}}$. In this context, SATNet provides
formulas and algorithms for computing the derivatives $\partial l/\partial
V_{\mathcal{I}}$ and $\partial l/\partial C$.

We recall the formulas for the derivatives. Let
$o_{1}<o_{2}<\ldots<o_{|\mathcal{O}|}$ be the sorting of the indices in
$\mathcal{O}$. Assume that SATNet was run until convergence, so that the
output columns in $V_{\mathcal{O}}$ are the fixed point of the coordinate
descent updates: for all $o\in\mathcal{O}$,
$g_{o}=\sum_{\begin{subarray}{c}j\in[n],j\neq o\end{subarray}}C_{o,j}v_{j}$
and $v_{o}=-\frac{g_{o}}{\|g_{o}\|}$. The formulas for $\partial l/\partial
V_{\mathcal{I}}$ and $\partial l/\partial C$ at
$(V_{\mathcal{I}},V_{\mathcal{O}},C)$ are defined in terms of the next
quantities:

$\displaystyle
C^{\prime},D^{\prime}\in\mathbb{R}^{|\mathcal{O}|\times|\mathcal{O}|},\qquad
P\in\mathbb{R}^{|\mathcal{O}|k\times|\mathcal{O}|k},\qquad
U\in\mathbb{R}^{|\mathcal{O}|k\times 1},\qquad
W\in\mathbb{R}^{|\mathcal{O}|k\times n^{2}}.$

They have the following definitions: for $i,j\in[|\mathcal{O}|]$, $p,q\in[k]$,
and $r,s\in[n]$,

$\displaystyle\begin{aligned} (C^{\prime})_{i,j}&=\begin{cases}0&\text{if
}i\,{=}\,j\\\ C_{o_{i},o_{j}}&\text{if }i\,{\neq}\,j,\end{cases}\\\
(D^{\prime})_{i,j}&=\begin{cases}\|g_{o_{i}}\|&\text{if }i\,{=}\,j\\\
0&\text{if }i\,{\neq}\,j,\end{cases}\\\ P_{\begin{subarray}{c}(i-1)k+p,\\\
(j-1)k+q\end{subarray}}&=\begin{cases}(I_{k}-v_{o_{i}}v_{o_{i}}^{T})_{p,q}&\text{if
}i=j\\\ 0&\text{if }i\neq j,\end{cases}\end{aligned}\qquad\begin{aligned}
&U=\left(P((D^{\prime}+C^{\prime})\otimes
I_{k})\right)^{\dagger}\left(\frac{\partial
l}{\partial\operatorname*{vec}(V_{\mathcal{O}})}\right)^{T},\\\
&W_{\begin{subarray}{c}(i-1)k+p,\\\
(r-1)n+s\end{subarray}}=\begin{cases}0&\text{if }r=o_{i}\text{ and }s=o_{i}\\\
V_{p,s}&\text{if }r=o_{i}\text{ and }s\neq o_{i}\\\ V_{p,s}&\text{if }r\neq
o_{i}\text{ and }s=o_{i}\\\ 0&\text{if }r\neq o_{i}\text{ and }s\neq
o_{i}.\end{cases}\end{aligned}$

Here $\otimes$ is the Kronecker product, $\\_^{\dagger}$ is Moore-Penrose
inverse (also known as pseudo inverse), and
$\operatorname*{vec}(V_{\mathcal{O}})$ is the vector obtained by stacking the
columns of $V_{\mathcal{O}}$. Let
$C_{\mathcal{O},\mathcal{I}}\in\mathbb{R}^{|\mathcal{O}|\times|\mathcal{I}|}$
be obtained by restricting $C$ to the indices $(o,i)$ with $o\in\mathcal{O}$
and $i\in\mathcal{I}$. Then,

$\displaystyle\partial l/\partial\operatorname*{vec}(C)=-U^{T}W,\qquad\partial
l/\partial\operatorname*{vec}(V_{\mathcal{I}})=-U^{T}(C_{\mathcal{O},\mathcal{I}}\otimes
I_{k}).$ (2)

SATNet computes the above derivative formulas efficiently by iterative
algorithms.

###  2.2 Symmetries and equivariant maps

By symmetries on a set $\mathcal{X}$, we mean a group $G$ that acts on
$\mathcal{X}$. The acting here refers to a function $\\_\cdot\\_$ from
$G\times\mathcal{X}$ to $\mathcal{X}$, called _group action_ , such that (i)
$e\cdot x=x$ for the unit $e\in G$ and any $x\in\mathcal{X}$, and (ii)
$(g\circ g^{\prime})\cdot x=g\cdot(g^{\prime}\cdot x)$ for all
$g,g^{\prime}\in G$ and $x\in\mathcal{X}$, where $\\_\circ\\_$ is the group
operator of $G$.

We use symmetries of permutations on a finite set. The set $\mathcal{X}$ in
our case is $\mathbb{R}^{k\times n}$, the space of the matrix $V$ in (1), and
$G$ is a subgroup of the group $\mathcal{S}_{n}$ of all permutations on $[n]$.
The group action $g\cdot V$ is then defined by permuting the columns of $V$ by
$g$: for all $i,j\in[n]$, $(g\cdot V)_{i,j}=V_{i,g^{-1}(j)}$. This group
action can be expressed compactly with the $n\times n$ permutation matrix
$P_{g^{-1}}$ where
$(P_{g^{-1}})_{i,j}=\operatorname*{\mathds{1}}_{\\{i=g^{-1}(j)\\}}$ and
$g\cdot V=VP_{g^{-1}}$ for the indicator function
$\operatorname*{\mathds{1}}$. Throughout the paper, we often equate each
element $g$ of $G$ with its permutation matrix $P_{g}$, and view $G$ itself as
the group of permutation matrices $P_{g}$ for $g\in G$ with the standard
matrix multiplication.

One important reason for considering symmetries is to study maps that preserve
these symmetries, called equivariant maps. Let $G$ be a group that acts on
sets $\mathcal{X}$ and $\mathcal{Y}$.

######  Definition 2.1.

A function $f:\mathcal{X}\to\mathcal{Y}$ is _$G$ -equivariant_ or
_equivariant_ if $f(g\cdot x)=g\cdot(f(x))$ for all $g\in G$ and
$x\in\mathcal{X}$. It is _$G$ -invariant_ or _invariant_ if $f(g\cdot x)=f(x)$
for all $g\in G$ and $x\in\mathcal{X}$.

The forms of equivariant maps have been studied extensively in the work on
equivariant neural networks and group representation theory [Cohen and
Welling, 2016, Maron et al., 2019, Wang et al., 2020]. In particular, when $f$
is linear, various representation theorems for different $G$’s describe the
matrix form of $f$. We use permutation groups defined inductively by the
following three operations.

######  Definition 2.2.

Let $G$ and $H$ be permutation groups on $[p]$ and $[q]$, with each group
element viewed as a $p\times p$ or $q\times q$ permutation matrix. The _direct
sum_ $G\oplus H$, the _direct product_ $G\otimes H$, and the _wreath product_
$H\wr G$ are the following groups of $(p+q)\times(p+q)$ or $pq\times pq$
permutation matrices with matrix multiplication as their composition:

$\displaystyle G\oplus H$ $\displaystyle=\\{g\oplus h:g\in G,\,h\in H\\};$
$\displaystyle G\otimes H$ $\displaystyle=\\{g\otimes h:g\in G,\,h\in H\\};$
$\displaystyle H\wr G$ $\displaystyle=\\{\operatorname*{wr}(\vec{h},g):g\in
G,\vec{h}\in H^{p}\\}.$

Here $g\oplus h$ is the block diagonal matrix with $p\times p$ matrix $g$ as
its upper-left corner and $q\times q$ matrix $h$ as the lower-right corner,
$g\otimes h$ is the Kronecker product of two matrices $g$ and $h$, and
$\operatorname*{wr}(\vec{h},g)$ is the $pq\times pq$ permutation matrix
defined by
$\operatorname*{wr}(\vec{h},g)_{(i-1)q+j,(i^{\prime}-1)q+j^{\prime}}=\operatorname*{\mathds{1}}_{\\{g_{i,i^{\prime}}=(h_{i})_{j,j^{\prime}}=1\\}}$
for all $i,i^{\prime}\in[p]$ and $j,j^{\prime}\in[q]$.

The next theorem specifies the representation of $G$-equivariant linear maps
for an inductively-constructed $G$, by describing a basis of those linear
maps. For a permutation group $G$ on $[m]$, let
$\mathcal{E}(G)=\left\\{M\in\mathbb{R}^{m\times m}:Mg=gM,\,\forall g\in
G\right\\}$, the vector space of $G$-equivariant linear maps, where each $g\in
G$ is regarded as a permutation matrix. See Appendix B for the proof of
Theorem 2.3.

######  Theorem 2.3.

Let $G,H$ be permutation groups on $[p]$ and $[q]$, and
$\mathcal{B}(G),\mathcal{B}(H)$ be some bases of $\mathcal{E}(G)$ and
$\mathcal{E}(H)$, respectively. Then, the following sets form bases for
$G\oplus H$, $G\otimes H$, and $H\wr G$:

$\displaystyle\mathcal{B}(G\,{\oplus}\,H)$
$\displaystyle=\\{A\,{\oplus}\,\mathbf{0}_{q},\mathbf{0}_{p}\,{\oplus}\,B\,{:}\,A\,{\in}\,\mathcal{B}(G),B\,{\in}\,\mathcal{B}(H)\\}$
$\displaystyle\qquad\cup\\{\mathbf{1}_{O\times(p+O^{\prime})},\mathbf{1}_{(p+O^{\prime})\times
O}\,{:}\,O\,{\in}\,\mathcal{O}(G),O^{\prime}\,{\in}\,\mathcal{O}(H)\\};$
$\displaystyle\mathcal{B}(G\,{\otimes}\,H)$
$\displaystyle=\\{A\,{\otimes}\,B\,{:}\,A\,{\in}\,\mathcal{B}(G),B\,{\in}\,\mathcal{B}(H)\\};\
$ $\displaystyle\mathcal{B}(H\,{\wr}\,G)$
$\displaystyle{}=\\{A\,{\otimes}\,\mathbf{1}_{O^{\prime}\times
O^{\prime\prime}}\,{:}\,A\,{\in}\,\mathcal{B}(G),A_{i,i}=0\ \text{ for
}i\in[p]\text{ and }O^{\prime},O^{\prime\prime}\,{\in}\,\mathcal{O}(H)\\}$
$\displaystyle\qquad{}\cup\\{I_{O}\,{\otimes}\,B\,{:}\,B\,{\in}\,\mathcal{B}(H),O\,{\in}\,\mathcal{O}(G)\\}.$

Here $\mathbf{0}_{m}$ is an everywhere-zero matrix in $\mathbb{R}^{m\times
m}$, and $\mathbf{1}_{R\times S}$, $I_{R}$ are the matrices defined by
$(\mathbf{1}_{R\times S})_{i,j}=\operatorname*{\mathds{1}}_{\\{i\in R,\,j\in
S\\}}$, $(I_{R})_{i,j}=\operatorname*{\mathds{1}}_{\\{i=j,\,i\in R\\}}$ whose
shapes are defined by the context in which they are used. Here,
$\mathbf{1}_{O\times(p+O^{\prime})},\mathbf{1}_{(p+O^{\prime})\times
O}\in\mathbb{R}^{(p+q)\times(p+q)}$, $\mathbf{1}_{O^{\prime}\times
O^{\prime\prime}}\in\mathbb{R}^{q\times q}$, $I_{O}\in\mathbb{R}^{p\times p}$.
Also, $\mathcal{O}(G)=\\{\\{g(i):g\in G\\}:i\in[p]\\}$ (i.e., the set of
$G$-orbits), and $p+O=\\{p+i:i\in O\\}$.

##  3 Symmetry-aware SATNet

In this section, we present SymSATNet, which abbreviates symmetry-aware
SATNet. This variant is designed to operate when symmetries of a learning task
are known a priori (via an algorithm or a domain expert). The proofs of the
theorem and the lemma in the section are in Appendix C.

SymSATNet solves the optimisation problem of SATNet, but under the following
assumptions:

######  Assumption 3.1.

The optimisation objective $\langle C,V^{T}V\rangle$ in (1) as a map on
$V=\mathbb{R}^{k\times n}$ is invariant under a permutation group $G$, whose
action is of type in Section 2.2 (i.e., each $g\in G$ acts as a permutation on
the columns of $V$).

Continuing our convention, we denote $P_{g}$ by $g$. One immediate consequence
of Assumption 3.1 is

$\displaystyle\langle C,V^{T}V\rangle$ $\displaystyle{}=\langle C,(g\cdot
V)^{T}(g\cdot V)\rangle{}=\langle
C,(Vg^{-1})^{T}(Vg^{-1})\rangle\quad\text{for all }g\in G\text{ and
}V\in\mathbb{R}^{k\times n}.$

The next theorem re-phrases this property of the optimisation objective as
equivariance of $C$:

######  Theorem 3.2.

Let $C$ be a symmetric $n\times n$ matrix. Then,

$\langle C,V^{T}V\rangle=\langle C,(Vg^{-1})^{T}(Vg^{-1})\rangle$ (3)

for all $V\in\mathbb{R}^{k\times n}$ and $g\in G$ if $C$ as a linear map on
$\mathbb{R}^{n}$ is $G$-equivariant, that is, $Cg=gC$ for all $g\in G$.
Furthermore, if $k=n$, the converse also holds.

This theorem lets us incorporate the symmetries into the objective of SATNet,
and leaves the handling of the diagonal constraints of SATNet. The next lemma
says that those constraints require no special treatment, though, since they
are already preserved by the action of any $g\in G$.

######  Lemma 3.3.

Let $V\in\mathbb{R}^{k\times n}$ and $g\in G$. Every column of $V$ has the
$L_{2}$-norm $1$ if and only if every column of $Vg^{-1}$ has the $L_{2}$-norm
$1$.

Recall $\mathcal{E}(G)=\\{M\in\mathbb{R}^{n\times n}:Mg=gM,\,\forall g\in
G\\}$ is the vector space of $G$-equivariant matrices. Let
$\mathcal{E}(G)_{s}$ be the subset of $\mathcal{E}(G)$ containing only
symmetric matrices. When $G$ is a permutation group constructed by direct sum,
direct product, and wreath product, we can generate a basis $\mathcal{B}(G)$
of $\mathcal{E}(G)$ automatically using Theorem 2.3. Then, we can convert
$\mathcal{B}(G)$ to an orthogonal basis of $\mathcal{E}(G)_{s}$ by applying
the Gram-Schmidt orthogonalisation to $\\{B+B^{T}:B\in\mathcal{B}(G)\\}$. Let
$\mathcal{B}(G)_{s}=\\{B_{1},\,\ldots,\,B_{d}\\}$ be such an orthogonal basis
of $\mathcal{E}(G)_{s}$.

SymSATNet is SATNet where the matrix $C$ in the optimisation objective has the
form:

$C=\sum_{\alpha=1}^{d}\theta_{\alpha}B_{\alpha}$ (4)

for some scalars $\theta_{1},\ldots,\theta_{d}\in\mathbb{R}$. Note that by
this condition on the form of $C$, SymSATNet has only $d$ parameters
$\theta_{1},\ldots,\theta_{d}$, instead of $n\times m$ for some $m$ in the
original formulation of SATNet. When the learning target has enough
symmetries, $d$ is usually far smaller than $n^{2}$ or even $n$, and this
reduction brings speed-up and improved generalisation.

The forward computation of SymSATNet is precisely that of SATNet, the repeated
coordinate-wise updates until convergence, and the backward computation is the
one of SATNet extended (by the chain rule) with a step backpropagating the
derivatives $\partial l/\partial C$ to each $\partial
l/\partial\theta_{\alpha}$ for $\alpha\in[d]$.333SymSATNet is implemented
based on the SATNet code [Wang et al., 2019] available under the MIT License.

We summarise SymSATNet below using the usual notation of SATNet
($\mathcal{I}$, $\mathcal{O}$, $V_{\mathcal{I}}$, $V_{\mathcal{O}}$, and $V$):

  * •

The input is $V_{\mathcal{I}}$, the matrix of the input columns of $V$.

  * •

The parameters are $(\theta_{1},\ldots,\theta_{d})\in\mathbb{R}^{d}$. They
define the matrix $C$ by (4).

  * •

The forward computation solves the following optimisation problem using
coordinate descent, and returns $V_{\mathcal{O}}$, the matrix of the output
columns of $V$:

$\displaystyle\operatorname*{argmin}_{V_{\mathcal{O}}\in\mathbb{R}^{k\times|\mathcal{O}|}}\,\langle
C,V^{T}V\rangle\quad\text{subject to}\ \|v_{o}\|=1\ \text{for}\
o\in\mathcal{O}.$

  * •

The backward computation computes $\partial l/\partial V_{\mathcal{I}}$ and
$\partial l/\partial\theta_{\alpha}$ by (2) and the chain rule:

$\displaystyle\partial l/\partial\theta_{\alpha}$ $\displaystyle=(\partial
l/\partial\operatorname*{vec}(C))\operatorname*{vec}(B_{\alpha})=-U^{T}W\operatorname*{vec}(B_{\alpha}).$
(5)

##  4 Discovery of Symmetries

One obstacle for using SymSATNet is that the user has to specify symmetries.
We now discuss how to alleviate this issue by presenting an algorithm for
discovering candidate symmetries automatically.

The goal of our algorithm, denoted by SymFind, is to find a permutation group
$G$ that captures the symmetries of an unknown learning target and is
expressible by the following grammar:
$G::=\,\mathcal{I}_{m}\,\mid\,\mathbb{Z}_{m}\,\mid\,\mathcal{S}_{m}\,\mid\,G\oplus
G\,\mid\,G\otimes G\,\mid\,G\wr G\ $ for $m\in\mathbb{N}$. The
$\mathcal{I}_{m}$ denotes the trivial group containing only the identity
permutation on $[m]$, and $\mathbb{Z}_{m}$ denotes the group of cyclic
permutations on $[m]$, each of which maps $i\in[m]$ to
$(i+n)\,\mathrm{mod}\,m$ for some $n$. The $\mathcal{S}_{m}$ is the group of
all the permutations on $[m]$. The last three cases are direct sum, direct
product, and wreath product (see Definition 2.2). They describe three ways of
decomposing a group $G$ into smaller parts. Having such a decomposition of $G$
brings the benefit to recursively and efficiently compute a basis of
$G$-equivariant linear maps.

The design of SymFind is based on our empirical observation that a softened
version of symmetries often emerges in the parameter matrix $C$ of the
original SATNet during training. Even in the early part of training, many
entries of $C$ share similar values, and there is a large-enough group $G$
with $Cg\approx gC$ for all $g\in G$, which intuitively means that $G$
captures symmetries of $C$. Furthermore, we observed, such $G$ often consists
of symmetries of the learning target. This observation suggests an algorithm
that takes $C$ as input and finds such $G$ in our grammar or its slight
extension.

Algorithm 1 SymFind with a threshold $\lambda>0$

1: Input: $M\in\mathbb{R}^{m\times m}$ Output: $(G,\sigma)$

2: if
$\|\operatorname*{prj}(\operatorname*{grp}(\mathcal{S}_{m},\mathrm{id}_{m}),M)-M\|_{F}\leq\lambda$
then

3: return $(\mathcal{S}_{m},\mathrm{id}_{m})$

4: end if

5: $\mathcal{A}\leftarrow\\{(\mathcal{I}_{m},\mathrm{id}_{m})\\}$

6: if
$\|\operatorname*{prj}(\operatorname*{grp}(\mathbb{Z}_{m},\mathrm{id}_{m}),M)-M\|_{F}\leq\lambda$
then

7:
$\mathcal{A}\leftarrow\mathcal{A}\cup\\{(\mathbb{Z}_{m},\mathrm{id}_{m})\\}$

8: end if

9: $(G^{\prime},\sigma^{\prime})\leftarrow\textsc{SumFind}(M)$;

10: $\mathcal{A}\leftarrow\mathcal{A}\cup\\{(G^{\prime},\sigma^{\prime})\\}$

11: for every divisor $p$ of $m$ do

12:
$(G^{\prime\prime},\sigma^{\prime\prime})\,{\leftarrow}\,\textsc{ProdFind}(M,p)$;

13:
$\mathcal{A}\,{\leftarrow}\,\mathcal{A}\,{\cup}\,\\{(G^{\prime\prime},\sigma^{\prime\prime})\\}$

14: end for

15:
$(G,\sigma)\leftarrow\operatorname*{argmin}_{(G,\sigma)\in\mathcal{A}}\dim(\mathcal{E}(\operatorname*{grp}(G,\sigma)))$

16: return $(G,\sigma)$

The input of SymFind is a matrix $M\in\mathbb{R}^{m\times m}$. As previously
explained, when SymFind is called at the top level, it receives as input the
parameter $C$ of SATNet learnt by a fixed number of training steps. However,
subsequent recursive calls to SymFind may have input $M$ different from $C$.
Then, SymFind returns a group $G$ in our grammar and a permutation $\sigma$ on
$[m]$, together defining a permutation group on $[m]$:

$\displaystyle\textsc{SymFind}(M)=(G,\sigma);$
$\displaystyle\operatorname*{grp}(G,\sigma)=\\{\sigma\circ
g\circ\sigma^{-1}:g\in G\\},$

where $\circ$ is the composition of permutations. When $G$ is decomposed into,
say $G_{1}\oplus G_{2}$, the $\sigma$ specifies which indices in $[m]$ get
permuted by $G_{1}$ and $G_{2}$. Once top-level SymFind returns $(G,\sigma)$,
we construct $\mathcal{B}(\operatorname*{grp}(G,\sigma))_{s}$, as in Section 3
with a minor adjustment with $\sigma$.444We construct
$\mathcal{B}(\operatorname*{grp}(G,\sigma))=\\{\sigma
B\sigma^{T}:B\in\mathcal{B}(G)\\}$, which is an orthogonal basis for
$\mathcal{E}(\operatorname*{grp}(G,\sigma))$.

Algorithm 1 describes SymFind, where $\mathrm{id}_{m}$ is the identity
permutation on $[m]$, $\|{}\cdot{}\|_{F}$ is the Frobenius norm,
$\dim(\mathcal{V})$ is the dimension of a vector space $\mathcal{V}$, and the
Reynolds operator $\operatorname*{prj}$ projects a matrix
$M\in\mathbb{R}^{m\times m}$ orthogonally to the subspace of $G$-equivariant
$m\times m$ matrices:

$\operatorname*{prj}(G,M)=\frac{1}{|G|}\sum_{g\in G}gMg^{T},$

so that $\|\operatorname*{prj}(G,M)-M\|_{F}$ computes the $L^{2}$ distance
between the matrix $M$ and the space $\mathcal{E}(G)$.

In the lines 2-4, the algorithm first checks whether $\mathcal{S}_{m}$ models
symmetries of the input $M$ accurately. If so, the algorithm returns
$(\mathcal{S}_{m},\mathrm{id}_{m})$. Otherwise, it assumes that an appropriate
group for $M$’s symmetries is one of the remaining cases in the grammar, and
constructs a list $\mathcal{A}$ of candidates intially containing the trivial
group $(\mathcal{I}_{m},\mathrm{id}_{m})$. In the lines 6-8, the algorithm
adds a pair $(\mathbb{Z}_{m},\mathrm{id}_{m})$ to $\mathcal{A}$ if it
approximates $M$’s symmetries well. In the lines 9-10, the algorithm calls the
subroutine SumFind which finds
$\operatorname*{grp}(G^{\prime},\sigma^{\prime})$ with
$G^{\prime}=\bigoplus_{i}G_{i}^{\prime}$ that approximates $M$’s symmetries
well. In the lines 11-14, the algorithm calls the other subroutine ProdFind
for every divisor $p$ of $m$. For each $p$, ProdFind finds
$\operatorname*{grp}(G^{\prime\prime},\sigma^{\prime\prime})$ with
$G^{\prime\prime}=G^{\prime\prime}_{1}\otimes G^{\prime\prime}_{2}$ (or
$G^{\prime\prime}_{2}\wr G^{\prime\prime}_{1}$) that approximates $M$’s
symmetries well, where $G^{\prime\prime}_{1}$ and $G^{\prime\prime}_{2}$ are
permutation groups on $[p]$ and $[m/p]$. Finally, in the line 15, SymFind
picks a pair $(G,\sigma)$ from the candidates $\mathcal{A}$ with the strongest
level of symmetries in the sense that the basis of
$\operatorname*{grp}(G,\sigma)$-equivariant matrices has the fewest elements.

The subroutine SumFind clusters entries of $M$ as blocks since block-shaped
clusters commonly arise in matrices equivariant with respect to a direct sum
of groups. The other subroutine ProdFind uses a technique [Van Loan and
Pitsianis, 1993] to exploit a typical pattern of Kronecker product of
matrices, and detects the presence of the pattern in $M$ by applying SVD to a
reshaped version of $M$. Each subroutine may call SymFind recursively. See
Appendix D for the details.

##  5 Experimental Results

We experimentally evaluated SymSATNet and the SymFind algorithm on the tasks
of learning rules of two problems, Sudoku and the completion problem of
Rubik’s cube. The original SATNet was used as a baseline, and both ground-
truth and automatically-discovered symmetries were used for SymSATNet. For
SymFind, we also tested its ability to recover known symmetries given randomly
generated equivariant matrices. We observed significant improvement of
SymSATNet over SATNet in various learning tasks, and also the promising
results and limitation of SymFind.

#### Sudoku problem

In Sudoku, we are asked to fill in the empty cells of a $9\times 9$ board such
that every row, every column, and each of nine $3\times 3$ blocks have all
numbers $1-9$. Let $A\in\\{0,1\\}^{9\times 9\times 9}$ be the encoding of a
full number assignment for the board where the $(i,j,k)$-th entry of $A$ is
$1$ if the $(i,j)$-th cell of the board contains $k$. In SATNet, we flatten
$A$ to the assignment on the $n=9^{3}$ boolean variables, and relax each
variable into $\mathbb{R}^{k}$, resulting $V\in\mathbb{R}^{k\times n}$ in the
objective of SATNet.

(a) Sudoku

(b) Rubik’s cube

Figure 1: Test accuracies over training epochs.

The rules of Sudoku have symmetries formalised by
$G=(\mathcal{S}_{3}\wr\mathcal{S}_{3})\otimes(\mathcal{S}_{3}\wr\mathcal{S}_{3})\otimes\mathcal{S}_{9}$.
Each of two $\mathcal{S}_{3}\wr\mathcal{S}_{3}$ refers to solution-preserving
permutations for rows and columns in Sudoku. The last $\mathcal{S}_{9}$ refers
to permutations of the assigned numbers $1-9$ in each cell. See Appendix F for
more information about the symmetry group for Sudoku.

To learn the rules of Sudoku using SymSATNet, we constructed a basis
$\mathcal{B}(G)_{s}$ as explained in Section 3. It has 18 elements, thus
SymSATNet has 18 parameters to learn.

We used 9K training and 1K test examples generated by the Sudoku generator
[Park, 2018]. Each example is a pair $(V_{\mathcal{I}},V_{\mathcal{O}})$ where
the input $V_{\mathcal{I}}$ assigns 31-42 cells (out of 81 cells) and the
output $V_{\mathcal{O}}$ specifies the remaining cells. SymSATNet was compared
with SATNet-Plain without auxiliary variables, and SATNet-300aux with 300
auxiliary variables. We used binary cross entropy loss and Adam optimizer
[Kingma and Ba, 2015], with the learning rate $\eta=2\times 10^{-3}$ for
SATNet-Plain and SATNet-300aux as the original work and $\eta=4\times 10^{-2}$
for SymSATNet. We measured test accuracy, the rate of the correctly-solved
Sudoku instances by the forward computations. We reported the average over 10
runs with $95\%$ confidence interval.

The results are in Figure 1 and Table 1. (SymSATNet-Auto refers to a variant
of SymSATNet that uses SymFind as a subroutine to find symmetries
automatically, and will be described later in this section.) Our SymSATNet
outperformed SATNet-Plain and SATNet-300aux. On average, its 100 epochs
finished $2$-$4\times$ faster in terms of wall clock than two alternatives.
due to the reduced number of iterations and avoidance of matrix operations.
See Appendix G for the efficiency of SymSATNet. Despite the speed up, the best
test accuracy of SymSATNet ($99.2\%$) was significantly better than SATNet-
Plain ($88.1\%$) and slightly better than SATNet-300aux ($97.9\%$).

Table 1: Best test accuracies during 100 epochs and average train times ($10^{2}$ sec). Additional times for automatic symmetry detection are also reported after $+$. Model | Sudoku | Cube
---|---|---
| Acc. | Time | Acc. | Time
SATNet-Plain | 88.1% | 48.0 | 55.7% | 1.8
|  $\pm$1.8% |  $\pm$0.17 |  $\pm$0.7% |  $\pm$0.01
SATNet-300aux | 97.9% | 90.3 | 56.5% | 14.0
|  $\pm$0.3% |  $\pm$0.68 |  $\pm$0.9% |  $\pm$0.12
SymSATNet | 99.2% | 25.6 | 66.9% | 1.1
|  $\pm$0.2% |  $\pm$0.14 |  $\pm$1.2% |  $\pm$0.00
SymSATNet-Auto | 99.5% | 22.7 | 68.1% | 3.4
|  $\pm$0.2% |  $+$0.14 $\pm$0.35 |  $\pm$2.8% |  $+$0.66 $\pm$0.19

#### Completion problem of Rubik’s cube

The Rubik’s cube is composed of $6$ faces, each of which has $9$ facelets. We
considered a constraint satisfaction problem where we are asked to complete
the missing facelets of the Rubik’s cube such that the resulting cube is
solvable; by moving the cube, we can make all facelets in each face have the
same colour, and no same colours appear in two faces. Let
$A\in\\{0,1\\}^{6\times 9\times 6}$ be a colour assignment of Rubik’s cube
where the $(i,j,k)$-th entry has $1$ if and only if the $j$-th facelet of the
$i$-th face has colour $k$. We formulate the optimisation objective of SATNet
for Rubik’s cube using the relaxation of $A$ to $V\in\mathbb{R}^{k\times n}$
for $n=6\times 9\times 6$.

This problem has symmetries formalised by
$G=\mathcal{R}_{54}\otimes\mathcal{R}_{6}$ on $[n]$. Here $\mathcal{R}_{54}$
and $\mathcal{R}_{6}$ are permutation groups on $[54]$ and $[6]$, each of
which captures the allowed moves of facelets, and the rotations of the whole
cube. If a colour assignment $A$ is solvable, so is the transformation of $A$
by any permutations in $G$. See Appendix F for more information about the
symmetries of this problem.

We generated a basis $\mathcal{B}(G)_{s}$ in three steps. We first created
$\mathcal{B}(\mathcal{R}_{54})$ and $\mathcal{B}(\mathcal{R}_{6})$ using the
generators of each group [Finzi et al., 2021]. Next, we combined them using
Theorem 2.3 to get $\mathcal{B}(G)$, which was converted to a symmetric
orthogonal basis $\mathcal{B}(G)_{s}$ via Gram-Schmidt. The final result has
$48$ basis elements.

We used a dataset of 9K training and 1K test examples generated by randomly
applying moves to the solution of the cube. Each example is a pair
$(V_{\mathcal{I}},V_{\mathcal{O}})$ where $V_{\mathcal{I}}$ assigns colours to
facelets except for two corner facelets, two edge facelets, and one center
facelet, and $V_{\mathcal{O}}$ specifies the colours of those five missing
facelets. In the test examples, only $V_{\mathcal{I}}$ is used. We trained
SymSATNet, SATNet-Plain, and SATNet-300aux for 100 epochs, under the same
configuration as in the Sudoku case.

The results appear in Figure 1 and Table 1. On average, the $100$-epoch
training of SymSATNet completed faster in the wall-clock time than those of
SATNet-Plain and SATNet-300aux. Also, it achieved better test accuracies
($66.9\%$) than these alternatives ($55.7\%$ and $56.5\%$). Note that unlike
Sudoku, the test accuracy of SATNet-300aux was only marginally better than
that of SATNet-Plain, which indicates that both suffered from the overfitting
issue. Note also the sharp increase in the training time of SATNet-300aux.
These two indicate that adding auxiliary variables is not so effective for the
completion problem for Rubik’s cube, while exploiting symmetries is still
useful.

#### Automatic discovery of symmetries

To test the effectiveness of SymFind, we tested whether SymFind could find
proper symmetries in Sudoku and Rubik’s cube. We applied SymFind to the
parameter $C$ of SATNet-Plain in $T$-th training epoch, where $T=10$ for
Sudoku and $T=20$ for Rubik’s cube. For Sudoku, SymFind always recovered the
full symmetries with
$G=(\mathcal{S}_{3}\wr\mathcal{S}_{3})\otimes(\mathcal{S}_{3}\wr\mathcal{S}_{3})\otimes\mathcal{S}_{9}$
in our 10 trials. For Rubik’s cube, the group of full symmetries is
$\operatorname*{grp}(G,\sigma)$ for
$G=((\mathcal{S}_{2}\wr\mathcal{S}_{3})\oplus(\mathcal{S}_{3}\wr\mathcal{S}_{8})\oplus(\mathcal{S}_{2}\wr\mathcal{S}_{12}))\otimes(\mathcal{S}_{2}\wr\mathcal{S}_{3})$.
SymFind recovered all the parts except $\mathcal{S}_{2}\wr\mathcal{S}_{12}$.
Instead of this, the algorithm found $\mathcal{S}_{12}\otimes\mathcal{S}_{2}$
or $\mathcal{S}_{3}\otimes\mathcal{S}_{8}$, or
$\mathcal{S}_{4}\otimes\mathcal{S}_{6}$ in our 10 trials. We manually observed
that the entries of $C$ in the corresponding part were difficult to be
clustered, violating the assumption of SymFind. This illustrates the
fundamental limitation of SymFind.

(a) Sudoku

(b) Rubik’s cube

Figure 2: Best test accuracies for noisy Sudoku and Rubik’s cube datasets.

To account for the limitation of SymFind, we refined the group $G$ of detected
symmetries to a subgroup in an additional validation step, before training
SymSATNet with those symmetries. In the validation step, we checked the
usefulness of each part $G_{i}$ of the expression of $G$ in our grammar.
Concretely, we rewrote $G$ only with the part $G_{i}$ in concern, where all
the other parts of $G$ were masked by the trivial groups $\mathcal{I}_{k}$.
After projecting $C$ with the masked groups using Reynolds operator, we
measured the improvement of accuracy of SATNet over validation examples.
Finally, we assembled only the parts $G_{i}$ that led to sufficient
improvement. For example, if
$G=(\mathbb{Z}_{3}\otimes\mathcal{S}_{4})\oplus\mathcal{S}_{5}$ is discovered
by SymFind, we consider the parts $G_{1}=\mathbb{Z}_{3}$,
$G_{2}=\mathcal{S}_{4}$, and $G_{3}=\mathcal{S}_{5}$. Then, we construct three
masked groups
$G_{1}^{\prime}=(\mathbb{Z}_{3}\otimes\mathcal{I}_{4})\oplus\mathcal{I}_{5}$,
$G_{2}^{\prime}=(\mathcal{I}_{3}\otimes\mathcal{S}_{4})\oplus\mathcal{I}_{5}$,
and
$G_{3}^{\prime}=(\mathcal{I}_{3}\otimes\mathcal{I}_{4})\oplus\mathcal{S}_{5}$,
and measure the accuracy of SATNet with $C$ projected by each $G_{i}^{\prime}$
over validation examples. If $G_{1}^{\prime}$ and $G_{3}^{\prime}$ show
accuracy improvements greater than a threshold, we combine $G_{1}$ and $G_{3}$
to form
$G^{\prime}=(\mathbb{Z}_{3}\otimes\mathcal{I}_{4})\oplus\mathcal{S}_{5}$,
which is then used to train SymSATNet.

We used 8K training, 1K validation, and 1K test examples to train SymSATNet
with symmetries found by SymFind and the validation step. We denote these runs
by SymSATNet-Auto. We took a group $G$ discovered by SymFind in $T$-th
training epoch (with the same $T$ above) and constructed its subgroup
$G^{\prime}$ via the validation step. SymSATNet was then trained after being
initialised by the projection of $C$ with $G^{\prime}$. The other
configurations are the same as before.

As shown in Figure 1 and Table 1, SymSATNet-Auto performed the best for Sudoku
($99.5\%$) and Rubik’s cube ($68.1\%$) better than even SymSATNet. During the
10 trials with Sudoku, SymSATNet-Auto was always given the full symmetries in
Sudoku. For Rubik’s cube, when SymSATNet-Auto was given correct subgroups
(e.g.,
$((\mathcal{S}_{2}\wr\mathcal{S}_{3})\oplus(\mathcal{S}_{3}\wr\mathcal{S}_{8})\oplus(\mathcal{I}_{4}\otimes\mathcal{I}_{6}))\otimes(\mathcal{S}_{2}\wr\mathcal{S}_{3})$,
$((\mathcal{S}_{2}\wr\mathcal{S}_{3})\oplus(\mathcal{S}_{3}\wr(\mathcal{S}_{2}\wr\mathcal{S}_{4}))\oplus(\mathcal{I}_{8}\otimes\mathcal{I}_{3}))\otimes(\mathcal{S}_{2}\wr\mathcal{S}_{3})$),
then it performed even better than SymSATNet. In two of the 10 trials,
slightly incorrect symmetries were exploited, but it outperformed SATNet-Plain
and SATNet-300aux. These results show the partial symmetries of subgroups
derived by the validation step are still useful, even when they are slightly
inaccurate.

#### Robustness to noise

We tested robustness of SymSATNet and SymSATNet-Auto to noise by noise-
corrupted datasets. We generated noisy Sudoku and Rubik’s cube datasets where
each training example is corrupted with noise; it alters the value of a random
cell or the colour of a random facelet to a random value other than the
original. We used noisy datasets with 0-3 corrupted instances to measure the
test accuracy, and tried 10 runs for each dataset to report the average and
$95\%$ confidence interval. All the other setups are the same as before.
Figure 2 shows the results. In both problems, SymSATNet was the most robust,
showing remarkably consistent accuracies. SymSATNet-Auto showed comparable
robustness to SATNet-300aux in noisy Sudoku, but outperformed the two
baselines in noisy Rubik’s cube.

Next, to show the robustness of SymFind, we applied it to restore permutation
groups $G$ from noise-corrupted $G$-equivariant symmetric matrices $M$. We
picked $(G_{i},\,\sigma_{i})$ for $i\in[4]$ where $\sigma_{1}$, $\sigma_{2}$,
$\sigma_{3}$ are random permutations on $[15]$, $[30]$, $[12]$, and
$\sigma_{4}$ is the identity permutation on $[8]$, and

$\displaystyle G_{1}$
$\displaystyle=\mathbb{Z}_{3}\oplus\mathbb{Z}_{3}\oplus\mathbb{Z}_{3}\oplus\mathbb{Z}_{3}\oplus\mathbb{Z}_{3},$
$\displaystyle G_{2}$ $\displaystyle=\mathcal{S}_{3}\wr\mathcal{S}_{10},$
$\displaystyle G_{3}$
$\displaystyle=(\mathcal{S}_{3}\wr\mathcal{S}_{3})\oplus\mathbb{Z}_{3},$
$\displaystyle G_{4}$
$\displaystyle=\mathcal{S}_{2}\otimes\mathcal{S}_{2}\otimes\mathcal{S}_{2}.$

Then, we generated $\operatorname*{grp}(G_{i},\sigma_{i})$-equivariant
symmetric matrices $M_{i}$ by projecting random matrices with standard normal
entries into the space
$\mathcal{E}(\operatorname*{grp}(G_{i},\sigma_{i}))_{s}$. Then, Gaussian
noises from $\mathcal{N}(0,\omega^{2})$ for $\omega=5\times 10^{-3}$ are added
to $M_{i}$’s entries, and the resulting matrix $M^{\prime}_{i}$ is given to
SymFind.

Table 2: Full accuracies and partial accuracies of SymFind for given groups over 1K runs. Group | Full Acc. | Partial Acc.
---|---|---
$\bigoplus_{i=1}^{5}\mathbb{Z}_{3}$ | 76.6% | 79.2%
$\mathcal{S}_{3}\wr\mathcal{S}_{10}$ | 60.3% | 79.9%
$(\mathcal{S}_{3}\wr\mathcal{S}_{3})\oplus\mathbb{Z}_{3}$ | 77.5% | 87.0%
$\mathcal{S}_{2}\otimes\mathcal{S}_{2}\otimes\mathcal{S}_{2}$ | 93.5% | 94.3%

For each $(G_{i},\sigma_{i})$, we repeatedly generated $M^{\prime}_{i}$ and
ran SymFind on $M^{\prime}_{i}$ for 1K times, and measured the portion where
SymFind recovered $(G_{i},\sigma_{i})$ exactly (full accuracy), and also the
portion of cases where SymFind returned a subgroup of $(G_{i},\sigma_{i})$
which is not the trivial group $\mathcal{I}_{m}$ (partial accuracy). As Table
2 shows, the measured full accuracies were in the range of $60.3-93.5\%$, and
the partial accuracies were in the range of $79.2-94.3\%$. These results show
the ability of SymFind to recover meaningful and sometimes full symmetries.

(a) Normal Sudoku $\rightarrow$ Hard Sudoku

(b) Normal Rubik’s cube $\rightarrow$ Hard Rubik’s cube

(c) Easy Sudoku $\rightarrow$ Hard Sudoku

(d) Easy Rubik’s cube $\rightarrow$ Hard Rubik’s cube

Figure 3: Transfer learning with various difficulties of training and test
examples. For both problems, normal and easy examples were used to train, and
hard examples were used to test each model.

#### Transfer learning

To test the transferability of SymSATNet, we generated Sudoku and Rubik’s cube
datasets with varying difficulties, where each dataset consisted of 9K
training and 1K test examples. For SymSATNet-Auto, we split the 9K training
examples into 8K training and 1K validation examples. We used three levels of
difficulties for Sudoku and Rubik’s cube (easy, normal, hard), based on the
number of missing cells for Sudoku or missing facelets for Rubik’s cube. The
input part of Sudoku examples was generated with 21 (easy), or 31 (normal), or
41 masked cells (hard), and the input part of Rubik’s cube examples was
generated with 3 (easy), or 4 (normal), or 5 missing facelets (hard). For both
problems, we used the training examples of easy or normal datasets for
training, and the test examples of hard datasets for testing. We repeated
every task in this experiment five times. Here we report the average test
accuracies and $95\%$ confidence interval.

Figure 3 shows test accuracies throughout 100 epochs in the four types of
transfer learning tasks. SymSATNet achieved the best result in the whole
tasks; it succeeded in solving hard problems after learning from easier
examples. As Figures 3(a) and 3(c) indicate, SymSATNet-Auto exploited the full
group symmetries in Sudoku. For Rubik’s cube, Figure 3(b) shows SymSATNet-Auto
achieved better performance over the baselines by finding partial symmetries.
These results show the promise of SymSATNet and SymSATNet-Auto to learn
transferable rules even from easier examples.

Note that for the easy Rubik’s cube dataset, SymSATNet-Auto showed poor
performance (Figure 3(d)). The poor performance comes from the violation of
the assumption of SymFind; the group symmetries sometimes did not emerge in
SATNet in this case. In three out of five trials with the easy Rubik’s cube
dataset, SATNet learnt nothing while producing the $0\%$ test accuracy, and
SymFind returned the trivial group which equated SymSATNet-Auto with SATNet-
Plain. In the remaining two trials, SATNet learnt correct rules, and SymFind
and the validation step found correct partial symmetries, which led to the
improved performance. These results exhibit the fundamental limitation of
SymSATNet-Auto, whose performance strongly depends on the original SATNet.

##  6 Conclusion

We presented SymSATNet, that is capable of exploiting symmetries of the rules
or constraints to be learnt by SATNet. We also described the SymFind algorithm
for automatically discovering symmetries from the parameter $C$ of the
original SATNet at a fixed training epoch, which is based on our empirical
observation that symmetries emerge during training as duplicated or similar
entries of $C$. Our experimental evaluations with two rule-learning problems
related to Sudoku and Rubik’s cube show the benefit of SymSATNet and the
promise and limitation of SymFind. Although components of SymFind are
motivated by the theoretical analysis of the space of equivariant matrices,
such as Theorem 2.2, SymFind lacks a theoretical justification on its overall
performance. One interesting future direction is to fill in this gap by
identifying when symmetries emerge during the training of SATNet and proving
probabilistic guarantees on when SymFind returns correct group symmetries.

## Acknowledgements

This work was supported by the Engineering Research Center Program through the
National Research Foundation of Korea (NRF) funded by the Korean Government
MSIT (NRF-2018R1A5A1059921) and also by the Institute for Basic Science
(IBS-R029-C1).

## References

  * Aloul et al. [2003] Fadi A Aloul, Arathi Ramani, Igor L Markov, and Karem A Sakallah.  Solving difficult instances of boolean satisfiability in the presence of symmetry.  _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , 22(9):1117–1137, 2003.
  * Aloul et al. [2006] Fadi A Aloul, Karem A Sakallah, and Igor L Markov.  Efficient symmetry breaking for boolean satisfiability.  _IEEE Transactions on Computers_ , 55(5):549–558, 2006.
  * Basu et al. [2021] Sourya Basu, Akshayaa Magesh, Harshit Yadav, and Lav R Varshney.  Group equivariant neural architecture search via group decomposition and reinforcement learning.  _arXiv e-prints_ , pages arXiv–2104, 2021.
  * Benhamou et al. [2010] Belaıd Benhamou, Tarek Nabhani, Richard Ostrowski, and Mohamed Réda Saıdi.  Dynamic symmetry breaking in the satisfiability problem.  In _Proceedings of the 16th international conference on Logic for Programming, Artificial intelligence, and Reasoning. LPAR-16, Dakar, Senegal (April 25-may 1, 2010)_ , 2010.
  * Cingillioglu and Russo [2019] Nuri Cingillioglu and Alessandra Russo.  Deeplogic: Towards end-to-end differentiable logical reasoning.  In _Proceedings of the AAAI 2019 Spring Symposium on Combining Machine Learning with Knowledge Engineering (AAAI-MAKE 2019)_ , 2019.
  * Cohen and Welling [2016] Taco S. Cohen and Max Welling.  Group equivariant convolutional networks.  In _Proceedings of ICML’16_ , pages 2990–2999, 2016.
  * Crawford [1992] James Crawford.  A theoretical analysis of reasoning by symmetry in first-order logic.  In _AAAI Workshop on Tractable Reasoning_ , pages 17–22. Citeseer, 1992.
  * Crawford et al. [1996] James M. Crawford, Matthew L. Ginsberg, Eugene M. Luks, and Amitabha Roy.  Symmetry-breaking predicates for search problems.  In _Proceedings of the Fifth International Conference on Principles of Knowledge Representation and Reasoning_ , pages 148–159. Morgan Kaufmann Publishers Inc., 1996.
  * Darga et al. [2004] Paul T. Darga, Mark H. Liffiton, Karem A. Sakallah, and Igor L. Markov.  Exploiting structure in symmetry detection for CNF.  In _Proceedings of the 41st Annual Design Automation Conference_ , pages 530–534. Association for Computing Machinery, 2004.
  * Darga et al. [2008] Paul T Darga, Karem A Sakallah, and Igor L Markov.  Faster symmetry discovery using sparsity of symmetries.  In _2008 45th ACM/IEEE Design Automation Conference_ , pages 149–154. IEEE, 2008.
  * Dehmamy et al. [2021] Nima Dehmamy, Robin Walters, Yanchen Liu, Dashun Wang, and Rose Yu.  Automatic symmetry discovery with Lie algebra convolutional network.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Devriendt et al. [2017] Jo Devriendt, Bart Bogaerts, and Maurice Bruynooghe.  Symmetric explanation learning: Effective dynamic symmetry handling for SAT.  In _International Conference on Theory and Applications of Satisfiability Testing_ , pages 83–100. Springer, 2017.
  * Dummit and Foote [1991] David S Dummit and Richard M Foote.  _Abstract algebra_ , volume 1999.  Prentice Hall Englewood Cliffs, NJ, 1991.
  * Evans and Grefenstette [2018] Richard Evans and Edward Grefenstette.  Learning explanatory rules from noisy data.  _J. Artif. Intell. Res._ , 61:1–64, 2018.
  * Finzi et al. [2021] Marc Finzi, Max Welling, and Andrew Gordon Wilson.  A practical method for constructing equivariant multilayer perceptrons for arbitrary matrix groups.  In _Proceedings of ICML’21_ , pages 3318–3328, 2021.
  * Kingma and Ba [2015] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _CoRR_ , abs/1412.6980, 2015.
  * Majumdar et al. [2012] Subrata Majumdar, Kalyan Kumar Dey, and Mohd Altab Hossain.  Direct product and wreath product of transformation semigroups.  _GANIT: Journal of Bangladesh Mathematical Society_ , 31:1–7, Apr. 2012.
  * Maron et al. [2019] Haggai Maron, Heli Ben-Hamu, Nadav Shamir, and Yaron Lipman.  Invariant and equivariant graph networks.  In _Proceedings of ICLR’19_ , 2019.
  * McKay [1976] Brendan D. McKay.  Backtrack programming and the graph isomorphism problem.  1976\.
  * McKay [1981] Brendan D. McKay.  Practical graph isomorphism.  _Congressus Numerantium_ , 30:45–87, 1981.
  * Park [2018] Kyubyong Park.  Can convolutional neural networks crack sudoku puzzles?  https://github.com/Kyubyong/sudoku, 2018.
  * Selsam et al. [2019] Daniel Selsam, Matthew Lamm, Benedikt Bünz, Percy Liang, Leonardo de Moura, and David L. Dill.  Learning a SAT solver from single-bit supervision.  In _International Conference on Learning Representations_ , 2019.
  * Topan et al. [2021] Sever Topan, David Rolnick, and Xujie Si.  Techniques for symbol grounding with SATNet.  In _Advances in Neural Information Processing Systems_ , volume 34, pages 20733–20744. Curran Associates, Inc., 2021.
  * Van Loan and Pitsianis [1993] Charles F Van Loan and Nikos Pitsianis.  Approximation with kronecker products.  In _Linear algebra for large scale and real-time applications_ , pages 293–314. Springer, 1993.
  * Wang et al. [2018] Po-Wei Wang, Wei-Cheng Chang, and J. Zico Kolter.  The mixing method: low-rank coordinate descent for semidefinite programming with diagonal constraints.  _CoRR_ , abs/1706.00476, 2018.
  * Wang et al. [2019] Po-Wei Wang, Priya L. Donti, Bryan Wilder, and J. Zico Kolter.  Satnet: Bridging deep learning and logical reasoning using a differentiable satisfiability solver.  In _Proceedings of ICML’19_ , pages 6545–6554, 2019.  URL https://github.com/locuslab/SATNet.git.
  * Wang et al. [2020] Renhao Wang, Marjan Albooyeh, and Siamak Ravanbakhsh.  Equivariant networks for hierarchical structures.  In _Proceedings of NeurIPS’20_ , pages 13806–13817, 2020.
  * Yang et al. [2017] Fan Yang, Zhilin Yang, and William W. Cohen.  Differentiable learning of logical rules for knowledge base reasoning.  In _Proceedings of NeurIPS’17_ , pages 2319–2328, 2017.

Figure 4: The overall procedure of SATNet, SymSATNet, and SymFind. The
original SATNet takes as input the assignment pairs in the dataset, and learns
the parameter matrix $C$ which describes the logical rules to solve. Our
SymFind algorithm receives the learnt parameter $C$ of SATNet, and uses the
subroutine algorithms SumFind and ProdFind and the Reynolds operator
$\operatorname*{prj}$ to find the candidate groups, and returns the strongest
group $G$ among the candidates. Finally, SymSATNet exploits the groups
symmetries in logical rules provided by domain experts, or automatically
detected group $G$ by SymFind.

##  Appendix A Notation

For natural numbers $n_{1}$ and $n_{2}$ with $n_{1}\leq n_{2}$, we write
$[n_{1}:n_{2}]$ for the set $\\{n_{1},n_{1}+1,\ldots,n_{2}\\}$. For an
$m\times m$ matrix $M$, $1\leq a\leq b\leq m$, and $1\leq c\leq d\leq m$, we
define $M[a:b,c:d]$ by a $(b-a+1)\times(d-c+1)$ matrix $M$:

$\displaystyle M[a:b,c:d]_{i,j}=M_{a+i-1,c+j-1}\quad\text{for
}i\in[b-a+1]\text{ and }j\in[d-c+1].$

Also, if $H$ is a subgroup of $G$, we denote it by $H\leq G$.

##  Appendix B Proof of Theorem 2.3

In this section, we prove Theorem 2.3. Here we prove the direct sum, direct
product, and wreath product cases by an argument similar to the previous work
[Wang et al., 2020]. For the wreath product case, we slightly generalise the
previous results [Wang et al., 2020], which considered only transitive group
actions. We will use the notation from Theorem 2.3. Let $G$ and $H$ be
permutation groups over $[p]$ and $[q]$, respectively. Also, let
$A\in\mathcal{B}(G)$, $B\in\mathcal{B}(H)$, $O\in\mathcal{O}(G)$, and
$O^{\prime},O^{\prime\prime}\in\mathcal{O}(H)$.

######  Claim B.1.

The matrices $A\oplus\mathbf{0}_{q}$, $\mathbf{0}_{p}\oplus B$,
$\mathbf{1}_{O\times(p+O^{\prime})}$, and $\mathbf{1}_{(p+O)\times
O^{\prime}}$ are $(G\oplus H)$-equivariant. Also, the matrices of these types
for all possible choices of $A$, $B$, $O$, and $O^{\prime}$ form a linearly
independent set.

###### Proof.

Pick arbitrary group elements $g\in G$ and $h\in H$. Let $P_{g}$ and $P_{h}$
be the permutation matrices corresponding to $g$ and $h$, respectively.

We prove the claimed equivariance below:

$\displaystyle(P_{g}\oplus P_{h})(A\oplus\mathbf{0}_{p})$
$\displaystyle=(P_{g}A)\oplus\mathbf{0}_{p}=(AP_{g})\oplus\mathbf{0}_{p}=(A\oplus\mathbf{0}_{p})(P_{g}\oplus
P_{h}),$ $\displaystyle(P_{g}\oplus P_{h})(\mathbf{0}_{q}\oplus B)$
$\displaystyle=\mathbf{0}_{q}\oplus(P_{h}B)=\mathbf{0}_{q}\oplus(BP_{h})=(\mathbf{0}_{q}\oplus
B)(P_{g}\oplus P_{h}),$ $\displaystyle(P_{g}\oplus
P_{h})\mathbf{1}_{O\times(p+O^{\prime})}$
$\displaystyle=\mathbf{1}_{O\times(p+O^{\prime})}=\mathbf{1}_{O\times(p+O^{\prime})}(P_{g}\oplus
P_{h}),$ $\displaystyle(P_{g}\oplus P_{h})\mathbf{1}_{(p+O^{\prime})\times O}$
$\displaystyle=\mathbf{1}_{(p+O^{\prime})\times
O}=\mathbf{1}_{(p+O^{\prime})\times O}(P_{g}\oplus P_{h}).$

The third and fourth lines use the fact that
$\mathbf{1}_{O\times(p+O^{\prime})}$ and $\mathbf{1}_{(p+O^{\prime})\times O}$
are invariant under the left or right multiplication of the permutation matrix
$P_{g}\oplus P_{h}$. This holds because the orbits of $G$ are preserved by any
permutation in $G$, and those of $H$ are preserved by all permutations in $H$,
so that

$\displaystyle P_{g}\mathbf{1}_{O\times O^{\prime}}=\mathbf{1}_{O\times
O^{\prime}}=\mathbf{1}_{O\times O^{\prime}}P_{h},\qquad
P_{h}\mathbf{1}_{O^{\prime}\times O}=\mathbf{1}_{O^{\prime}\times
O}=\mathbf{1}_{O^{\prime}\times O}P_{g}.$ (6)

Next, we show the claimed linear independence by analysing the indices of the
nonzero entries of the four types of matrices in the claim.

  1. 1.

The matrices of the type $A\oplus\mathbf{0}_{q}$ for some $A\in\mathcal{B}(G)$
are linearly independent, since their $A$ parts are linearly independent.

  2. 2.

The matrices of the type $\mathbf{0}_{p}\oplus B$ for some
$B\in\mathcal{B}(H)$ are also linearly independent by similar reason.

  3. 3.

Different matrices of the type $\mathbf{1}_{O\times(p+O^{\prime})}$ for some
$O$ and $O^{\prime}$ do not share an index of a nonzero entry, since different
orbits of a group are disjoint. Thus, the matrices of this type are linearly
independent.

  4. 4.

Different matrices of the type $\mathbf{1}_{(p+O)\times O^{\prime}}$ for some
$O$ and $O^{\prime}$ do not share an index of a nonzeron entry. Thus, the
matrices of this form are linearly independent.

Also, any matrices of the above four types form a linearly independent set
because those linear combinations do not share any indices of nonzero entries.
From this and the above reasoning for each of the four types of matrices it
follows that the matrices of those four types are linearly independent, as
claimed. ∎

######  Claim B.2.

The matrix $A\otimes B$ is $(G\otimes H)$-equivariant. Also, the matrices of
this shape for all possible choices of $A$ and $B$ form a linearly independent
set.

###### Proof.

Pick arbitrary group elements $g\in G$ and $h\in H$. Let $P_{g}$ and $P_{h}$
be the permutation matrices corresponding to $g$ and $h$, respectively. Then,

$\displaystyle(P_{g}\otimes P_{h})(A\otimes B)$
$\displaystyle=(P_{g}A)\otimes(P_{h}B)=(AP_{g})\otimes(BP_{h})=(A\otimes
B)(P_{g}\otimes P_{h}).$

For linear independence, we can prove it using the fact $\langle A\otimes
B,\,A^{\prime}\otimes B^{\prime}\rangle=\langle
A,\,A^{\prime}\rangle\cdot\langle B,\,B^{\prime}\rangle$. ∎

######  Claim B.3.

Recall $A\in\mathcal{B}(G)$ and $B\in\mathcal{B}(H)$. The matrices
$A\otimes\mathbf{1}_{O^{\prime}\times O^{\prime\prime}}$ such that $A_{i,i}=0$
for $i\in[p]$ and $I_{O}\otimes B$ are $(H\wr G)$-equivariant. Also, the
matrices of these types for all possible choices of $A$, $B$, $O$,
$O^{\prime}$, and $O^{\prime\prime}$ form a linearly independent set.

###### Proof.

Pick arbitrary group elements $g\in G$ and
$\vec{h}=(h_{1},\,\ldots,\,h_{p})\in H^{p}$. Let $P_{g^{-1}}$ and $P_{h_{i}}$
be the permutation matrices corresponding to $g^{-1}$ and $h_{i}$ for
$i\in[p]$. We can express $\operatorname*{wr}(\vec{h},g)$ and $P_{g^{-1}}$ by

$\displaystyle\operatorname*{wr}(\vec{h},g)=\sum_{i=1}^{p}\mathbf{1}_{\\{i\\}\times\\{g(i)\\}}\otimes
P_{h_{i}},\qquad
P_{g^{-1}}=\sum_{i=1}^{p}\mathbf{1}_{\\{i\\}\times\\{g(i)\\}}.$

Then, we prove the following equivariance:

$\displaystyle\operatorname*{wr}(\vec{h},g)\left(A\otimes\mathbf{1}_{O^{\prime}\times
O^{\prime\prime}}\right)$
$\displaystyle=\left(\sum_{i=1}^{p}\mathbf{1}_{\\{i\\}\times\\{g(i)\\}}\otimes
P_{h_{i}}\right)\left(A\otimes\mathbf{1}_{O^{\prime}\times
O^{\prime\prime}}\right)$
$\displaystyle=\sum_{i=1}^{p}\left(\mathbf{1}_{\\{i\\}\times\\{g(i)\\}}A\otimes
P_{h_{i}}\mathbf{1}_{O^{\prime}\times O^{\prime\prime}}\right)$
$\displaystyle=\left(\sum_{i=1}^{p}\mathbf{1}_{\\{i\\}\times\\{g(i)\\}}\right)A\otimes\mathbf{1}_{O^{\prime}\times
O^{\prime\prime}}$ (7)
$\displaystyle=P_{g^{-1}}A\otimes\mathbf{1}_{O^{\prime}\times
O^{\prime\prime}}=AP_{g^{-1}}\otimes\mathbf{1}_{O^{\prime}\times
O^{\prime\prime}}$
$\displaystyle=A\left(\sum_{i=1}^{p}\mathbf{1}_{\\{i\\}\times\\{g(i)\\}}\right)\otimes\mathbf{1}_{O^{\prime}\times
O^{\prime\prime}}$
$\displaystyle=\sum_{i=1}^{p}\left(A\,\mathbf{1}_{\\{i\\}\times\\{g(i)\\}}\otimes\mathbf{1}_{O^{\prime}\times
O^{\prime\prime}}P_{h_{i}}\right)$ (8)
$\displaystyle=\left(A\otimes\mathbf{1}_{O^{\prime}\times
O^{\prime\prime}}\right)\left(\sum_{i=1}^{p}\mathbf{1}_{\\{i\\}\times\\{g(i)\\}}\otimes
P_{h_{i}}\right)$ $\displaystyle=\left(A\otimes\mathbf{1}_{O^{\prime}\times
O^{\prime\prime}}\right)\operatorname*{wr}(\vec{h},g),$
$\displaystyle\operatorname*{wr}(\vec{h},g)\left(I_{O}\otimes B\right)$
$\displaystyle=\left(\sum_{i=1}^{p}\mathbf{1}_{\\{i\\}\times\\{g(i)\\}}\otimes
P_{h_{i}}\right)\left(I_{O}\otimes B\right)$
$\displaystyle=\sum_{i=1}^{p}\left(\mathbf{1}_{\\{i\\}\times\\{g(i)\\}}I_{O}\otimes
P_{h_{i}}B\right)$
$\displaystyle=\sum_{i=1}^{p}\left(\mathbf{1}_{\\{i\\}\times(\\{g(i)\\}\cap
O)}\otimes P_{h_{i}}B\right)$
$\displaystyle=\sum_{i=1}^{p}\left(\mathbf{1}_{(\\{i\\}\cap
O)\times\\{g(i)\\}}\otimes BP_{h_{i}}\right)$ (9)
$\displaystyle=\sum_{i=1}^{p}\left(I_{O}\mathbf{1}_{\\{i\\}\times\\{g(i)\\}}\otimes
BP_{h_{i}}\right)$ $\displaystyle=\left(I_{O}\otimes
B\right)\left(\sum_{i=1}^{p}\mathbf{1}_{\\{i\\}\times\\{g(i)\\}}\otimes
P_{h_{i}}\right)$ $\displaystyle=\left(I_{O}\otimes
B\right)\operatorname*{wr}(\vec{h},g).$

Here, (7) and (8) use the same argument in (6), and (9) uses the fact that
$i\in O\iff g(i)\in O$ for any $g\in G$ and $O\in\mathcal{O}(G)$. For linear
independence, we again use the fact $\langle A\otimes B,\,A^{\prime}\otimes
B^{\prime}\rangle=\langle A,\,A^{\prime}\rangle\cdot\langle
B,\,B^{\prime}\rangle$ to show the orthogonality of all possible matrices of
types $A\otimes\mathbf{1}_{O^{\prime}\times O^{\prime\prime}}$ and
$I_{O}\otimes B$. ∎

######  Claim B.4.

The number of basis elements of $G\oplus H$, $G\otimes H$, and $H\wr G$ can be
computed as follows:

$\displaystyle|\mathcal{B}(G\oplus H)|$
$\displaystyle=|\mathcal{B}(G)|+|\mathcal{B}(H)|+2|\mathcal{O}(G)||\mathcal{O}(H)|,$
$\displaystyle|\mathcal{B}(G\otimes H)|$
$\displaystyle=|\mathcal{B}(G)||\mathcal{B}(H)|,$
$\displaystyle|\mathcal{B}(H\wr G)|$
$\displaystyle=|\mathcal{O}(G)||\mathcal{B}(H)|+\left(|\mathcal{B}(G)|-|\mathcal{O}(G)|\right)|\mathcal{O}(H)|^{2}.$

###### Proof.

First, we derive some general fact on a finite permutation group. Let
$\mathbf{G}$ be a permutation group on $[r]$ for some $r$. Consider the action
of $\mathbf{G}$ on the space of $r$-dimensional vectors
$\mathcal{X}=\mathbb{R}^{r}$ by the row permutation action $g\cdot v=P_{g}v$
for any $g\in\mathbf{G}$. Define
$\mathcal{F}(\mathbf{G})=\\{v\in\mathcal{X}:g\cdot v=v,\,\forall
g\in\mathbf{G}\\}$. Consider the following linear operator

$\displaystyle\phi_{\mathbf{G}}$ $\displaystyle:\mathcal{X}\to\mathcal{X},$
$\displaystyle\phi_{\mathbf{G}}(v)$
$\displaystyle=\dfrac{1}{|\mathbf{G}|}\sum\limits_{g\in\mathbf{G}}g\cdot v,$

which can also be represented as the following matrix:

$\phi_{\mathbf{G}}=\dfrac{1}{|\mathbf{G}|}\sum\limits_{g\in\mathbf{G}}P_{g}.$

The operator $\phi_{\mathbf{G}}$ is a projection map, since

$\displaystyle\phi_{\mathbf{G}}(\phi_{\mathbf{G}}(v))$
$\displaystyle=\dfrac{1}{|\mathbf{G}|^{2}}\sum\limits_{g_{1}\in\mathbf{G}}\sum\limits_{g_{2}\in\mathbf{G}}g_{1}\cdot(g_{2}\cdot
v)$
$\displaystyle=\dfrac{1}{|\mathbf{G}|^{2}}\sum\limits_{g_{1}\in\mathbf{G}}\sum\limits_{g\in\mathbf{G}}g\cdot
v$ $\displaystyle=\dfrac{1}{|\mathbf{G}|}\sum\limits_{g\in G}g\cdot v.$

Also, we have $\operatorname*{im}(\phi_{\mathbf{G}})=\mathcal{F}(\mathbf{G})$
where $\operatorname*{im}(f)$ is the image of the function $f$. Now, by noting
that a linear projection map has only eigenvalues $0$ and $1$, the dimension
of $\mathcal{F}(\mathbf{G})$ can be computed by the sum of eigenvalues of
$\phi_{\mathbf{G}}$, i.e., $\operatorname*{tr}(\phi_{\mathbf{G}})$. Also,
using Burnside’s lemma, we can count the orbits of $\mathbf{G}$ (acting on
$[r]$) by

$\displaystyle|\mathcal{O}(\mathbf{G})|$
$\displaystyle=\dfrac{1}{|\mathbf{G}|}\sum\limits_{g\in\mathbf{G}}\operatorname*{tr}(P_{g})$
$\displaystyle=\operatorname*{tr}\left(\dfrac{1}{|\mathbf{G}|}\sum\limits_{g\in\mathbf{G}}P_{g}\right)$
$\displaystyle=\operatorname*{tr}(\phi_{\mathbf{G}}).$

Putting together, we get

$\dim\mathcal{F}(\mathbf{G})=\operatorname*{tr}(\phi_{\mathbf{G}})=|\mathcal{O}(\mathbf{G})|.$
(10)

Next, we instantiate what we have just shown above for the following case that
$\mathbf{G}$ is the following group:

$\mathbf{G}_{0}^{\otimes 2}=\\{g\otimes g:g\in\mathbf{G}_{0}\\}$

for some permutation group $\mathbf{G}_{0}$ on $[n]$. Note that $\mathbf{G}$
is a permutation group on $[n^{2}]$. As explained above,
$\mathbf{G}_{0}^{\otimes 2}$ can act on the space of $n^{2}$-dimensional
vectors $\mathbb{R}^{n^{2}}$. By vectorizing matrices, we can express the
space $\mathcal{E}(\mathbf{G}_{0})$ of $\mathbf{G}_{0}$-equivariant linear
maps on $\mathbb{R}^{n}$ by

$\displaystyle\operatorname*{vec}(\mathcal{E}(\mathbf{G}_{0}))$
$\displaystyle=\left\\{\operatorname*{vec}(M):P_{g}MP_{g}^{T}=M,\,\forall
g\in\mathbf{G}_{0}\right\\}$
$\displaystyle=\left\\{\operatorname*{vec}(M):(P_{g}\otimes
P_{g})\operatorname*{vec}(M)=\operatorname*{vec}(M),\,\forall
g\in\mathbf{G}_{0}\right\\}$
$\displaystyle=\left\\{\operatorname*{vec}(M):(g\otimes
g)\cdot\operatorname*{vec}(M)=\operatorname*{vec}(M),\,\forall
g\in\mathbf{G}_{0}\right\\}$
$\displaystyle=\mathcal{F}(\mathbf{G}_{0}^{\otimes 2}).$

Thus,
$|\mathcal{B}(\mathbf{G}_{0})|=\dim\operatorname*{vec}(\mathcal{E}(\mathbf{G}_{0}))=\dim\mathcal{F}(\mathbf{G}_{0}^{\otimes
2})$. We now calculate the dimension of $\mathcal{F}(\mathbf{G}_{0}^{\otimes
2})$ using the relationship in (10):

$\displaystyle|\mathcal{B}(\mathbf{G}_{0})|$
$\displaystyle=\dim\mathcal{F}(\mathbf{G}_{0}^{\otimes
2})=\operatorname*{tr}(\phi_{\mathbf{G}_{0}^{\otimes 2}})$
$\displaystyle=\dfrac{1}{|\mathbf{G}_{0}^{\otimes 2}|}\sum\limits_{{g\otimes
g}\in\mathbf{G}_{0}^{\otimes 2}}\operatorname*{tr}(P_{g\otimes
g})=\frac{1}{|\mathbf{G}_{0}|}\sum\limits_{g\in\mathbf{G}_{0}}\operatorname*{tr}(P_{g})^{2}.$
(11)

Finally, we complete the proof by calculating (11) for $\mathbf{G}_{0}=G\oplus
H$, $\mathbf{G}_{0}=G\otimes H$, and $\mathbf{G}_{0}=H\wr G$:

$\displaystyle|\mathcal{B}(G\oplus H)|$ $\displaystyle=\dfrac{1}{|G\oplus
H|}\sum\limits_{g\oplus h\in G\oplus H}\operatorname*{tr}(P_{g\oplus h})^{2}$
$\displaystyle=\dfrac{1}{|G||H|}\sum\limits_{g\in G}\sum\limits_{h\in
H}(\operatorname*{tr}(P_{g})+\operatorname*{tr}(P_{h}))^{2}$
$\displaystyle=\dfrac{1}{|G|}\sum\limits_{g\in
G}\operatorname*{tr}(P_{g})^{2}+\dfrac{1}{|H|}\sum\limits_{h\in
H}\operatorname*{tr}(P_{h})^{2}+\dfrac{2}{|G||H|}\left(\sum\limits_{g\in
G}\sum\limits_{h\in
H}\operatorname*{tr}(P_{g})\operatorname*{tr}(P_{h})\right)$
$\displaystyle=|\mathcal{B}(G)|+|\mathcal{B}(H)|+2|\mathcal{O}(G)||\mathcal{O}(H)|,$
$\displaystyle|\mathcal{B}(G\otimes H)|$ $\displaystyle=\dfrac{1}{|G\otimes
H|}\sum\limits_{g\otimes h\in G\otimes H}\operatorname*{tr}(P_{g\otimes
h})^{2}$ $\displaystyle=\dfrac{1}{|G||H|}\sum\limits_{g\in G}\sum\limits_{h\in
H}\operatorname*{tr}(P_{g})^{2}\operatorname*{tr}(P_{h})^{2}$
$\displaystyle=|\mathcal{B}(G)||\mathcal{B}(H)|,$
$\displaystyle|\mathcal{B}(H\wr G)|$ $\displaystyle=\dfrac{1}{|H\wr
G|}\sum\limits_{\operatorname*{wr}(\vec{h},g)\in H\wr
G}\operatorname*{tr}(\operatorname*{wr}(\vec{h},g))^{2}$
$\displaystyle=\dfrac{1}{|G||H|^{p}}\sum\limits_{g\in
G,(h_{1},\ldots,h_{p})\in
H^{p}}\operatorname*{tr}\left(\sum\limits_{i=1}^{p}\mathbf{1}_{\\{i\\}\times\\{g(i)\\}}\otimes
P_{h_{i}}\right)^{2}$ $\displaystyle=\dfrac{1}{|G||H|^{p}}\sum\limits_{g\in
G,(h_{1},\ldots,h_{p})\in
H^{p}}\left(\sum\limits_{i=1}^{p}\operatorname*{tr}\left(\mathbf{1}_{\\{i\\}\times\\{g(i)\\}}\otimes
P_{h_{i}}\right)\right)^{2}$
$\displaystyle=\dfrac{1}{|G||H|^{p}}\sum\limits_{g\in
G,(h_{1},\ldots,h_{p})\in
H^{p}}\left(\sum\limits_{i=1}^{p}\operatorname*{\mathds{1}}\nolimits_{\\{i=g(i)\\}}\operatorname*{tr}(P_{h_{i}})\right)^{2}$
$\displaystyle=\dfrac{1}{|G||H|^{p}}\sum\limits_{g\in
G,(h_{1},\ldots,h_{p})\in
H^{p}}\left(\sum\limits_{i=1}^{p}\operatorname*{\mathds{1}}\nolimits_{\\{i=g(i)\\}}\operatorname*{tr}(P_{h_{i}})^{2}+\sum\limits_{i\neq
j}\operatorname*{\mathds{1}}\nolimits_{\\{i=g(i),j=g(j)\\}}\operatorname*{tr}(P_{h_{i}})\operatorname*{tr}(P_{h_{j}})\right)$
$\displaystyle=\dfrac{1}{|G||H|^{p}}\sum\limits_{g\in
G}\left\\{\sum\limits_{i=1}^{p}\operatorname*{\mathds{1}}\nolimits_{\\{i=g(i)\\}}\left(\sum\limits_{(h_{1},\ldots,h_{p})\in
H^{p}}\operatorname*{tr}(P_{h_{i}})^{2}\right)\right.$ $\displaystyle\hskip
62.59605pt+\left.\sum\limits_{i\neq
j}\operatorname*{\mathds{1}}\nolimits_{\\{i=g(i),j=g(j)\\}}\left(\sum\limits_{(h_{1},\ldots,h_{p})\in
H^{p}}\operatorname*{tr}(P_{h_{i}})\operatorname*{tr}(P_{h_{j}})\right)\right\\}$
$\displaystyle=\dfrac{1}{|G||H|^{p}}\sum\limits_{g\in
G}\left\\{\left(\sum\limits_{i=1}^{p}\operatorname*{\mathds{1}}\nolimits_{\\{i=g(i)\\}}\right)|H|^{p}|\mathcal{B}(H)|+\left(\sum\limits_{i\neq
j}\operatorname*{\mathds{1}}\nolimits_{\\{i=g(i),j=g(j)\\}}\right)|H|^{p}|\mathcal{O}(H)|^{2}\right\\}$
$\displaystyle=\dfrac{1}{|G|}\sum\limits_{g\in
G}\left\\{\left(\sum\limits_{i=1}^{p}\operatorname*{\mathds{1}}\nolimits_{\\{i=g(i)\\}}\right)|\mathcal{B}(H)|+\left(\sum\limits_{i=1}^{p}\operatorname*{\mathds{1}}\nolimits_{\\{i=g(i)\\}}\sum\limits_{j=1}^{p}\operatorname*{\mathds{1}}\nolimits_{\\{j=g(j)\\}}-\sum\limits_{i=1}^{p}\operatorname*{\mathds{1}}\nolimits_{\\{i=g(i)\\}}\right)|\mathcal{O}(H)|^{2}\right\\}$
$\displaystyle=\dfrac{1}{|G|}\sum\limits_{g\in
G}\left(\operatorname*{tr}(P_{g})|\mathcal{B}(H)|+\left(\operatorname*{tr}(P_{g})^{2}-\operatorname*{tr}(P_{g})\right)|\mathcal{O}(H)|^{2}\right)$
$\displaystyle=|\mathcal{O}(G)||\mathcal{B}(H)|+\left(|\mathcal{B}(G)|-|\mathcal{O}(G)|\right)|\mathcal{O}(H)|^{2}.$

∎

#### Proof of Theorem 2.3

Claims B.1, B.2, and B.3 show that each set of matrices on the right-hand side
of the equations in Theorem 2.3 consists of linearly independent matrices, and
it is contained in the corresponding space of equivariant linear maps.
Furthermore, Claim B.4 shows that the matrices in the set span the whole of
the space of equivariant linear maps, since their number coincides with the
dimension of the space. Hence, these three claims complete the proof of the
theorem.

##  Appendix C Proofs of Theorem 3.2 and Lemma 3.3

#### Proof of Theorem 3.2

Assume that $C$ is $G$-equivariant. Pick any $g\in G$. Then, $C$ preserves the
standard action of $g$ on $\mathbb{R}^{n}$ via its permutation matrix. Thus,
$Cg=gC$, where $g$ denotes the permutation matrix $P_{g}$. This equation is
equivalent to

$g^{T}Cg=C,$ (12)

since $g^{T}=g^{-1}$. Using this equality, we show that the optimisation
objective is invariant with respect to $g$’s action:

$\displaystyle\langle C,(Vg^{-1})^{T}(Vg^{-1})\rangle$
$\displaystyle{}=\langle C,(g^{-1})^{T}(V^{T}V)g^{-1}\rangle$
$\displaystyle{}=\operatorname*{tr}(C^{T}g(V^{T}V)g^{T}))$
$\displaystyle{}=\operatorname*{tr}((g^{T}C^{T}g)(V^{T}V))$
$\displaystyle{}=\operatorname*{tr}((g^{T}Cg)^{T}(V^{T}V))$
$\displaystyle{}=\langle g^{T}Cg,V^{T}V\rangle$ $\displaystyle{}=\langle
C,V^{T}V\rangle.$

Here $\operatorname*{tr}$ is the usual trace operator on matrices, the third
equality uses the cyclic property of $\operatorname*{tr}$, and the last
equality uses (12).

We move on to the proof of the converse. Assume $k=n$ and the equation (3)
holds for every $V\in\mathbb{R}^{k\times n}$ and $g\in G$. Pick any $g\in G$.
Then, for all $V\in\mathbb{R}^{k\times n}$,

$\displaystyle\langle C,V^{T}V\rangle$ $\displaystyle{}=\langle
C,(Vg^{-1})^{T}(Vg^{-1})\rangle{}=\langle g^{T}Cg,V^{T}V\rangle.$

The first equality follows from (3), and the second from our derivation above.
But the space of matrices of the form $V^{T}V$ is precisely that of $n\times
n$ symmetric positive semi-definite matrices, which is known to contain
$n(n+1)/2$ independent elements. Since $n(n+1)/2$ is precisely the dimension
of the space of $n\times n$ symmetric matrices, the above equality for every
$V\in R^{k\times n}$ implies that $C=g^{T}Cg$, that is, $Cg=gC$. This gives
the $G$-equivariance of $C$, as desired.

#### Proof of Lemma 3.3

Pick arbitrary $V\in R^{k\times n}$ and $g\in G$. Let $v_{1},\ldots,v_{n}$ be
the columns of $V$ in order, and $v^{\prime}_{1},\ldots,v^{\prime}_{n}$ be
those of $Vg^{-1}$ in order. Then, $\|v_{i}\|=1$ for every $i\in[n]$ if and
only if $\|v_{g^{-1}(i)}\|=1$ for all $i\in[n]$. The latter is equivalent to
the statement that $\|v^{\prime}_{i}\|=1$ for every $i\in[n]$.

##  Appendix D Subroutines of SymFind

In this section, we describe the subroutines SumFind and ProdFind of our
SymFind algorithm.

###  D.1 Detection of direct sum

We recall the distributivity law of group constructors [Dummit and Foote,
1991, Majumdar et al., 2012]:

$\displaystyle G\otimes(H\oplus H^{\prime})$ $\displaystyle=(G\otimes
H)\oplus(G\otimes H^{\prime}),$ $\displaystyle(G\oplus G^{\prime})\otimes H$
$\displaystyle=(G\otimes H)\oplus(G^{\prime}\otimes H),$ $\displaystyle
G\wr(H\oplus H^{\prime})$ $\displaystyle=(G\wr H)\oplus(G\wr H^{\prime}),$
$\displaystyle(G\oplus G^{\prime})\wr H$ $\displaystyle=(G\wr
H)\oplus(G^{\prime}\wr H).$

By this law and the fact that
$\mathcal{I}_{m}=\bigoplus\limits_{i=1}^{m}\mathbb{Z}_{1}$ for all $m$, any
permutation group $G$ definable in our grammar can be expressed as a direct
sum

$\displaystyle G=\bigoplus_{i=1}^{N}{H_{(i,1)}\diamond\cdots\diamond
H_{(i,n_{i})}}$ (13)

where each $H_{(i,j)}\in\\{\mathbb{Z}_{m(i,j)},\,\mathcal{S}_{m(i,j)}\\}$ and
$\diamond$ is either $\otimes$ or $\wr$. Let

$H_{i}={H_{(i,1)}\diamond\cdots\diamond H_{(i,n_{i})}},$

and $p_{i}$ be the natural number such that $H_{i}$ is a permutation group on
$[p_{i}]$. Since both $\mathbb{Z}_{m}$ and $\mathcal{S}_{m}$ have only one
orbit for all $m$, each $H_{i}$ also has only one orbit. Hence, by Theorem
2.3,

$\mathcal{B}(G)=\bigcup_{i,j\in[N]}B_{ij},$

where

$\displaystyle B_{ii}$
$\displaystyle=\\{\mathbf{0}_{p_{1}+\cdots+p_{i-1}}\oplus
A\oplus\mathbf{0}_{p_{i+1}+\cdots+p_{N}}:A\in\mathcal{B}(H_{i})\\}\text{ for
}i\in[N],\text{ and}$ $\displaystyle B_{ij}$
$\displaystyle=\\{\mathbf{1}_{(p_{1}+\cdots+p_{i-1}+[p_{i}])\times(p_{1}+\cdots+p_{j-1}+[p_{j}])}\\}\text{
for }i,j\in[N]\text{ with }i\neq j.$

This means that all off-diagonal blocks of any $G$-equivariant matrices are
constant matrices, just like the matrix shown in (b) of Figure 5.

Given a matrix $M\in\mathbb{R}^{m\times m}$, the subroutine SumFind aims at
finding a permutation $\sigma:[m]\to[m]$ and the split $m=p_{1}+\ldots+p_{N}$
for $p_{1},\ldots,p_{N}>0$ such that $P_{\sigma}^{T}MP_{\sigma}$ is
approximately $G$-equivariant for some permutation group $G$ on $[m]$ and this
group $G$ has the form in (13) where $H_{i}={H_{(i,1)}\diamond\cdots\diamond
H_{(i,n_{i})}}$ is a permutation group on $[p_{i}]$. The result of the
subroutine is $(G,\sigma)$.

(a) Original matrix $M$

(b) Permuted matrix $M_{\sigma^{-1}}=P_{\sigma}^{T}MP_{\sigma}$

Figure 5: Block pattern which commonly appears in the equivariant matrices
under a direct-sum group. In this example,
$M_{\sigma^{-1}}\in\mathcal{E}(\mathbb{Z}_{4}\oplus\mathbb{Z}_{4}\oplus\mathbb{Z}_{4})$.
The red and blue "L-shaped" clusters have the same values.

SumFind achieves its aim using two key observations. The first is an important
implication of the $G$-equivariance condition on
$M_{\sigma^{-1}}=P_{\sigma}^{T}MP_{\sigma}$ that we mentioned above: when the
split $m=p_{1}+\ldots+p_{N}$ is used to group entries of $M_{\sigma^{-1}}$
into blocks, the $G$-equivalence of $M_{\sigma^{-1}}$ implies that the off-
diagonal blocks of $M_{\sigma^{-1}}$ are constant matrices. This implication
suggests one approach: for every permutation $\sigma$ on $[m]$, construct the
matrix $M_{\sigma^{-1}}$, and find a split $m=p_{1}+\ldots+p_{N}$ for some $N$
such that off-diagonal blocks of $M_{\sigma^{-1}}$ from this split are
constant matrices. Note that this approach is not a feasible option in
practice, though, since there are exponentially many permutations $\sigma$.
The second observation suggests a way to overcome this exponential blow-up
issue of the approach. It is that when $M_{\sigma^{-1}}$ is $G$-equivariant,
in many cases its diagonal blocks are cyclic in the sense that the $(i,j)$-th
entry of a block is the same as the $(i+1,j+1)$-th entry of the block. This
pattern can be used to search for a good permutation $\sigma$ efficiently.

SumFind works as follows. It first finds a permutation $\sigma:[m]\to[m]$ and
a split $m=p_{1}+\ldots+p_{N}$ using the process that we will explain shortly.
Figure 5 illustrates the input matrix $M$, and its permuted $M_{\sigma^{-1}}$
that has nine blocks induced by the split of $m=p_{1}+p_{2}+p_{3}$. Then,
SumFind calls SymFind recursively on each diagonal block of $M_{\sigma^{-1}}$
and gets, for every $i\in[N]$,

$(H_{i},\sigma_{i})=\textsc{SymFind}(M_{\sigma^{-1}}[p_{i-1}+1:p_{i},p_{i-1}+1:p_{i}])$

where $p_{0}=0$. Finally, SumFind returns

$\displaystyle\left(\bigoplus_{i=1}^{N}H_{i},\,\sigma\circ\left(\bigoplus_{i=1}^{N}\sigma_{i}\right)\right)$

where

$\displaystyle\sigma_{i}\oplus\sigma_{j}$
$\displaystyle:[p_{i}+p_{j}]\to[p_{i}+p_{j}],$
$\displaystyle(\sigma_{i}\oplus\sigma_{j})(a)$
$\displaystyle=\begin{cases}\sigma_{i}(a)&\text{if}\ a\in[p_{i}],\\\
\sigma_{j}(a-p_{i})+p_{i}&\text{otherwise}.\end{cases}$

We now explain the first part of SumFind that finds a permutation $\sigma$ and
a split $m=p_{1}+\ldots+p_{N}$. For simplicity, we ignore the issue of noise,
and present a simpler version that uses equality instead of approximate
equality (i.e., being close enough). We start by initialising
$\sigma=\sigma_{I}$, the identity permutation, and $M_{\sigma^{-1}}=M$. Then,
we pick an index $k_{1}$ from the set
$K=\\{k_{1}\in[2:m]:(M_{\sigma^{-1}})_{k_{1},k_{1}}=(M_{\sigma^{-1}})_{1,1}\\}$.
Then, we locate $(M_{\sigma^{-1}})_{1,k_{1}}$ in the $(1,\,2)$-th entry by
swapping the indices $2$ and $k_{1}$, which can be done by updating

$\displaystyle M_{\sigma^{-1}}$ $\displaystyle\leftarrow
P_{\sigma_{(2,k_{1})}}M_{\sigma^{-1}}P_{\sigma_{(2,k_{1})}}^{T},$
$\displaystyle\sigma$
$\displaystyle\leftarrow\sigma\circ\sigma_{(2,k_{1})}^{-1}.$

Here $\sigma_{(a,b)}(a)=b$, $\sigma_{(a,b)}(b)=a$, and $\sigma_{(a,b)}(c)=c$
for $c\neq a,\,b$. Next, we find a new index $k_{2}\in[3:m]$ which is suitable
to swap with the index $3$ in the updated $M_{\sigma^{-1}}$. As Figure 5
shows, we need to find $k_{2}$ such that
$(M_{\sigma^{-1}})_{2,k_{2}}=(M_{\sigma^{-1}})_{1,2}$,
$(M_{\sigma^{-1}})_{k_{2},k_{2}}=(M_{\sigma^{-1}})_{2,2}$, and
$(M_{\sigma^{-1}})_{k_{2},2}=(M_{\sigma^{-1}})_{2,1}$. Once we find such
$k_{2}$, we swap the indices $3$ and $k_{2}$ by updating

$\displaystyle M_{\sigma^{-1}}$ $\displaystyle\leftarrow
P_{\sigma_{(3,k_{2})}}M_{\sigma^{-1}}P_{\sigma_{(3,k_{2})}}^{T},$
$\displaystyle\sigma$
$\displaystyle\leftarrow\sigma\circ\sigma_{(3,k_{2})}^{-1}.$

We repeat this process to find an index $k_{l}\in[l+1:m]$ and the "L-shaped"
entries which preserve the cyclic pattern. If we cannot find such $k_{l}$, we
stop the process, and check (i) whether all rows in
$(M_{\sigma^{-1}})[1:l,l+1:m]$ are the same and also (ii) whether all columns
in $(M_{\sigma^{-1}})[l+1:n,1:l]$ are the same. If the answers for both (i)
and (ii) are yes, we move on to find the next diagonal block in
$M_{\sigma^{-1}}$ in the same manner.

If (i) or (ii) has a negative answer, we go back to the step right before
choosing $k_{1}$ from $K$ and resetting $\sigma$ and $M$ to the values at that
step. Then, we repeat the above process with a new choice of $k_{1}$ from $K$.
If no choice of $k_{1}\in K$ leads to the situation where both (i) and (ii)
have positive answers, then we conclude that $p_{1}=1$, and move on to find
the next diagonal block of $M$ starting from the index $2$.

###  D.2 Detection of direct product

Let $M$ be a $m\times m$ matrix. Assume that we are given $p,q\in\mathbb{N}$
with $pq=m$. If for some permutation groups $H\leq\mathcal{S}_{p}$ and
$K\leq\mathcal{S}_{q}$, the matrix $M$ lies in $\mathcal{E}(H\otimes K)$
(i.e., $M$ is $H\otimes K$-equivariant), then $M$ can be represented as a
linear combination of Kronecker products:

$\displaystyle M=\sum\limits_{i=1}^{\gamma}(X_{i}\otimes Y_{i})$ (14)

where $X_{i}\in\mathcal{E}(H)$, $Y_{i}\in\mathcal{E}(K)$, and $\gamma$ is a
natural number.

The representation suggests the following strategy of finding a group $G$ of
symmetries of $M$ that is the direct product of two permutation groups on
$[p]$ and $[q]$. First, we express $M$ as a linear combination of Kronecker
products $X_{i}\otimes Y_{i}$ for $X_{i}\in\mathbb{R}^{p\times p}$ and
$Y_{i}\in\mathbb{R}^{q\times q}$. Second, we pick a pair $X_{i}$ and $Y_{i}$
in the linear combination. Third, we call SymFind recursively on $X_{i}$ and
$Y_{i}$ to get group-permutation pairs $(H,\sigma_{H})$ and $(K,\sigma_{K})$.
Finally, we return $(H\otimes K,\sigma_{H}\otimes\sigma_{K})$ where

$\displaystyle\sigma_{H}\otimes\sigma_{K}:[m]\to[m],$
$\displaystyle(\sigma_{H}\otimes\sigma_{K})((a-1)q+b)=(\sigma_{H}(a)-1)q+\sigma_{K}(b)\quad\text{for
all $a\in[p]$ and $b\in[q]$}.$

Our ProdFind is the implementation of the strategy just described. The only
non-trivial steps of the strategy are the first two, namely, to find the
representation of $M$ in (14), and to pick good $X_{i}$ and $Y_{i}$. ProdFind
also has to account for the fact that the representation in (14) holds only
approximately at best in our context.

For the representation finding, ProdFind uses the technique [Van Loan and
Pitsianis, 1993] that solves the following optimisation problem for given
$\gamma$:

$\displaystyle\operatorname*{argmin}_{X_{i}\in\mathbb{R}^{p\times
p},Y_{i}\in\mathbb{R}^{q\times q}}\left\lVert
M-\sum\limits_{i=1}^{\gamma}(X_{i}\otimes Y_{i})\right\rVert_{F}.$ (15)

The technique performs the singular value decomposition of the rearranged
matrix $\hat{M}\in\mathbb{R}^{p^{2}\times q^{2}}$ of $M$ defined by

$\displaystyle\hat{M}_{(i^{\prime}-1)p+i,(j^{\prime}-1)q+j}=M_{(i-1)q+j,(i^{\prime}-1)q+j^{\prime}}$
(16)

for $i,i^{\prime}\in[p]$ and $j,j^{\prime}\in[q]$. To understand why it does
so, note that the optimisation problem in (15) is equivalent to the following
problem:

$\displaystyle\operatorname*{argmin}_{X_{i}\in\mathbb{R}^{p\times
p},Y_{i}\in\mathbb{R}^{q\times
q}}\left\|\hat{M}-\sum\limits_{i=1}^{\gamma}\operatorname*{vec}(X_{i})\operatorname*{vec}(Y_{i})^{T}\right\|_{F}.$
(17)

This new optimisation problem can be solved using SVD. Concretely, if we have
the SVD of $M$, namely,

$\displaystyle\hat{M}=\sum\limits_{i=1}^{\operatorname*{rank}(\hat{M})}s_{i}(u_{i}v_{i}^{T})$
(18)

where $\operatorname*{rank}(\hat{M})$ is the rank of $\hat{M}$, $s_{i}$ is the
$i$-th largest singular value, and $u_{i}$ and $v_{i}$ are the corresponding
left and right singular vectors, we can solve the minimisation problem (17) by
reshaping each term of (18) corresponding to the $\gamma$ largest singular
values, i.e.,

$\displaystyle\hat{M}\approx\sum\limits_{i=1}^{\gamma}\operatorname*{vec}(X_{i})\operatorname*{vec}(Y_{i})^{T}\quad\text{where
}\operatorname*{vec}(X_{i})=\sqrt{s_{i}}u_{i}\text{ and
}\operatorname*{vec}(Y_{i})=\sqrt{s_{i}}v_{i}\text{ for every
$i\in[\gamma]$}.$

The first step of ProdFind is to apply the above technique [Van Loan and
Pitsianis, 1993]. ProdFind performs the SVD of the matrix $\hat{M}$ and gets
$s_{i},u_{i},v_{i}$ in (18) where $s_{1}\geq s_{2}\geq\ldots\geq s_{l}$ for
some $l$. Then, it sets $L=s_{1}/5$, and picks $\gamma$ so that all $i$’s with
$s_{i}\geq L$ are included when we approximate $\hat{M}$.

The second step of ProdFind is to pick $X_{i}$ and $Y_{i}$. ProdFind simply
picks $X_{1}$ and $Y_{1}$ that correspond to matrices (reshaped from $u_{1}$
and $v_{1}$) for the largest singular value $s_{1}$. We empirically observed
that the matrices $X_{1}$ and $Y_{1}$ are most informative about the
symmetries of $M$, and are least polluted by noise from the training of
SATNet.

The third and fourth steps of ProdFind are precisely the last two steps of the
strategy that we described above.

Our description of SymFind in the main text says that ProdFind gets called for
all the divisors $p$ of $n$ with $q$ being $n/p$, since we do not know the
best divisior $p$ a priori, and that each invocation generates a new candidate
in the candidate list $\mathcal{A}$. However, in practice, we ignore some
divisors that violates our predefined criteria. In particular, we set an upper
bound $U$ on $\gamma$, and use a divisor $p$ only when ProdFind for $p$ picks
$\gamma$ with $\gamma\leq U$, i.e., not too many terms are considered in the
linear combination of Kronecker products. This criterion can be justified by
the following lemma:

######  Lemma D.1.

Let $G=H\otimes K$ be the direct product of permutation groups $H$ and $K$ on
$[p]$ and $[q]$, respectively. Consider a $G$-equivariant matrix
$M\in\mathcal{E}(G)$ and its rearranged version $\hat{M}$ in (16). Then,
$\operatorname*{rank}(\hat{M})\leq\min\left(|\mathcal{B}(H)|,\,|\mathcal{B}(K)|\right)$.

###### Proof.

Let $M\in\mathcal{E}(G)$. The matrix $M$ can be expressed as a linear
combination of the basis elements in $\mathcal{B}(G)$:

$\displaystyle M$
$\displaystyle=\sum\limits_{i=1}^{|\mathcal{B}(H)|}\sum\limits_{j=1}^{|\mathcal{B}(K)|}\alpha_{ij}(A_{i}\otimes
B_{j})$
$\displaystyle=\sum\limits_{j=1}^{|\mathcal{B}(K)|}\left(\sum\limits_{i=1}^{|\mathcal{B}(H)|}\alpha_{ij}A_{i}\right)\otimes
B_{j}$ (19)
$\displaystyle=\sum\limits_{i=1}^{|\mathcal{B}(H)|}A_{i}\otimes\left(\sum\limits_{j=1}^{|\mathcal{B}(K)|}\alpha_{ij}B_{j}\right)$
(20)

where $A_{i}\in\mathcal{B}(H)$ and $B_{j}\in\mathcal{B}(K)$. We can also
rewrite (19) and (20) with equations for $\hat{M}$:

$\displaystyle\begin{aligned}
\hat{M}&=\sum\limits_{j=1}^{|\mathcal{B}(K)|}\operatorname*{vec}(A_{j}^{\prime})\operatorname*{vec}(B_{j})^{T}=\sum\limits_{i=1}^{|\mathcal{B}(H)|}\operatorname*{vec}(A_{i})\operatorname*{vec}(B_{i}^{\prime})^{T}\\\
\\\ &\qquad\qquad\qquad\text{where
}A_{j}^{\prime}=\left(\sum\limits_{i=1}^{|\mathcal{B}(H)|}\alpha_{ij}A_{i}\right)\text{
and
}B_{i}^{\prime}=\left(\sum\limits_{j=1}^{|\mathcal{B}(K)|}\alpha_{ij}B_{j}\right).\end{aligned}$

Since the summands
$\operatorname*{vec}(A_{j}^{\prime})\operatorname*{vec}(B_{j})^{T}$ and
$\operatorname*{vec}(A_{i})\operatorname*{vec}(B_{i}^{\prime})^{T}$ are
rank-$1$ matrices, by the subadditivity of rank, i.e.,
$\operatorname*{rank}(X+Y)\leq\operatorname*{rank}(X)+\operatorname*{rank}(Y)$,
we have
$\operatorname*{rank}(\hat{M})\leq\min(|\mathcal{B}(H)|,\,|\mathcal{B}(K)|)$.
∎

Having fewer basis elements is generally better because it leads to a small
number of parameters to learn in SymSATNet. The above lemma says that a
divisor $p$ with large $\gamma$ (which roughly corresponds to the large rank
of $\hat{M}$ in the lemma) leads to large $|\mathcal{B}(H)|$ and
$|\mathcal{B}(K)|$. Our criterion is designed to avoid such undesirable cases.

###  D.3 Detection of wreath product

Assume that the permutation groups $G$ and $H$ act transitively on $[p]$ and
$[q]$ (i.e., for all $i,j\in[p]$ and $i^{\prime},j^{\prime}\in[q]$, there are
$g\in G$ and $h\in H$ such that $g(i)=j$ and $h(i^{\prime})=j^{\prime}$). We
recall that in this case, every equivariant matrix $M$ under the wreath
product of groups $H\wr G$ can be expressed as follows [Wang et al., 2020]:

$\displaystyle M=A\otimes\mathbf{1}_{q}+I_{p}\otimes B$ (21)

for some $A\in\mathcal{E}(G)$ and $B\in\mathcal{E}(H)$, where
$\mathbf{1}_{m}$, $I_{m}$ are everywhere-one, identity matrices in
$\mathbb{R}^{m\times m}$. We use this general form of $\mathcal{B}(H\wr
G)$-equivariant matrices and make ProdFind detect the case that symmetries are
captured by a wreath product.

Our change in ProdFind is based on a simple observation that the form in (21)
is exactly the one in (14) with $\gamma=2$. We change ProdFind such that if we
get $\gamma=2$ while runnning $\textsc{ProdFind}(M,p)$, we check whether $M$
can be written as the form in (21). This checking is as follows. Using given
$p$ and $q$, we create blocks

$\displaystyle M^{(i,j)}=M[(i-1)\times q+1:i\times q,(j-1)\times q+1:j\times
q]$

of $q\times q$ submatrices in $M$ for all $i,j\in[p]$. Then, we test whether
all $M^{(i,i)}-M^{(j,j)}$ and $M^{(i,j)}$ for $i\neq j$ are constant matrices
(i.e., matrices of the form $\alpha\mathbf{1}_{q}$ for some
$\alpha\in\mathbb{R}$). If this test passes, $M$ has the desired form. In that
case, we compute $A$ and $B$ from $M$, and recursively call SymFind on $A$ and
$B$.

We can also make $\textsc{SumFind}(M)$ detect wreath product. The required
change is to perform a similar test on each diagonal block of the final
$M_{\sigma^{-1}}$ computed by $\textsc{SumFind}(M)$. The only difference is
that when the test fails, we look for a rearrangement of the rows and columns
of the diagonal block that makes the test succeed. If the test on a diagonal
block succeeds (after an appropriate rearrangement), the group corresponding
to this block, which is one summand of a direct sum, becomes a wreath product.

##  Appendix E Computational Complexity of SymFind

In this section, we analyse the computational complexity of our SymFind
algorithm in terms of the dimension of the input matrix. Suppose that an
$n\times n$ matrix is given to SymFind as an input. The complexity of SymFind
shown in Algorithm 1 is $O(n^{3+\epsilon})$ for any arbitrarily small
$\epsilon>0$.

We can show the claimed complexity of SymFind using induction. When SumFind is
called in the line 9 of Algorithm 1, it first spends $O(n^{3})$ time for
clustering, and then makes recursive calls to SymFind with $n_{i}\times n_{i}$
submatrices for $\sum n_{i}=n$. Each of these calls takes
$O(n_{i}^{3+\epsilon})$ time by induction, and the total cost of all the calls
is $\sum_{i}O(n_{i}^{3+\epsilon})=O(n^{3+\epsilon})$. Also, when ProdFind is
called in the lines 11-12 of Algorithm 1, for each divisor $p$ of $n$,
ProdFind performs SVD in $O(n^{3})$ steps, and makes recursive calls to
SymFind. The recursive calls here together cost
$O(p^{3+\epsilon/2}+(n/p)^{3+\epsilon/2})=O(n^{3+\epsilon/2})$ by induction.
Thus, the loop in lines 11-14 costs
$O(n^{3+\epsilon/2}d(n))=O(n^{3+\epsilon})$, where $d(n)$ is the number of
divisors of $n$, since $d(n)=O(n^{\epsilon/2})$ for any arbitrarily small
$\epsilon>0$. It remains to show that applying the Reynolds operator (in the
lines 2 and 6) and counting the number of basis elements (in the line 15) take
$O(n^{3+\epsilon})$ steps. The former costs $O(n^{2})$ since it is actually an
average pooling operation when we consider permutation groups. The latter
costs at most $O(n)$ (when it is given
$\mathcal{B}(\bigoplus_{i=1}^{n}\mathbb{Z}_{1})$).

##  Appendix F Symmetries in Sudoku and Rubik’s Cube Problems

In this section, we describe the group symmetries in Sudoku and the completion
problem of Rubik’s cube. In both problems, we can find certain group
$G\leq\mathcal{S}_{n}$ acting on $\mathbb{R}^{k\times n}$ such that for any
valid assignment $V\in\mathbb{R}^{k\times n}$ of problem, $g\cdot V$ is also
valid for any $g\in G$.

In $9\times 9$ Sudoku, we denote the first three rows of a Sudoku board by the
first band, and the next three rows by the second band, and the last three
rows by the third band. In the same manner, we denote each three columns by a
stack. For Sudoku problem, one can observe that any row permutations within
each band preserve the validity of the solutions. Also, the permutations of
bands preserve the validity. We can represent this type of hierarchical
actions by wreath product groups. In this case, 3 bands are permuted in a
higher level, and 3 rows in each band are permuted in a lower level, which
corresponds to the group action of $\mathcal{S}_{3}\wr\mathcal{S}_{3}$. The
same process can be applied to the stacks and the columns, which also
corresponds to the group $\mathcal{S}_{3}\wr\mathcal{S}_{3}$. Furthermore, the
permutations of number occurences in Sudoku (e.g., switching all the
occurences of 3’s and 9’s) also preserves the validity of the solutions. In
this case, any permutations over $[9]$ are allowed, thus $\mathcal{S}_{9}$
represents this permutation action. Overall, we have three permutation groups,
$G_{1}=\mathcal{S}_{3}\wr\mathcal{S}_{3}$,
$G_{2}=\mathcal{S}_{3}\wr\mathcal{S}_{3}$, and $G_{3}=\mathcal{S}_{9}$. These
three groups are in different levels; $G_{1}$ is at the outermost level, and
$G_{3}$ is at the innermost level. Note that these three types of actions are
commutative, i.e., the order does not matter if we apply the actions in
different levels. This relation forms a direct product between each level, and
we conclude
$G=(\mathcal{S}_{3}\wr\mathcal{S}_{3})\otimes(\mathcal{S}_{3}\wr\mathcal{S}_{3})\otimes\mathcal{S}_{9}$.

For Rubik’s cube, we first introduce some of the conventional notations. We
are given a $3\times 3\times 3$ Rubik’s cube composed of 6 faces and 9
facelets in each face. There are also 26 cubies in the Rubik’s cube, and they
are classified into 8 corners, 12 edges, and 6 centers. To change the color
state of Rubik’s cube, we can move the cubies with 9 types of rotations, which
are denoted by $U$, $D$, $F$, $B$, $L$, $R$, $M$, $E$, $S$. We can represent
each color state of the Rubik’s cube by a function from $F=[6]\times[9]$ to
$C=[6]$, where $F$ encodes the face and facelets, and $C$ encodes the colors.
The initial state of Rubik’s cube is $s:F\to C$ satisfying $s(i,j)=i$ for all
$(i,\,j)\in F$. If there exists a sequence of rotations that transforms $s$ to
the initial state, then $s$ is solvable. Our completion problem of Rubik’s
cube is to find a color assignment such that the assignment is solvable. For
example, if a corner cubie contains 2 same colors, it is not solvable because
these adjacent same colors cannot be separated by the 9 types of rotations.

We can observe the following group symmetries in the completion problem of
Rubik’s cube: any solvable color state is still solvable after being
transformed by the 9 types of rotations. These 9 rotations generate a
permutation group $\mathcal{R}_{54}$ acting on the $6\times 9$ facelets of a
Rubik’s cube, which is called the Rubik’s cube group. Furthermore, like
Sudoku, we can also consider the permutations of color occurences. In this
case, if we assume the colors $0$ and $5$ are initially on the opposite sides,
then the solvability can be broken by switching the colors $0$ and $1$ since
the colors $0$ and $5$ are then no longer on the opposite sides. Considering
this, computing all possible permutations acting on the color occurences is
not trivial. Instead, we can generate such group $\mathcal{R}_{6}$ by 3
elements, each of which corresponds to the $90^{\circ}$ rotation in one of the
three axes in 3-dimension. Finally, these two groups $\mathcal{R}_{54}$ and
$\mathcal{R}_{6}$ form different levels of actions, and actions from different
levels commute each other. Therefore, we conclude
$G=\mathcal{R}_{54}\otimes\mathcal{R}_{6}$.

##  Appendix G Efficient Implementation of SymSATNet

Algorithm 2 Forward pass of SATNet

1: Input: $V_{\mathcal{I}}$

2: init random unit vectors $v_{o}$, $\forall o\in\mathcal{O}$

3: compute $\Omega=VS^{T}$

4: repeat

5: for $o\in\mathcal{O}$ do

6: compute $g_{o}=\Omega s_{o}-\|s_{o}\|^{2}v_{o}$

7: compute $v_{o}=-g_{o}/\|g_{o}\|$

8: update $\Omega=\Omega+(v_{o}-v_{o}^{\text{prev}})s_{o}^{T}$

9: end for

10: until not converged

11: Output: $V_{\mathcal{O}}$

Algorithm 3 Forward pass of SymSATNet

1: Input: $V_{\mathcal{I}}$

2: init random unit vectors $v_{o}$, $\forall o\in\mathcal{O}$

3: repeat

4: for $o\in\mathcal{O}$ do

5: compute $g_{o}=Vc_{o}-C_{o,o}v_{o}$

6: compute $v_{o}=-g_{o}/\|g_{o}\|$

7: update $V=V+(v_{o}-v_{o}^{\text{prev}})\mathds{1}_{o}^{T}$

8: end for

9: until not converged

10: Output: $V_{\mathcal{O}}$

Our SymSATNet includes forward-pass and backward-pass algorithms as described
in Section 3. Here, we describe slight changes of the forward-pass and
backward-pass algorithms in SymSATNet, and the improvement of efficiency
obtained by these changes. Algorithms 2 and 3 are the forward pass algorithms
of SATNet and SymSATNet, and Algorithms 4 and 5 are the backward pass
algorithms of SATNet and SymSATNet. In those algorithms, we let $c_{o}$ be the
$o$-th column of the matrix $C$, and $v_{o}^{\text{prev}}$ and
$u_{o}^{\text{prev}}$ be the previous $o$-th columns of $V$ and $U$ before the
update. Also, let $P_{o}=I_{k}-v_{o}v_{o}^{T}$, and $\mathds{1}_{o}$ be the
n-dimensional one-hot vector whose only $o$-th element is $1$. Note that our
algorithms are only slightly different from the ones of original SATNet, and
the only difference is that the inner products $s_{i}^{T}s_{j}$ are
substituted by $C_{i,j}$ which can be directly derived from $C=S^{T}S$. This
allows us to implement the forward-pass and the backward-pass algorithm of
SymSATNet in the same manner as the original SATNet.

The above changes bring a small difference in the computational complexity. In
SATNet, new matrices $\Omega$ and $\Psi$ are required for the rank-$1$ update
in each loop. Recall that $V,U\in\mathbb{R}^{k\times n}$ and
$S\in\mathbb{R}^{m\times n}$ imply $\Omega,\Psi\in\mathbb{R}^{k\times m}$.
SATNet costs $O(nmk)$ twice every iteration, in the lines 6, 8 of Algorithm 2
and in the lines 5, 7 of Algorithm 4. SymSATNet costs $O(n^{2}k)$ once every
iteration, in the line 5 of Algorithm 3 and in the line 5 of Algorithm 5.
Another slight difference is in the first line of the body of each for loop,
where SATNet computes $\|s_{o}\|^{2}$, but SymSATNet uses the constant
$C_{o,o}$. If we denote the required number of iterations by $t$, then SATNet
costs $O(nmk\cdot 2t)$, and SymSATNet costs $O(n^{2}k\cdot t)$ inside the
loop. These make one of the major differences between the runtime of SATNet
and SymSATNet, since these are scaled by $t$, the number of iterations until
convergence. Outside the loop, SATNet additionally costs $O(nmk)$ in the line
3 of Algorithm 2, and SymSATNet costs $O(n^{2}d)$ in the forward computation
of $C$ in (4), and in the backward computation with the chain rule in (5).
These costs are not significant because each of these occurs only once every
gradient update in SATNet and SymSATNet.

We also inspected the values of $t$ in SATNet and SymSATNet in a training run
for Sudoku and Rubik’s cube problem in the configuration we used in our
experiments. For Sudoku, SATNet repeated the loop $t=21.39$ times on average,
in the range of $15\leq t\leq 32$, while SymSATNet finished the loop in
$t=19.10$ on average, in the scope of $3\leq t\leq 26$. For Rubik’s cube,
SATNet completes the loop in $t=9.25$ on average, whose scope was $3\leq t\leq
25$, while SymSATNet iterates $t=7.30$ times on average, with the range $2\leq
t\leq 26$. These results show that the component $t$ is reduced in SymSATNet
for those two problems, and bring further improvement of efficiency in
practice.

##  Appendix H Hyperparameters

In this section, we specify the hyperparameters to run SymFind algorithm and
the validation steps in SymSATNet-Auto. Each validation step requires a
threshold of improvement of validation accuracy to measure the usefulness of
each smaller part in a given group, i.e., only the smaller parts showing
improvement greater than the threshold were considered to be useful, and were
combined to construct the whole group. We determined this threshold by the
number of corruption in the dataset; $0.05$ for the dataset with 0 and 1
corruption, $0.15$ for 2 corruption, and $0.2$ for 3 corruption. Also, our
SymFind algorithm receives $\lambda$ as an input, which determines the degree
of sensitivity of SymFind. In practice, this tolerance was implemented by two
hyperparameters, say $\lambda_{1}$ and $\lambda_{2}$. $\lambda_{1}$ was used
to compute a threshold to decide whether a pair of entries have the same value
within our tolerance so that they have to be clustered. We computed this
threshold by the multiplication of $\lambda_{1}$ and the norm of input matrix,
and $\lambda_{1}$ was fixed by $0.1$ for both problems. $\lambda_{2}$ played a
role of a threshold to conclude that the input matrix truly has the symmetries
under the discovered group by SymFind. $\lambda_{2}$ was used in the last step
of SymFind, where the largest group is selected among the candidates. With the
Reynolds operator, we computed the distance between the input matrix and the
equivariant space under each candidate group, and filtered out the groups
whose distance is greater than $\lambda_{2}$. We determined $\lambda_{2}$ by
the number of corruption in the dataset; $0.4$ for the dataset with 0
corruption, $0.5$ for 1 corruption, $0.55$ for 2 corruption, and $0.6$ for 3
corruption.

Algorithm 4 Backward pass of SATNet

1: Input: $\\{\partial\ell/\partial v_{o}:o\in\mathcal{O}\\}$

2: init $U_{\mathcal{O}}=0$ and $\Psi=U_{\mathcal{O}}S_{\mathcal{O}}^{T}=0$

3: repeat

4: for $o\in\mathcal{O}$ do

5: compute $dg_{o}=\Psi s_{o}-\|s_{o}\|^{2}u_{o}-\partial\ell/\partial v_{o}$

6: compute $u_{o}=-P_{o}dg_{o}/\|g_{o}\|$

7: update $\Psi=\Psi+(u_{o}-u_{o}^{\text{prev}})s_{o}^{T}$

8: end for

9: until not converged

10: Output: $U_{\mathcal{O}}$

Algorithm 5 Backward pass of SymSATNet

1: Input: $\\{\partial\ell/\partial v_{o}:o\in\mathcal{O}\\}$

2: init $U_{\mathcal{O}}=0$

3: repeat

4: for $o\in\mathcal{O}$ do

5: compute $dg_{o}=Uc_{o}-C_{o,o}u_{o}-\partial\ell/\partial v_{o}$

6: compute $u_{o}=-P_{o}dg_{o}/\|g_{o}\|$

7: update $U=U+(u_{o}-u_{o}^{\text{prev}})\mathds{1}_{o}^{T}$

8: end for

9: until not converged

10: Output: $U_{\mathcal{O}}$

##  Appendix I GPU Usage

The implementation of SymSATNet is based on the original SATNet code, and thus
the calculations in the forward-pass and the backward-pass algorithms of
SymSATNet can be accelerated by GPU. We used GeForce RTX 2080 Ti for every
running of SATNet and SymSATNet.

##  Appendix J Ablation Study of Validation Steps

Figure 6: Best test accuracies of SymSATNet-Auto with or without the
validation step using noisy Rubik’s cube datasets.

In this section, we present the results of ablation study for the additional
validation step in SymSATNet-Auto. All the configurations are the same, except
that the discovered group symmetries by SymFind are directly exploited by
SymSATNet without the subsequent validation step. We used noisy Sudoku and
Rubik’s cube datasets, and repeated our experiment with each dataset ten
times. Here we report the average test accuracies with $95\%$ confidence
interval.

In Sudoku, our SymFind algorithm always succeeded in finding the full group
symmetries, and these correct symmetries were preserved by the validation step
since they always significantly improved the validation accuracies after
projecting the parameter $C$ of SATNet. For this reason, omitting the
validation step did not change the performance of SymSATNet-Auto. Thus, we do
not plot the results for the Sudoku case. For Rubik’s cube, the results are
shown in Figure 6. Without the validation step, SymSATNet-Auto often overfit
into the wrong group symmetries induced by the noise in the parameter matrix
$C$ of SATNet, resulting in lower improvement of performance. These results
show that the additional validation step is useful not only when the training
and the test examples share the same distribution, but also when the
distribution of training examples is the perturbation of that of test examples
by noise.

##  Appendix K Emergence of Symmetries in SATNet

(a) Sudoku

(b) Rubik’s cube

Figure 7: Projection errors of the SATNet’s parameter $C$ during 100 epochs.
We repeated 30 training trials for each problem to report the projection
errors.

In this section, we report our experiment for detecting the appearance and
disappearance of symmetries in the SATNet’s parameter matrix $C$ throughout
the training epochs. We measured how close $C$ is to $\mathcal{E}(G)$, the
space of equivariant matrices under $G$, by computing
$\|\mathrm{prj}(G,C)-C\|_{F}$ (denoted by projection error below) where $G$ is
the ground-truth permutation group for a given learning problem (e.g., Sudoku
or Rubik’s cube). We evaluated the projection errors throughout 100 training
epochs. Our experiment for each problem was repeated 30 times.

Figure 7 shows the results. For Sudoku, in multiple cases out of 30 trials,
the projection error hit the lowest point around the 5-10th epochs, and after
these epochs, the projection error started to increase until certain epochs,
so that the training ended with a high projection error. We did not observe
any clear difference in the training loss (or the test loss) between these
cases with high projection errors and the other cases (which showed low
projection errors). This result suggests possible overfitting from the
perspective of symmetry discovery in the original SATNet. For Rubik’s cube, in
all of our trials, the projection error always decreased throughout the
epochs, and no sign of overfitting was detected. Answering why this is the
case is an interesting topic for future research. Also, in the very beginning
of the training for both problems, the projection errors were not sufficiently
small. By choosing proper stopping epochs in SymSATNet-Auto, we can avoid high
projection errors, as we did (i.e., we picked the 10th epoch for Sudoku, and
the 20th epoch for Rubik’s cube). Finally, we point out that there may be
factors other than the projection error that influence the performance of our
symmetry-discovery algorithm. If found, those factors would become useful
tools for analysing the theoretical guarantees of SymFind algorithm, and we
leave it to future work.

##  Appendix L SymSATNet with Auxiliary Variables

Table 3: Best train and test accuracies during 100 epochs and average total train times ($10^{2}$ sec). Each experiment was repeated 10 times and its average and 95% confidence interval are reported. Additional times for automatic symmetry detection (SymFind and validation) are reported after $+$. Model | Sudoku | Cube
---|---|---
| Train Acc. | Test Acc. | Time | Train Acc. | Test Acc. | Time
SATNet-Plain | 93.5% | 88.1% | 48.0 | 99.4% | 55.7% | 1.8
|  $\pm$1.2% |  $\pm$1.3% |  $\pm$0.17 |  $\pm$0.1% |  $\pm$0.5% |  $\pm$0.01
SATNet-300aux | 99.8% | 97.9% | 90.3 | 99.8% | 56.5% | 14.0
|  $\pm$0.0% |  $\pm$0.2% |  $\pm$0.68 |  $\pm$0.0% |  $\pm$0.6% |  $\pm$0.12
SymSATNet-Plain | 98.8% | 99.2% | 25.6 | 67.1% | 66.9% | 1.1
|  $\pm$0.1% |  $\pm$0.2% |  $\pm$0.14 |  $\pm$0.2% |  $\pm$0.9% |  $\pm$0.00
SymSATNet-300aux | 99.2% | 99.3% | 53.6 | 69.6% | 67.6% | 9.5
|  $\pm$0.1% |  $\pm$0.1% |  $\pm$0.56 |  $\pm$0.4% |  $\pm$0.5% |  $\pm$0.29
SymSATNet-Auto | 99.3% | 99.5% | 22.7 | 70.2% | 68.1% | 3.4
|  $\pm$0.1% |  $\pm$0.2% |  $+$0.14 $\pm$0.35 |  $\pm$3.1% |  $\pm$2.0% |  $+$0.66 $\pm$0.19

(a) Sudoku

(b) Rubik’s cube

Figure 8: Test accuracies throughout the first 100 training epochs. Each
training run was repeated 10 times to compute the average and the 95%
confidence interval.

By introducing auxiliary variables, the original SATNet can achieve higher
expressive power and better performance while trading-off the runtime
efficiency. The similar process can be done in SymSATNet. If a permutation
group $G$ captures group symmetries of a problem to be solved without
auxiliary variables, then $G\oplus\mathcal{I}_{m}$ represents the extension of
the same symmetries with the $m$ auxiliary variables. The singleton group
$\mathcal{I}_{m}$ here acts trivially (i.e., permutes nothing) on those
auxiliary variables.

We now report the findings of our experiments with SymSATNet with 300
auxiliary variables (SymSATNet-300aux) that retains the above extension of
group symmetries. We compared the performance of SymSATNet-300aux with the
other four models that we used before, on the 0-corrupted Sudoku and Rubik’s
cube datasets. We denote SymSATNet without auxiliary variables by SymSATNet-
Plain to distinguish it from SymSATNet-300aux. We repeated the training for
each dataset 10 times. Here we report the average and best test accuracies,
the average training times, and $95\%$ confidence interval. We trained
SymSATNet-300aux with the learning rate $\eta=4\times 10^{-2}$, which is the
same as SymSATNet-Plain. All the experimental setups were the same as before.

Figure 8 and Table 3 show the results. In both tasks, SymSATNet-300aux
outperformed all the other models except SymSATNet-Auto. Also, for Rubik’s
cube, SymSATNet-300aux showed the fastest convergence in epoch among all the
models. Note that both the test accuracies and the training times of
SymSATNet-300aux were remarkably improved over those of SATNet baselines.
These results demonstrate that SATNet with auxiliary variables still enjoyed
the benefits by exploiting a minor extension of group symmetries, although the
benefits were not as significant as in the case of SymSATNet-Plain.

##  Appendix M Loss Curves

(a) Training losses in Sudoku with 0 corruption

(b) Test losses in Sudoku with 0 corruption

(c) Training losses in Sudoku with 1 corruption

(d) Test losses in Sudoku with 1 corruption

(e) Training losses in Sudoku with 2 corruptions

(f) Test losses in Sudoku with 2 corruptions

(g) Training losses in Sudoku with 3 corruptions

(h) Test losses in Sudoku with 3 corruptions

Figure 9: Training and test loss curves for Sudoku. Each loss is averaged over
10 trials and each $95\%$ confidence interval is included.

(a) Training losses in Rubik’s cube with 0 corruption

(b) Test losses in Rubik’s cube with 0 corruption

(c) Training losses in Rubik’s cube with 1 corruption

(d) Test losses in Rubik’s cube with 1 corruption

(e) Training losses in Rubik’s cube with 2 corruptions

(f) Test losses in Rubik’s cube with 2 corruptions

(g) Training losses in Rubik’s cube with 3 corruptions

(h) Test losses in Rubik’s cube with 3 corruptions

Figure 10: Training and test loss curves for the Rubik’s cube problem. Each
loss is averaged over 10 trials and each $95\%$ confidence interval is
included.