# Chemo-dynamical substructure in the M31 inner halo globular clusters:
Further evidence for a recent accretion event

Geraint F. Lewis1, Brendon J. Brewer2, Dougal Mackey3, Annette M. N.
Ferguson4, Yuan (Cher) Li2 & Tim Adams1
1Sydney Institute for Astronomy, School of Physics, A28, The University of
Sydney, NSW 2006, Australia
2Department of Statistics, The University of Auckland Private Bag 92019,
Auckland 1142, New Zealand
3Research School of Astronomy and Astrophysics, Australian National
University, Canberra, ACT 2611, Australia
4Institute for Astronomy, University of Edinburgh Royal Observatory, Blackford
Hill, Edinburgh EH9 3HJ, UK E-mail<EMAIL_ADDRESS>(GFL)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Based upon a metallicity selection, we identify a significant sub-population
of the inner halo globular clusters in the Andromeda Galaxy which we name the
Dulais Structure. It is distinguished as a co-rotating group of 10-20 globular
clusters which appear to be kinematically distinct from, and on average more
metal-poor than, the majority of the inner halo population. Intriguingly, the
orbital axis of this Dulais Structure is closely aligned with that of the
younger accretion event recently identified using a sub-population of globular
clusters in the outer halo of Andromeda, and this is strongly suggestive of a
causal relationship between the two. If this connection is confirmed, a
natural explanation for the kinematics of the globular clusters in the Dulais
Structure is that they trace the accretion of a substantial progenitor ($\sim
10^{11}M_{\odot}$) into the halo of Andromeda during the last few billion
years, that may have occurred as part of a larger group infall.

###### keywords:

Local Group $--$ globular clusters $--$ Andromeda galaxy

††pubyear: 2022††pagerange: Chemo-dynamical substructure in the M31 inner halo
globular clusters: Further evidence for a recent accretion event–A

##  1 Introduction

Within our cosmological paradigm, large galaxies like the Milky Way have grown
over time through the accretion of smaller stellar systems. The rate of this
accretion is expected to be a relatively steady rain over cosmic time,
potentially punctuated by more substantial accretion events (Wechsler et al.,
2002). Since the discovery of the disrupting Sagittarius Dwarf Galaxy (Ibata
et al., 1994), extensive surveys of the halo of the Milky Way have revealed a
wealth of substructure in the form of almost one hundred stellar streams
(e.g., Ibata et al., 2020; Li et al., 2022; Mateu, 2022), as well as more
ancient extended structures (e.g., Myeong et al., 2018; Belokurov et al.,
2018; Helmi et al., 2018; Kruijssen et al., 2020). Furthermore, deep imaging
of the Andromeda Galaxy (M31) has uncovered large-scale stellar substructure
strewn throughout the halo, similarly revealing an extensive history of
accretion (e.g., McConnachie et al., 2018). Recently, the identification of
kinematically distinct populations of globular clusters orbiting in the outer
halo of the Andromeda Galaxy, distinguished by their relationship to various
underlying stellar substructures, revealed two pronounced accretion episodes,
one ancient and the other more recent (Mackey et al., 2019b).

The discovery of these distinct subgroups in the outer halo globular cluster
population motivates us to revisit the kinematic properties of the inner
globular clusters of Andromeda. Being the nearest large galaxy to our own,
there has been significant interest in similarities and differences between
the Milky Way and Andromeda globular cluster populations (e.g. van den Bergh,
1969; Huchra et al., 1982; Hartwick & Sargent, 1974; Freeman, 1983; Burstein
et al., 1984; Crampton et al., 1985; Elson & Walterbos, 1988), identifying a
bimodel metallicity distribution (e.g. Barmby et al., 2000), and potential
small-scale kinematic and chemical substructures (Perrett et al., 2002;
Perrett et al., 2003). Recently, Caldwell & Romanowsky (2016) undertook an
extensive spatial-kinematic-chemical analysis of the inner population of
Andromeda’s globular clusters based upon a homogeneous spectroscopic data set.
Splitting their sample into a metal-rich, intermediate- and poor- sub-
populations, they conclude that the kinematics of the metal-rich components
are consistent with galactic rotation, whilst more metal-poor sub-components
show a higher rotation and velocity dispersion and are more weakly
concentrated with regards to galactic disk, and hence appear to represent a
transition to the larger halo population of globular clusters. However, we
note that the kinematic analysis of Caldwell & Romanowsky (2016) was
undertaken with respect to rotation about the minor axis of Andromeda, without
exploration of the potential orientation of the rotation.

Hence, here we undertake a more expansive exploration of the kinematic and
chemical distributions of the inner globular cluster populations in Andromeda,
including additional freedom to explore the axes of rotation. As detailed
below, the results of this analysis is the identification of a pronounced
substructure within the inner globular cluster population, which we name the
Dulais Structure111Pronounced ‘dill’+‘ice’, after the Welsh river name meaning
‘black stream’.; intriguingly the kinematic and spatial properties causally
link it to the recent accretion events in the outer halo identified by Mackey
et al. (2019b). The layout of this paper is as follows: In Section 2 we
outline our approach, including the data employed and the statistical methods
used in model assessment. We discuss the results of this study in 3 and
present our conclusions in Section 4.

##  2 Approach

In our previous exploration of the kinematic properties of the globular
clusters in the outer halo of M31 (i.e., at projected galactocentric radii
$>25$ kpc), we defined two distinct populations according to whether or not
the clusters were associated with underlying stellar substructure (Mackey et
al., 2019a). This indicated whether the globular clusters were accreted
relatively recently, such that their progenitors are still undergoing tidal
disruption, or whether they potentially trace more ancient accretions which
have been completely dispersed. However, it is impossible to identify faint
stellar substructures against the stellar disk and extensive tidal debris
within the more central regions of M31 (McConnachie et al., 2018), meaning
that this experiment cannot easily be repeated for the inner halo population.
Instead, our present study focuses upon splitting the inner clusters according
to their chemical enrichment, or metallicity, in order to assess whether there
are kinematically distinct sub-groups inhabiting the inner regions of
Andromeda.

We adopt a Bayesian modelling approach, providing an assessment of Bayesian
evidence (marginal likelihood) and allowing a determination of the relative
efficacy of various models in fitting the data. The goal of our analysis is to
compare whether a single component model, or a dual component model separating
the globular clusters based upon their metallicity, provides a more plausible
explanation of the data.

As detailed below, we assume that the kinematics of a given globular cluster
population can be represented by a rotational component coupled with a
velocity dispersion, and consider three distinct forms for the rotation. In
the first stage of the analysis, we fit models to the total globular cluster
population, calculating the Bayesian evidence for each rotational form. We
then repeated the process to explore whether two distinct populations, split
by metallicity and each with their own rotational component and velocity
dispersion, represented a better fit of the data. Given the relatively large
uncertainties on the available metallicity measurements, we adopted a
probabilistic approach to the two-population "split", weighting each globular
cluster’s contribution to a kinematic component according to its metallicity
and associated uncertainty relative to a metallicity threshold (treated as a
free parameter). For the sake of clarity, the following discussion refers to
globular clusters with metallicites lower than the metallicity cut as metal-
poor, whilst those above are referred to as metal-rich.

###  2.1 Globular Cluster Data

The data underpinning our analysis is the publicly-available catalogue of
Caldwell & Romanowsky (2016), originally obtained through spectroscopic
observations for nearly all of the known globular clusters in the inner
regions of Andromeda using the Hectospec multi-fibre spectrograph (Fabricant
et al., 2005) on the 6.5m MMT Observatory telescope in Arizona. Caldwell &
Romanowsky (2016) used these spectra to provide a set of uniformly-derived
line-of-sight velocities and metallicity estimates for approximately $94\%$ of
globular clusters with projected galactocentric radii inside $21$ kpc. Our
previous analysis of globular clusters in Andromeda’s outer halo (Mackey et
al., 2019a, b) had an inner cut-off of $25$ kpc, so we searched the
literature, and in particular the cluster catalogue from the Pan-Andromeda
Archaeological Survey (PAndAS; Veljanoski et al., 2014), to ensure we were not
missing any objects in the range $21-25$ kpc. We also updated the velocities
from Caldwell & Romanowsky (2016) with more precise PAndAS measurements for a
few clusters (G268, H14, and SK255B). In summary, the total base catalogue
comprises 345 globular clusters extending to a projected galactocentric radius
of $25$ kpc.

Caldwell & Romanowsky (2016) convincingly established that the most metal-rich
globular clusters in the central regions of M31 are kinematically and
spatially associated with its disk. Since our analysis is focused on
Andromeda’s inner halo, we removed this component by excluding all clusters
with metallicity $FeH\geq-0.4$222Throughout we use $FeH$ to represent
$[\textrm{Fe}/\textrm{H}]$.. We further removed all globular clusters for
which the observations were of insufficient quality to determine a reliable
metallicity estimate.

Models $(\rm{I}\ \&\ \rm{II})$ | Functional form
---|---
$V$ | $A_{k}\sin(\theta-\phi_{k})$
$S$ | $A_{k}\left(x\sin\phi_{k}-y\cos\phi_{k}\right)$
$F$ | $A_{k}\tanh\left((x\sin\phi_{k}-y\cos\phi_{k})/L_{k}\right)$
$\top$ | Disjunction of all of the above $(\rm{I}\ \&\ \rm{II})$
Table 1: The three different rotational model components considered in our
kinematic models, providing the line-of-sight velocity as a function of
position on the sky. The positions are expressed in either Cartesian
coordinates $(x,y)$ or plane polar coordinates $(r,\theta)$ defined as
specified in the text. A single roman numeral subscript (e.g., $V_{\rm I}$)
corresponds to a single-population model whose parameters are represented with
a single subscript $k=0$ (e.g., $A_{0}$). For models with two populations,
(e.g., $V_{\rm II}$), the globular cluster subgroup below the metallicity cut
is represented by subscript $k=0$, whereas the globular cluster subgroup above
the metallicity cut is represented by a subscript $k=1$. The model $\top$ is
the proposition that one of the three rotational forms is correct.

The final sample for consideration contains $n=278$ globular clusters. Each
cluster possesses the following data:
$(x_{i},y_{i},v_{i},\sigma_{v_{i}},\widehat{FeH_{i}},\sigma_{FeH_{i}})$ –
namely the Cartesian tangent-plane position centred on Andromeda,
$(x_{i},y_{i})$, the Andromeda-centric line-of-sight velocity and associated
uncertainty, $(v_{i},\sigma_{v_{i}})$, and the reported metallicity and
associated uncertainty, $(\widehat{FeH_{i}},\sigma_{FeH_{i}})$; the
metallicity and velocities of the globular cluster population, coupled with
their uncertainties, is presented graphically in the Appendix as Figure 7.
Cartesian positions correspond to tangent-plane distances centred on M31;
these coordinates, and the Andromeda-centric velocities, were calculated
according to the prescription outlined in Mackey et al. (2019a). Note that the
transformation of Cartesian coordinates into regular polar coordinates is
defined as:

$r_{i}=\sqrt{x_{i}^{2}+y_{i}^{2}}\ \ \ \ \ \ \ \ \ \ \ \
\theta_{i}=\arctan{2(y_{i},x_{i})}.$ (1)

Here, $\theta_{i}$ is not the position angle $(PA)$ on the sky that is usually
quoted in astronomy, which is measured from North through East, but instead is
defined from East through North. For convenience, the results of this study
are converted into astronomical position angles at the completion of the
analysis.

The typical individual $1\sigma$ uncertainty on $\widehat{FeH_{i}}$ is
$\sigma_{FeH_{i}}=0.2$ dex, while that on the Andromeda-centric velocity
$v_{i}$ is $\sigma_{v_{i}}=7.5$ km/s. This latter value is dominated by the
measurement uncertainty (typically $\approx 6$ km/s), with a small
contribution due to the assumed systemic motion of M31 (Mackey et al., 2019a).

Parameter | Prior
---|---
$FeH_{\rm cut}$ | Uniform(-3, -1)
$A_{k}$ | Uniform(0, 800) km/s
$\phi_{k}$ | Uniform$(-\pi,\pi)$ radians
$L_{l}$ | Uniform$(0,2)$ degrees
$\Sigma_{k0}$ | Uniform(0, 400) km/s
$s_{k}$ | Uniform(-4, 4) (degrees)-1
Table 2: Prior distributions for the unknown parameters in our models. The $L$
parameter only applies to the model with asymptotically flat rotation ($F$).

###  2.2 Kinematic Modelling

As with previous studies (c.f. Mackey et al., 2019b), we explored three
different kinematic models. Each of these is a nonlinear regression model that
predicts line-of-sight velocity $v_{i}$ as a function of position
$(x_{i},y_{i})$ or $(r_{i},\theta_{i})$, and potentially (when there are two
populations) the metallicity $FeH_{i}$. The three models are denoted $V,S,F$
and possess different forms for the rotational component as presented in Table
1; note that these models will be subscripted with ${\rm I}$ and ${\rm II}$ if
they represent models with a single or double kinematic components. $V$ is a
simple sinusoidal function previously used in a variety of studies of the
rotation of the globular cluster population in the halo of Andromeda (e.g.,
Veljanoski et al., 2013, 2014; Caldwell & Romanowsky, 2016), whereas $S$
corresponds to a solid-body rotation. Finally, $F$, represents an
asymptotically-flat rotation curve which smoothly transitions from one side of
the galaxy to the other. For a single-population model (e.g. $V_{\rm I}$), the
index $k$ is fixed at $0$. When the metallicity cut is applied and the model
as two kinematic components (e.g. $V_{\rm II}$), the globular cluster
population below the metallicity cut is represented by subscript $0$, whereas
the globular cluster population above the cut is represented by subscript $1$.

As well as rotation, it is also assumed that each globular cluster population
possesses a velocity dispersion. We adopt the functional form of the velocity
dispersion over the populations to be

$\Sigma_{k}(r)=\Sigma_{k0}\exp{\left(s_{k}r\right)}$ (2)

where $\Sigma_{k0}$ and $s_{k}$ are parameters to be determined; clearly if
$s_{k}=0$ then the velocity dispersion is spatially constant over the
population.

###  2.3 Metallicity Selection

The models considering two populations assume that there exists a metallicity
threshold $FeH_{\rm cut}$ that cleaves the globular cluster sample into two
dependent subgroups based upon whether an individual globular cluster’s
metallicity, ${FeH_{i}}$, is above or below this limit. However, as the
metallicity uncertainties, $\sigma_{FeH_{i}}$, can be significant, this is
analogous to a curve fitting problem with error bars on the $x$-values (Hogg
et al., 2010). Thus, we included all $n=278$ true metallicities as extra
parameters to be inferred. We take the prior for each true metallicity to be a
normal distribution centered at the measured value and with the appropriate
standard deviation:

$\displaystyle FeH_{i}$
$\displaystyle\sim\textnormal{Normal}\left(\widehat{FeH_{i}},\sigma_{FeH_{i}}^{2}\right),$
(3)

where $\widehat{FeH_{i}}$ is the given measurement and $\sigma_{FeH_{i}}$ the
reported uncertainty.

For simplicity we did not use a hierarchical model, so the metallicities do
not ’borrow strength’ from each other (i.e., there is no tendency for an
unknown metallicity value to be dragged towards the preponderance of the other
metallicities).

Figure 1: Corner plot for the best-fit model as outlined in the text. Here,
the blue distributions correspond to the metal-poor component of the globular
cluster population, whilst the red represents the metal-rich component. The
green distributions show the relationship of the parameters of the two
populations against each other. The final row presents the metallicity cut.
Note that the $F_{\rm II}$ kinematic model has two additional scale-length
parameters $L_{k}$, which are constrained to be relatively small and are
omitted from the plot for clarity. The best-fit model parameters are tabulated
in Table 4.

###  2.4 Likelihood Calculation

For a given rotational model $M_{j}=\\{V_{j},\,S_{k},\,{\rm or}\ F_{j}\\}$
(where $j={\rm I\ or\ II}$) which are dependent upon a set of parameters
$\Theta_{k}$, and assuming a normal distribution for the velocity
uncertainties on the globular clusters, we define the corresponding likelihood
to be

${\cal L}_{k}({\cal
D}_{i})=\frac{1}{\sqrt{2\pi}\sigma_{ki}}\exp{\left(-\frac{(M_{j}({\cal
D}_{i}|\Theta_{k})-v_{i})^{2}}{2\sigma^{2}_{ki}}\right)}\,,$ (4)

where the velocity dispersion and the measurement uncertainty are added in
quadrature such that

$\sigma_{ki}^{2}=\sigma_{v_{i}}^{2}+\Sigma_{k}^{2}({\cal
D}_{i}|\Theta_{k})\,.$ (5)

When considering a single-population model, such that $k=0$, the total
likelihood is given by

${\cal L}=\prod_{i=1}^{n}{\cal L}_{0}({\cal D}_{i}|\Theta_{0})\,.$ (6)

However, when considering a model with two populations split by metallicity
such that $k=0,1$, it is appropriate to select the corresponding kinematic
component for the likelihood term in the product.

${\cal L}=\prod_{i=1}^{n}\left\\{\begin{array}[]{lr}{\cal L}_{0}({\cal
D}_{i}|\Theta_{0}),&FeH_{i}<FeH_{\rm cut}\\\ {\cal L}_{1}({\cal
D}_{i}|\Theta_{1}),&FeH_{i}\geq FeH_{\rm cut}\\\ \end{array}\right.$ (7)

As typical in numerical calculations involving likelihood, practical
considerations require that the log of the likelihood is computed.

For each model, we determine the posterior distribution for the parameters
given the data, and the marginal likelihood of the model as a whole, using
Diffusive Nested Sampling (Brewer et al., 2011; Brewer & Foreman-Mackey,
2018), a variant of Nested Sampling (Skilling, 2004, 2006).

###  2.5 Priors

For simplicity, we employed uniform priors for all parameters, with
astronomically informed limits; these are given in Table 2. In Bayesian model
comparison, the priors can have a significant effect on the marginal
likelihoods. As a consequence, the marginal likelihoods we compute in this
study are sensitive to the prior widths. However, this is not a major concern
for two reasons. Firstly, we used the same priors for all parameters shared in
common between models so they have similar predictive ability, and secondly,
the main caveat to Bayesian model comparison with uniform priors is that very
wide limits might overly penalise the more complex model. Since we found that
the two-component model was favoured despite the wide uniform priors, a more
sophisticated choice of priors would likely strengthen this result, rather
than weakening it.

Model | $\ln Z$ | $Z/Z_{\rm max}$
---|---|---
$V_{\rm I}$ | -1780.76 | 0.051
$S_{\rm I}$ | -1785.82 | 0.00032
$F_{\rm I}$ | -1780.27 | 0.083
$V_{\rm II}$ | -1778.39 | 0.54
$S_{\rm II}$ | -1784.79 | 0.00090
$F_{\rm II}$ | -1777.78 | 1
$\top$ | -1779.05 |
Table 3: Marginal likelihoods or evidences for the three single-component
models followed by the three two-component models. The flat kinematics $F$
predicts the data best, followed closely by the $V$ form. The solid-body $S$
kinematic model is heavily disfavoured. Two components are favoured over a
single component model in every case, by Bayes factors of around
$\exp(2)\approx 10$. Finally, we present the marginal likelihood for the top
statement $\top=V_{\rm I}\vee S_{\rm I}\vee F_{\rm I}\vee V_{\rm II}\vee
S_{\rm II}\vee F_{\rm II}$, which is the logical disjunction of all six
models, assuming $1/6$ prior probability each. This is the evidence for the
statement that one of these six models is true.

##  3 Results and Discussion

Figure 2: The best-fit model for the inner globular cluster population of the
Andromeda Galaxy (c.f. Figure 1). In each panel, the globular clusters are
represented as filled circles colour-coded with their metallicity in the left-
hand panel, and with velocity in the right-hand panel. The size of each circle
reflects the metallicity of each globular cluster, with more metal-poor
clusters being larger (scaling is set by a logistic function centred at
$FeH=-2.2$ and a width of 0.2 in metallicity). The transparency of each
globular cluster symbol in each panel has been adjusted based upon
metallicity, such that the left-hand panel emphasises the metal-rich
component, whereas the right-hand panel emphasises the metal-poor population.
The disk of the Andromeda Galaxy is represented as ellipses. In the left-hand
panel, the orientation of the two rotational components are presented, with
the metal-poor axis of rotation in blue and the metal-richer component in red,
with the width of the distributions representing the uncertainty on these
quantities from the Bayesian analysis. In the right-hand panel, the kinematic
model is represented, with the blue arrow representing the angular momentum of
the population. The red arrow is the angular momentum of the outer halo
globular cluster population that is associated with substructure and
represents a recent accretion into Andromeda. Note that, for the purposes of
this figure, the tangent plane coordinates are labelled with $(\xi,\eta)$ to
aid comparison to previous work. However, for the analysis presented in this
paper, $(x,y)$ are used to represent these coordinates.

###  3.1 Bayesian Evidence

The Bayesian evidences from our nested sampling exploration of the posterior
probability spaces are presented in Table 3. One immediate result is that,
irrespective of the rotational model under consideration, the Bayesian
evidence explicitly favours the two component rotational model over the single
component, with a separation of the populations at metallicity $FeH\sim-2.2$
consistent across the models. For the best-fitting model, $\\{F_{\rm II}\\}$,
the ratio of the Bayesian evidence between the single and two components fits,
known as the Bayes factor, is 12.06, with it taken as "strong" evidence that
the two component model is the better description of the data (Kass & Raftery,
1995). We summarise the properties of the best-fit model in the corner plot
presented in Figure 1 and tabulated in Table 4; again, we emphasise that the
characteristics of this best-fit are seen across all of the two component
fits. Essentially, the the metal-poor component, which consists of $\sim
10-20$ globular clusters, has a rotational velocity amplitude of $\sim
340\textrm{km/s}$, substantially higher than the more metal-rich population
(with a rotational amplitude of $\sim 55\textrm{km/s}$), blackalthough we note
that the large rotational amplitude is potentially driven by a small number of
globular clusters at the velocity extremes. Furthermore, the velocity
dispersion of the metal-poor component, at $\sim 200\textrm{km/s}$, is higher
than that for the more metal-rich component $(\sim 170\textrm{km/s})$. Whilst
these values refer to the $50^{\rm th}$ percentile of the posterior
probability distributions, it should be noted that some of these show
substantial skewness, and we note that with this the velocity dispersions are
roughly consistent within the uncertainties. blackWe comment on this
potentially large velocity dispersion in the conclusions to this paper.

Figure 3: blackSummary of the analysis of the entire globular cluster population, including those with $FeH\geq-0.4$ (see Section 2.1). The left-hand panel presents a two component model where the prior range of the metallicity cut was expanded to encompass the entire span of metallicity, showing distinct peaks at $FeH\sim-2$ and $FeH\sim-0.4$. The right-hand panel presents the distribution of rotational axes from a three component model with hard metallicity cuts at $FeH=-2.0$ and $FeH=-0.4$, with $FeH<-2.0$ in blue, $-2.0\leq FeH\leq-0.4$ in red and $FeH>-0.4$ in green. Clearly, the rotational axes of the sample with $FeH>-2.0$ are aligned with the rotational axis of the Andromeda Galaxy, whereas the more metal-poor sample possesses a distinctly offset rotational axis; it is this we identify as the Dulais Structure.  Parameter | Value | Units
---|---|---
$A_{0}$ | $340.66$ | ${}^{+269.83}_{-128.20}$ | km/s
$\phi_{0}$ | $1.41$ | ${}^{+0.36}_{-0.70}$ | radians
$\Sigma_{0}$ | $204.34$ | ${}^{+112.00}_{-81.89}$ | km/s
$s_{0}$ | $-0.80$ | ${}^{+1.32}_{-1.84}$ | (degrees)-1
$L_{0}$ | $0.77$ | ${}^{+0.79}_{-0.62}$ | degrees
$A_{1}$ | $57.56$ | ${}^{+10.40}_{-9.77}$ | km/s
$\phi_{1}$ | $2.72$ | ${}^{+0.13}_{-0.11}$ | radians
$\Sigma_{1}$ | $168.93$ | ${}^{+12.57}_{-11.50}$ | km/s
$s_{1}$ | $-0.53$ | ${}^{+0.13}_{-0.14}$ | (degrees)-1
$L_{1}$ | $0.13$ | ${}^{+0.09}_{-0.07}$ | degrees
$FeH_{cut}$ | $-2.29$ | ${}^{+0.20}_{-0.22}$ |
Table 4:  The best-fit parameter values for the favoured $F_{\rm II}$ model,
shown graphically in the corner plot in Figure 1. The values reflect the
$50^{th}$ quantile of the posterior distribution, with the uncertainty
determined from the $16^{th}$ and $84^{th}$ quantiles.

blackBefore closing this section, we return to the fact that metallicity cut
was made to the initial sample at $FeH\geq-0.4$ to remove the galactic
component identified by Caldwell & Romanowsky (2016, see Section 2.1). To
ensure this selection has not biased the findings in this paper, several
additional tests were undertaken which are summarised here. Essentially,
single and multiple rotational components, using model $F$ (see Table 1), were
fit for each of the components, and the Bayesian evidence calculated. In
fitting a single rotational component, the resultant Bayesian evidence
$Z_{1}=−2165.62$. The second model comprised of two components, but the prior
on the metallicity cut was expanded to cover the entire metallicity range. The
resultant evidence for this was $Z_{2}=−2162.58$, demonstrating that this is
favoured over the single component model. However, the posterior for the
metallicity cut, presented in the left-hand panel of Figure 3, possesses two
distinct peaks, one at $FeH\sim-2$ and the other at $FeH\sim-0.4$, indicating
that the modelling suggests that the globular cluster populations could be
split into three. Fixing the metallicity cut at $FeH=-2.0$ recovers the
results presented in this paper, namely that there are two kinematically
distinct populations with significantly different rotational axes. However,
fixing the metallicity cut at $FeH=-0.4$ shows that both populations share a
common axis of rotation, with those at $FeH>-0.4$ possessing a rotation
amplitude of $\sim 55$ km/s, whilst the lower metallicity population possesses
a rotational amplitude of $\sim 160$ km/s. Finally, a three component model
was fit to the entire population, with fixed metallicity cuts at $FeH=-2.0$
and $FeH=-0.4$. The resulting Bayesian evidence was $Z_{3}=−2156.5$,
demonstrating that this is the favoured model. In the right-hand panel, we
present the rotational axes of the three components from the posterior
exploration, with $FeH<-2.0$ in blue, $-2.0\leq FeH\leq-0.4$ in red and
$FeH>-0.4$ in green. The distribution of red and green axes overlap and are
aligned with rotational axis of Andromeda, whilst the most metal-poor
component in blue is distinctly offset from these. We conclude that
disregarding globular clusters with $FeH>-0.4$ has not biased the findings in
this paper.

###  3.2 Exploring the Two Component Model

If we assume that the inner globular clusters of Andromeda represent a mix of
those formed in-situ coupled with phase-mixed accretion events, a natural
explanation for the properties of the metal-poor globular clusters is that it
is simply a kinematically hotter component of the overall population. However,
if this was the correct interpretation, we would expect the orientation of the
rotation of the globular cluster population to reflect that of the overall
Andromeda Galaxy. To explore this, we present the orbital orientation of the
two components for the best-fit model is displayed in Figure 2. The globular
clusters are shown shown as filled circles colour-coded with their metallicity
(left-hand panel) and velocity (right-hand panel). In both panels, more metal-
poor clusters being larger. In the left-hand panel, the more metal-rich
globular clusters are emphasised through their opacity, whereas it is the more
metal-poor clusters emphasised in the right-hand panel. Underlying the
globular clusters are solid lines denoting the major and minor axes of the
Andromeda Galaxy, emphasised with overlying ellipses.

In the left-hand panel of Figure 2 the axes of rotation of the metal-poor
component (blue) and metal-rich component (red) are presented underlying the
globular cluster population; again these are taken from the $50^{\rm th}$
percentile of the posterior probability distribution, with one $\sigma$
uncertainties corresponding to the $16^{\rm th}$ and $84^{\rm th}$ percentiles
(see Figure 1). The rotational axis of the metal-rich component is reasonably
aligned with the minor axis of the Andromeda Galaxy, and so we conclude that
their kinematic properties are related to the larger-scale galactic rotation.
However, the rotational axis of the metal-poor component is substantially
offset from the minor axis of the Andromeda Galaxy by $\sim 75\degr$. This is
emphasised in the angular momentum vectors of the two components, with the
metal-rich population at $PA\sim 115\degr$ whereas the metal-poor component is
at $PA\sim 190\degr$333 The best-fit values for the orientation of the
rotation are $\phi_{0}=1.41$ and $\phi_{1}=2.72$ radians for the Dulais
Structure and dominant galactic population respectively (Table 4). However,
given the form of the underlying model (Table 1), the direction of the
rotational vector (given the right-hand rule) is these values $\pm\pi$. This
corresponds to $\theta_{0}\sim-100\degr$ and $\theta_{1}\sim-25\degr$ in the
model coordinates (Equation 1). Given that $\theta$ and $PA$ are related
through $\theta+PA=90\degr$, this gives the quoted $PA$s for the orientations
of the spin vectors of the two components.  Hence, we concluded that the
metal-poor globular clusters reveal the presence of a kinematically distinct
population from the overall inner globular clusters; it is this that we
identify as the Dulais Structure. We do, however, add a word of caution in
interpreting these findings on the metallicity of the Dulais Structure as it
is unlikely that globular cluster populations can cleanly be separated by
metallicity alone and some of the low metallicity clusters are likely to be
drawn from the galactic population, whereas the Dulais Structure will have
globular clusters above $FeH_{\rm cut}$. What the results of this study show
is that the kinematics of the Dulais Structure dominate at low metallicity,
whereas it is lost against the kinematics of the galactic population at higher
metallicities. This does indicate that the Dulais Structure is, overall, more
metal-poor than the galactic population of globular clusters. We explore this
in more detail below.

Figure 4: The ratio of the likelihoods, $R$, drawn from the two kinematic
components of the best-fit model for the globular clusters in our sample,
presented as a function of metallicity. The smaller the value of $R$, the more
likely the globular cluster is drawn from the galactic populations. Those with
$R<2$ are coloured blue, whereas those with $R>2$ are red.

The right-hand panel of Figure 2 presents the best-fit rotational model for
the Dulais Structure component, colour-coded for velocity; clearly the
direction of rotation is distinctly different to the rotation of Andromeda,
which is aligned along the minor axis. Also in this figure, the blue arrow
denotes the direction of the angular momentum vector for this rotation of the
Dulais Structure component. The red arrow, however, presents the direction of
the angular momentum for the outer globular cluster population that was
associated with underlying stellar substructure, and therefore indicative of a
relatively recent accretion event into the halo of the Andromeda Galaxy
(Mackey et al., 2019b). Clearly, given that the direction and rotational sense
of these inner and outer globular cluster populations are extremely similar,
it is natural to suggest that the two components are related. One obvious
conclusion is that the Dulais Structure represents the remnants of a stellar
system accreted as part of the larger-scale accretion seen in the outer halo
of Andromeda.

Figure 5: The total angular momentum, $L_{T}$, (top panel) and x-component of
the angular momentum, $L_{x}$ (bottom panel) of the globular cluster
population, colour-coded as described in Figure 4. The upper panel reveals
that those globular clusters more likely to be members of the Dulais Structure
(coloured red) are at larger absolute angular momentum than the overall
population. The lower panel illustrates that those globular clusters more
likely members of the Dulais Structure are clustered at positive values and
are distinct from the overall galactic population. We note that the
significant clustering in $L_{x}$ is apparent as the rotational axis of the
Dulais Structure is aligned with the y-axis of the tangent plane coordinates.

###  3.3 Posterior Analysis

As discussed previously, given the limitations of the data the sub-populations
of globular clusters were defined solely on a metallicity cut. However, care
must be taken in interpreting the metallicity distribution of the Dulais
Structure based upon the metallicity cut as the dominant population is likely
to extend below the cut, and conversely, there are likely to be globular
clusters in the Dulais Structure which are more metal-rich than the
metallicity cut. Hence, the metallicity cut is more likely to be indicative of
where the Dulais Structure population of globular clusters comes to dominate
over the larger galactic population. Given the best-fit model presented in
Section 3.2, we here undertake an analysis of the posterior distributions to
provide further insights into metallicity distribution of the Dulais
Structure, although we readily note that this is indicative only, and more
detailed modelling with more complete data is required to make robust
statements.

For each globular cluster in our sample, we assesses which population they are
more likely drawn from by considering the ratio of the likelihoods, $R$, of
the two models. We present this in Figure 4, noting that small $R$ values
correspond to globular clusters being more likely members of the more dominant
galactic population rather than the Dulais Structure. For the sake of clarity
in the coming analysis, we take an (arbitrary) line at $R=2$ and colour-code
this sub-population of globular clusters in red. As well as a low metallicity
tail, we see this sub-population has members that extend to higher
metallicity.

With this separation of the populations, it is interesting to explore whether
those more likely to be members of Dulais Structure possess clustering across
any other properties. Of course, without full 3D positions and velocities, we
are working within a limited representation of phase space. However, we define
two quantities representing the projected angular momentum of the form;

$L_{T}=rv\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ L_{x}=xv$ (8)

where ${x,r}$ are the tangent plane coordinates and $v$ is line-of-sight
velocity as the described in Section 2.1. These two quantities are presented
in Figure 5 as a function of the metallicity, incorporating the same colour-
coding with respect to $R$ as shown in Figure 4. In the upper panel, it is
immediately apparent that the globular clusters are more likely members of the
Dulais Structure are found at higher total angular momentum than the dominant
galactic population. This is reinforced in the lower panel of Figure 5 which
presents $L_{x}$ as a function of metallicity, which delineates the potential
members of the Dulais Structure from the dominant population by having large
positive values of this quantity. We note that this clustering is significant
in this quantity due to the alignment of the rotational axis of the Dulais
Structure with the $y$-coordinate of the tangent plane. It is again worth
emphasising that whilst this posterior analysis cannot be considered a robust
method of separating the dominant and Dulais Structure populations, such
clustering in terms of angular momentum adds credence to there being two
distinct sub-populations of globular clusters present.

Finally, in Figure 6, we present the spatial distribution of the globular
cluster population (c.f. Figure 2) colour-coded with the relative likelihood
parameter, $R$, with blue representing those most likely drawn from the
dominant population, whilst those from the Dulais Structure are in red. Here,
there is no obvious division between the two populations, and more in-depth
chemical and kinematic analysis will be needed to cleanly separate the
dominant population from the Dulais Structure.

##  4 Conclusions

In this paper we have presented a Bayesian evidence based analysis of the
inner globular cluster population of the Andromeda galaxy, finding that two
kinematic component models represent a better fit to the observed population
properties than single component models. The two sub-populations are defined
in relation to a metallicity cut, with the more metal-rich component being
consistent with galactic rotation. However, the orbital axis of the more
metal-poor component is highly offset with regards to this rotation, but is
consistent with the kinematics of the outer halo globular cluster population.
This inner substructure, which we name Dulais Structure, is therefore
potentially part of this large-scale structure projected onto the inner part
of the Andromeda Galaxy. Mackey et al. (2019b) conclude that the outer halo
globular clusters which are still associated with stellar substructure were
accreted in the last few billion years, and hence the Dulais Structure
represents another feature of that accretion. If the 10-20 globular clusters
of the Dulais Structure do represent a single accreted object, then this would
correspond to a progenitor halo mass of $\sim 10^{11}M_{\odot}$ (c.f. Harris
et al., 2015), although we caution that the full dynamical state and
membership of this structure is essentially unknown, blackalthough the
potentially large velocity dispersion found for Dulais Structure are
suggestive of the accretion of a galaxy group into the halo of Andromeda.

Figure 6: The spatial distribution of the globular clusters under
consideration, with coordinates presented as $(\xi,\eta)$ for comparison with
Figure 2. Again, the globular clusters have been coloured with respect to the
relative likelihood ratio, $R$, of belonging to the galactic populations
(blue) or Dulais Structure (red), with the division between the two sub-
populations presented in Figure 4.

Before closing, we note that, in terms of accretion, the inner regions of
Andromeda are dominated by the presence of the Giant Stellar Stream (GSS)
which snakes its way around the stellar disk of the galaxy and is responsible
for the stellar shells which are apparent (Ibata et al., 2007). Numerical
simulations of the tidal formation of the GSS place the present day location
of the remnants of the progenitor in front of the stellar disk of Andromeda
(Fardal et al., 2013; Kirihara et al., 2017), similar to where we find the
Dulais Structure; blackintriguingly, an estimate of the projected total
angular momentum, $L_{T}$ (see Figure 5) for the GSS is of order $\pm 150$ deg
km/s (see Figure 8 of Fardal et al., 2013), similar to that seen for the
Dulais Structure although a more detail dynamical exploration is needed to
identify whether this is more than coincidence. If the Dulais Structure is the
remains of the progenitor of the GSS, then the angular alignment with the
outer halo globular cluster population presents the intriguing possibility
that the inner and outer stellar debris represents a single accretion event of
a group of galaxies in the relatively recent history of the Andromeda Galaxy.
This will be explored in a future contribution.

## Acknowledgements

blackWe thank Jorge Peñarrubia and the anonymous referee for useful comments
and insights that improved the quality of the paper. This project was
originally conceived in discussions between GFL, DM and AF, and the Bayesian
exploration was planned in consultation with BJB. The initial kinematic
investigations were undertaken as part of two separate honours projects by TA
and YCL under the supervision of GFL and BJB respectively. GFL has received no
funding to support this research.

## Data Availability

The data and code used in this study will be made available with a reasonable
request to the authors. The globular cluster data that formed the basis of the
analysis was obtained from MMT spectroscopy, published by Caldwell &
Romanowsky (2016), and are available at
https://www.cfa.harvard.edu/oir/eg/m31clusters/M31_Hectospec.html.

## References

  * Barmby et al. (2000) Barmby P., Huchra J. P., Brodie J. P., Forbes D. A., Schroder L. L., Grillmair C. J., 2000, AJ, 119, 727
  * Belokurov et al. (2018) Belokurov V., Erkal D., Evans N. W., Koposov S. E., Deason A. J., 2018, MNRAS, 478, 611
  * Brewer & Foreman-Mackey (2018) Brewer B., Foreman-Mackey D., 2018, Journal of Statistical Software, Articles, 86, 1
  * Brewer et al. (2011) Brewer B. J., Pártay L. B., Csányi G., 2011, Statistics and Computing, 21, 649
  * Burstein et al. (1984) Burstein D., Faber S. M., Gaskell C. M., Krumm N., 1984, ApJ, 287, 586
  * Caldwell & Romanowsky (2016) Caldwell N., Romanowsky A. J., 2016, ApJ, 824, 42
  * Crampton et al. (1985) Crampton D., Cowley A. P., Schade D., Chayer P., 1985, ApJ, 288, 494
  * Elson & Walterbos (1988) Elson R. A., Walterbos R. A. M., 1988, ApJ, 333, 594
  * Fabricant et al. (2005) Fabricant D., et al., 2005, PASP, 117, 1411
  * Fardal et al. (2013) Fardal M. A., et al., 2013, MNRAS, 434, 2779
  * Freeman (1983) Freeman K. C., 1983, in Athanassoula E., ed., 0 Vol. 100, Internal Kinematics and Dynamics of Galaxies. pp 359–364
  * Harris et al. (2015) Harris W. E., Harris G. L., Hudson M. J., 2015, ApJ, 806, 36
  * Hartwick & Sargent (1974) Hartwick F. D. A., Sargent W. L. W., 1974, ApJ, 190, 283
  * Helmi et al. (2018) Helmi A., Babusiaux C., Koppelman H. H., Massari D., Veljanoski J., Brown A. G. A., 2018, Nature, 563, 85
  * Hogg et al. (2010) Hogg D. W., Bovy J., Lang D., 2010, arXiv e-prints, p. arXiv:1008.4686
  * Huchra et al. (1982) Huchra J., Stauffer J., van Speybroeck L., 1982, ApJ, 259, L57
  * Ibata et al. (1994) Ibata R. A., Gilmore G., Irwin M. J., 1994, Nature, 370, 194
  * Ibata et al. (2007) Ibata R., Martin N. F., Irwin M., Chapman S., Ferguson A. M. N., Lewis G. F., McConnachie A. W., 2007, ApJ, 671, 1591
  * Ibata et al. (2020) Ibata R., Bellazzini M., Thomas G., Malhan K., Martin N., Famaey B., Siebert A., 2020, ApJ, 891, L19
  * Kass & Raftery (1995) Kass R. E., Raftery A. E., 1995, Journal of the American Statistical Association, 90, 773
  * Kirihara et al. (2017) Kirihara T., Miki Y., Mori M., Kawaguchi T., Rich R. M., 2017, MNRAS, 464, 3509
  * Kruijssen et al. (2020) Kruijssen J. M. D., et al., 2020, MNRAS, 498, 2472
  * Li et al. (2022) Li T. S., et al., 2022, ApJ, 928, 30
  * Mackey et al. (2019a) Mackey A. D., et al., 2019a, MNRAS, 484, 1756
  * Mackey et al. (2019b) Mackey D., et al., 2019b, Nature, 574, 69
  * Mateu (2022) Mateu C., 2022, arXiv e-prints, p. arXiv:2204.10326
  * McConnachie et al. (2018) McConnachie A. W., et al., 2018, ApJ, 868, 55
  * Myeong et al. (2018) Myeong G. C., Evans N. W., Belokurov V., Sanders J. L., Koposov S. E., 2018, ApJ, 863, L28
  * Perrett et al. (2002) Perrett K. M., Bridges T. J., Hanes D. A., Irwin M. J., Brodie J. P., Carter D., Huchra J. P., Watson F. G., 2002, AJ, 123, 2490
  * Perrett et al. (2003) Perrett K. M., Stiff D. A., Hanes D. A., Bridges T. J., 2003, ApJ, 589, 790
  * Skilling (2004) Skilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds, American Institute of Physics Conference Series Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering. pp 395–405, doi:10.1063/1.1835238
  * Skilling (2006) Skilling J., 2006, Bayesian Analysis, 1, 833
  * Veljanoski et al. (2013) Veljanoski J., et al., 2013, ApJ, 768, L33
  * Veljanoski et al. (2014) Veljanoski J., et al., 2014, MNRAS, 442, 2929
  * Wechsler et al. (2002) Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, ApJ, 568, 52
  * van den Bergh (1969) van den Bergh S., 1969, ApJS, 19, 145

##  Appendix A Catalogue of Globular Clusters

Table 5: The 22 globular clusters with $R>2$, identified as potential members of the Dulais Structure (Section 3.3). For completeness, the metallicities and velocities for the entire sample of globular clusters are presented in Figure 7 with the associated error bars, demonstrating that there are no systematic issues in these quantities for those identified as being part of the Dulais Structure. Note that the tangent plane coordinates, $(\xi,\eta)$, correspond to the $(x,y)$ coordinates used in the paper.  Name | RA | Dec | $\xi\ (\degr)$ | $\eta\ (\degr)$ |  $v$ (km/s) | $\sigma_{v}$ | $\widehat{FeH_{i}}$ | $\sigma_{FeH}$ | $R$
---|---|---|---|---|---|---|---|---|---
B298-G021 | 9.50092 | 40.73217 | -0.8970 | -0.5311 | -244.1 | 17.6 | -1.8 | 0.3 | 10.8
B311-G033 | 9.89050 | 40.52075 | -0.6037 | -0.7459 | -222.1 | 20.3 | -1.9 | 0.2 | 5.3
B336-G067 | 10.19833 | 42.14533 | -0.3606 | 0.8772 | -306.0 | 27.4 | -2.5 | 0.6 | 72.0
B017-G070 | 10.20304 | 41.20197 | -0.3623 | -0.0663 | -233.9 | 3.0 | -0.8 | 0.2 | 2.8
B248 | 10.28308 | 40.88361 | -0.3036 | -0.3850 | -275.3 | 3.1 | -1.4 | 0.2 | 3.0
B020D-G089 | 10.32179 | 41.13586 | -0.2732 | -0.1328 | -257.1 | 3.1 | -2.1 | 0.3 | 2.5
B032-G093 | 10.33962 | 41.29172 | -0.2592 | 0.0230 | -224.7 | 7.2 | -0.7 | 0.2 | 2.2
B034-G096 | 10.36717 | 40.89711 | -0.2399 | -0.3717 | -255.5 | 10.2 | -0.6 | 0.2 | 2.1
B036 | 10.38679 | 41.43475 | -0.2233 | 0.1659 | -208.6 | 3.0 | -0.8 | 0.2 | 2.1
B198-G249 | 10.95879 | 41.53128 | 0.2053 | 0.2623 | 298.4 | 8.3 | -1.0 | 0.2 | 2.4
B217-G269 | 11.04417 | 41.39756 | 0.2697 | 0.1288 | 284.6 | 3.0 | -0.8 | 0.2 | 3.0
B229-G282 | 11.14096 | 41.64125 | 0.3411 | 0.3729 | 280.8 | 4.5 | -2.1 | 0.3 | 3.6
B233-G287 | 11.17550 | 41.73183 | 0.3664 | 0.4636 | 228.6 | 11.6 | -1.1 | 0.2 | 2.7
B283-G296 | 11.23071 | 41.28339 | 0.4104 | 0.0154 | 211.0 | 3.0 | -0.8 | 0.2 | 2.7
B235-G297 | 11.24137 | 41.49000 | 0.4171 | 0.2221 | 207.8 | 3.0 | -0.9 | 0.2 | 2.6
B237-G299 | 11.28842 | 41.37628 | 0.4531 | 0.1086 | 199.2 | 3.1 | -1.8 | 0.2 | 2.5
B238-G301 | 11.31113 | 41.32697 | 0.4705 | 0.0594 | 257.9 | 3.0 | -0.7 | 0.2 | 4.8
B240-G302 | 11.35433 | 41.10614 | 0.5047 | -0.1612 | 243.6 | 3.6 | -1.5 | 0.2 | 5.4
B270D | 11.45504 | 41.03033 | 0.5812 | -0.2364 | 338.9 | 17.8 | -2.1 | 0.4 | 28.8
B381-G315 | 11.52725 | 41.34967 | 0.6326 | 0.0835 | 217.5 | 4.4 | -1.1 | 0.2 | 4.1
B387-G323 | 11.63963 | 40.73706 | 0.7237 | -0.5283 | 310.4 | 14.3 | -2.2 | 0.2 | 83.4
B397-G336 | 11.86346 | 41.20289 | 0.8870 | -0.0604 | 179.8 | 14.7 | -1.2 | 0.2 | 2.3
Figure 7: The globular cluster metallicities and velocities, with associated
error bars, used in this study. Those with $R>2$, and hence identified as
being most likely a member of the Dulais Structure in a posterior analysis
(Section 3.3). It should be noted that there are no significant systematics
with regards to the error bars for those globular clusters identified as part
of the Dulais Structure.


# Deeply Shared Filter Bases for Parameter-Efficient Convolutional Neural
Networks

Woochul Kang
Department of Embedded Systems Engineering
Incheon National University
Yeonsu-gu, South Korea 22012
<EMAIL_ADDRESS>
Daeyeon Kim
Dain Technology, Inc.
<EMAIL_ADDRESS>
W. Kang and D. Kim are contributed equally. This work was conducted at Incheon
National University.

###### Abstract

Modern convolutional neural networks (CNNs) have massive identical convolution
blocks, and, hence, recursive sharing of parameters across these blocks has
been proposed to reduce the amount of parameters. However, naive sharing of
parameters poses many challenges such as limited representational power and
the vanishing/exploding gradients problem of recursively shared parameters. In
this paper, we present a recursive convolution block design and training
method, in which a recursively shareable part, or a filter basis, is separated
and learned while effectively avoiding the vanishing/exploding gradients
problem during training. We show that the unwieldy vanishing/exploding
gradients problem can be controlled by enforcing the elements of the filter
basis orthonormal, and empirically demonstrate that the proposed
_orthogonality regularization_ improves the flow of gradients during training.
Experimental results on image classification and object detection show that
our approach, unlike previous parameter-sharing approaches, does not trade
performance to save parameters and consistently outperforms overparameterized
counterpart networks. This superior performance demonstrates that the proposed
recursive convolution block design and the orthogonality regularization not
only prevent performance degradation, but also consistently improve the
representation capability while a significant amount of parameters are
recursively shared.

##  1 Introduction

Modern convolutional neural networks (CNNs) such as ResNets have massive
identical convolution blocks and recent analytic studies (Jastrzebski et al.,
2018) show that these blocks perform mostly iterative refinement of features
rather than learning new features. Inspired by these results, recursive
sharing of weights has been studied as a promising direction to parameter-
efficient CNNs (Jastrzebski et al., 2018; Guo et al., 2019; Savarese & Maire,
2019). However, naive sharing of parameters across many convolution layers
incurs several problems. First of all, recursive sharing of parameters can
result in the vanishing and the exploding gradients problem, which is one of
the main reasons that recurrent neural networks (RNNs) are so hard to train
properly (Pascanu et al., 2013; Jastrzebski et al., 2018). Another problem is
that overall representation power can be limited by iterative sharing of
parameters. Due to these challenges, most compression approaches based on
parameter-sharing suffer from performance degradation.

In this work, we conjecture that convolution layers or blocks can be separated
into inherently shareable parts and non-shareable parts, and can be trained
effectively by avoiding the vanishing and the exploding gradients problem. To
achieve this, for a full convolution operator, we first replace it with a
factorized version that splits the convolution operator into two separate
operators; one operator with inherently shareable filters, called a _filter
basis_ , and the other operator with non-shared filters, called
_coefficients_. When successive convolution blocks share a common filter
basis, they are positioned in the same vector subspace. However, their
representation capability is retained through non-shared coefficients that
learn diverse features by linearly combining the shared filter basis.

By separating shareable parts from non-shareable parts, we can impose
desirable properties on the shared parameters. To avoid performance
degradation from recursive sharing of parameters, we propose the
_orthogonality regularization_ , in which the vanishing/exploding gradients
problem is controlled by enforcing the elements of a shared filter basis
orthonormal during training. We both theoretically and empirically show that
the proposed orthogonality regularization improves the flow of the gradients
during training and reduces the redundancy in parameters effectively.

For efficient CNNs such as MobileNets (Howard et al., 2017), we do not need to
factorize convolution operators to uncover a shared filter basis since these
networks already have factorized convolution block structures for
computational efficiency. For such networks, our approach can be applied
simply by identifying one or two convolution operators of repeating
convolution blocks as a filter basis that shares weights across the repeating
blocks. Other convolution operators in each convolution block become block-
specific non-shared coefficients.

Since our focus is not on pushing the state-of-the-art performance, we
demonstrate the effectiveness of our work using widely-used models as base
models on image classification and object detection. Without bells and
whistles, simply applying the proposed convolution block design and the
orthogonality regularization saves a significant amount of parameters while
consistently outperforming over-parameterized counterpart networks. For
example, our method can save up to 46.0% of parameters of ResNets while
consistently achieving lower test errors. Even in compact models, such as
MobileNetV2, our approach can achieve further 8-21% parameter savings while
outperforming the original models. This superior performance demonstrates that
the proposed recursive convolution blocks and orthogonality regularization
enables effective learning of better feature representations while a
significant amount of parameters are shared recursively.

##  2 Related Work

Recursive networks and parameter sharing: Recurrent neural networks (RNNs)
(Graves et al., 2013) have been well-studied for temporal and sequential data.
As a generalization of RNNs, recursive variants of CNNs are used extensively
for visual tasks (Socher et al., 2011; Liang & Hu, 2015; Xingjian et al.,
2015; Kim et al., 2016; Zamir et al., 2017). For instance, Eigen et al. (2014)
explore recursive convolutional architectures that share filters across
multiple convolution layers. They show that recurrence with deeper layers
tends to increase performance. However, their recursive architecture shows
worse performance than independent convolution layers due to overfitting. In
contrast, we share only fundamentally shareable parts, or filter bases, that
are separated from conventional convolution operators. By separating shareable
parts from non-shareable parts, we can impose desirable properties on the
shared parameters that prevents the vanishing/exploding gradients problems
without damaging representational capability of the networks.

More recently, Jastrzebski et al. (2018) show that iterative refinement of
features in ResNets suggests that deep networks can potentially leverage
intensive parameter sharing. Guo et al. (2019) introduce a gate unit to
determine whether to jump out of the recursive loop of convolution blocks to
save computational resources. These works show that training recursive
networks with naively shared blocks leads to bad performance due to the
problem of gradient explosion and vanish like RNN (Pascanu et al., 2013;
Vorontsov et al., 2017). In order to deal with the problem of gradient
explosion and vanish, they suggest unshared batch normalization strategy. In
our work, we propose the orthogonality regularization of shared filter bases
to further address this problem.

A few recent works generate convolution filters by combining shared building
blocks (Bhalgat et al., 2020; Savarese & Maire, 2019; Yang et al., 2019; Qiu
et al., 2018). Although these approaches are different in details, they are
similar to our work since they all share a set of filters and combine them to
build layer-specific filters. However, they save parameters at the cost of
accuracy loss. For example, in ResNet50 on ImageNet, Yang et al. (2019)’s work
incurs 0.9% accuracy drop, while our method achieves performance improvement.
We believe that this performance gap comes from the gradients issue of shared
parameters and these previous works could benefit more from the improved flow
of gradients of our approach.

Model compression and efficient convolution block design: Reducing storage and
inference time of CNNs has been an important research topic for both resource
constrained mobile/embedded systems and energy-hungry data centers. A number
of research techniques have been developed such as filter pruning (LeCun et
al., 1990; Polyak & Wolf, 2015; Li et al., 2017; He et al., 2017b), low-rank
factorization (Denton et al., 2014; Jaderberg et al., 2014), and quantization
(Han et al., 2016), to name a few. Several model compression techniques
factorize trained filters to reduce the computation complexity (Zhang et al.,
2015; Li et al., 2019a). Our proposed block structure also benefit from the
reduced computational complexity too by factorizing convolution filters.
However, most previous compression techniques have been suggested as post-
processing steps that are applied after initial training. Therefore, their
accuracy is usually bounded by the approximated base models. In contrast, our
primary goal in factorizing convolution operators is to find a common filter
basis shared by iterative convolution layers, so the weights of the shared
filter basis are learned from scratch, not from decomposing pretrained
filters. Hence, unlike previous compression techniques, the performance of our
approach is not limited by the original models, and the experimental results
show that our parameter-sharing approach consistently outperforms the
overparameterized counterpart networks on various datasets and tasks.

##  3 Proposed Method

Figure 1: A full convolution operator in (a) with a filter $W$ can be replaced
by two separate convolution operators, as in (b): one operator with a _filter
basis_ ($W_{basis}$) and the other operator with _coefficients_ ($\alpha$s). A
filter basis is inherently shareable for convolution operators in the same
vector subspace, and, hence, a recursive architecture is suggested, as in (c).
The filter basis shared by recursive convolution layers can be learned through
typical gradient-based training while avoiding potential vanishing/exploding
gradients problem by enforcing its elements orthonormal.

In this section, we describe the design of recursive convolution blocks and
discuss how to train them to deal with the vanishing/exploding gradients
problem.

###  3.1 Recursive Sharing of a Filter Basis

Naive sharing of convolution filters across repetitive layers degrades the
overall performance due to limited representation power of individual layers.
We argue that convolution filters can be separated into fundamentally
shareable components and non-shareable components, and the performance
degradation can be prevented by sharing only fundamentally shareable parts.

More formally, we consider a convolution layer, shown in Figure 1 (a), that
has $S$ input channels, $T$ output channels, and a set of filters $W=\\{W_{t}\
\in R^{k\times k\times S},t\in[1..T]\\}$. Each filter $W_{t}$ can be
decomposed into a filter basis $W_{basis}$ and coefficients $\mathbf{\alpha}$:

$\centering W_{t}=\sum_{r=1}^{R}\alpha^{r}_{t}W^{r}_{basis},\@add@centering$
(1)

where $W_{basis}=\\{W^{r}_{basis}\in\mathbb{R}^{k\times k\times
S},r\in[1..R]\\}$ is a filter basis, and
$\mathbf{\alpha}=\\{\alpha_{t}^{r}\in\mathbb{R},r\in[1..R],t\in[1..T]\\}$ is
scalar coefficients. In Equation 1, $R$ is the rank of the basis. In a typical
convolution layer, output feature maps $V_{t}\in\mathbb{R}^{w\times h\times
T},t\in[1..T]$ are obtained by the convolution between input feature maps
$U\in\mathbb{R}^{w\times h\times S}$ and filters $W_{t},t\in[1..T]$. With
Equation 1, this convolution can be rewritten as follows:

$V_{t}=U\ast
W_{t}=U*\sum_{r=1}^{R}\alpha^{r}_{t}W^{r}_{basis}=\sum_{r=1}^{R}\alpha^{r}_{t}(U*W^{r}_{basis}),\textrm{
where}\ t\in[1..T].$ (2)

In Equation 2, the order of the convolution operation and the linear
combination of filter basis is reordered according to the linearity of
convolution operators. Therefore, a conventional convolution layer can be
replaced with two successive convolution layers as shown in Figure 1-(b).

This factorized convolution filters suggest a recursive architecture, shown in
Figure 1-(c), in which a common filter basis is used for iterative convolution
layers. Since all filters in the recursive loop have a common filter basis,
they are positioned in the same vector subspace. However, their representation
capability is still retained since their coordinates in the subspace are
diversified by their respective non-shared coefficients.

###  3.2 Orthonormality of Shared Filter Bases

For recursive sharing of a filter basis, as in Figure 1-(c), we need to find
an optimal filter basis that can be used by iterative convolution layers
without loss of performance. Although this optimal filter basis can be
searched by typical gradient-based optimization such as SGD, one major problem
is that the exploding/vanishing gradients problem of recursively shared filter
bases can prevent proper search of optimization space (Pascanu et al., 2013).

More formally, we consider a series of $N$ factorized convolution blocks, in
which a filter basis $W_{basis}$ is shared $N$ times, as in Figure 1-(c). Let
$\mathbf{x}^{i}$ be the input of the $i$-th convolution block, and $a^{i+1}$
be the output of the convolution of $\mathbf{x}^{i}$ with the filter basis
$W_{basis}$

$\centering
a^{i}(\mathbf{x}^{i-1})=W_{basis}^{\top}\mathbf{x}^{i-1}.\@add@centering$ (3)

In Equation 3, $W_{basis}\in\mathbb{R}^{k^{2}S\times R}$ is a reshaped filter
basis that has basis elements at its columns. We assume that input x is
properly adapted (e.g., with im2col operations) to express convolutions using
a matrix-matrix multiplication. Since $W_{basis}$ is shared by $N$ recursive
convolution blocks, the gradient of $W_{basis}$ for some loss function $L$ is:

$\centering\frac{\partial L}{\partial W_{basis}}=\sum_{i=1}^{N}\frac{\partial
L}{\partial a^{N}}\prod_{j=i}^{N-1}\left(\frac{\partial a^{j+1}}{\partial
a^{j}}\right)\frac{\partial a^{i}}{\partial W_{basis}},\@add@centering$ (4)

, where

$\centering\frac{\partial a^{j+1}}{\partial a^{j}}=\frac{\partial
a^{j+1}}{\partial\mathbf{x}^{j}}\frac{\partial\mathbf{x}^{j}}{\partial
a^{j}}=W_{basis}\frac{\partial\mathbf{x}^{j}}{\partial a^{j}}\@add@centering$
(5)

If we plug $W_{basis}\frac{\partial\mathbf{x}^{j}}{\partial a^{j}}$ in
Equation 5 into Equation 4, we can see that $\prod\frac{\partial
a^{j+1}}{\partial a^{j}}$ is the term that makes gradients unstable since
$W_{basis}$ is multiplied many times. This exploding/vanishing gradients can
be controlled to a large extent by keeping $W_{basis}$ close to orthogonal
(Vorontsov et al., 2017). For instance, if $W_{basis}$ admits
eigendecomposition, $[W_{basis}]^{N}$ can be rewritten as follows:

$\centering[W_{basis}]^{N}=[Q\Lambda
Q^{-1}]^{N}=Q\Lambda^{N}Q^{-1},\@add@centering$ (6)

where $\Lambda$ is a diagonal matrix with the eigenvalues placed on the
diagonal and $Q$ is a matrix composed of the corresponding eigenvectors. If
$W_{basis}$ is orthogonal, $[W_{basis}]^{N}$ neither explodes nor vanishes,
since all the eigenvalues of an orthogonal matrix have absolute value of 1.
This result shows that a shared filter basis with orthonormal elements ensures
that forward and backward signals neither explode nor vanish.

Based on this result, we propose the _orthogonality regularization_
111Technically, the columns (or rows) of an orthogonal matrix form an
_orthonormal basis_. However, we use the term _orthogonality regularization_
since orthonormal bases have no proper term for their matrix-form. to enforce
orthonormality to shared filter bases during training. For instance, when
convolution operators in each residual block group of a ResNet shares a filter
basis, the objective function $L_{R}$ can be defined to have the orthogonality
regularization term in addition to the original loss $L$:

$\centering L_{R}=L+\lambda\sum_{g}^{G}{\|W_{basis}^{(g)}{}^{\top}\cdot
W_{basis}^{(g)}-I\|^{2}},\@add@centering$ (7)

where $W_{basis}^{(g)}$ is a shared filter basis for $g$-th residual block
group and $\lambda$ is a hyperparameter.

In Equations 4 and 5, we also need to ensure that the norm of
$\frac{\partial\mathbf{x}^{j}}{\partial a^{j}}$ is bounded for stability
during forward and backward passes (Pascanu et al., 2013). It is shown that
batch normalization after non-linear activation at each convolution layer
ensures healthy norms (Ioffe & Szegedy, 2015; Guo et al., 2019; Jastrzebski et
al., 2018). In Section 4.4, we empirically show that the proposed
orthogonality regularization and batch normalization similarly improve the
flow of gradients during training.

###  3.3 Enhancing Representation Capability

Figure 2: The structure of residual blocks with a shared filter basis. For
ResNet34, two full convolution operators are factorized to uncover a shared
filter basis. For ResNets with bottleneck blocks (ResNet50/101) and
MobileNetV2, one or two computationally expensive convolution operators of
each convolution block are selected as shared filter bases and the other
convolution operators are designated as non-shared coefficients. In (a), the
filter basis has non-shared elements to further increase representation
capability.

When convolution operators share a common filter basis, they are all in the
same vector subspace. Therefore, if the rank of the filter basis is too low,
it might limit the representation capability of the convolution operators
sharing the filter basis. Conversely, if the rank of the shared filter basis
is too high (e.g., $R\geq T$), the computational gain of factorized structure
is mitigated. One way to increase the representational power of each
convolution operator, while still maintaining its computational complexity
low, is placing the convolution operators in different subspace by adding a
small number of non-shared elements to the filter basis. For instance, we
build a filter basis $W_{basis}$ by combining shared elements and non-shared
elements:

$\centering W_{basis}=W_{bs\\_shared}\cup W_{bs\\_unique},\@add@centering$ (8)

where $W_{bs\\_shared}=\\{W^{r}_{bs\\_shared}\in\mathbb{R}^{k\times k\times
S},r\in[1..n]\\}$ are shared filter basis elements, and
$W_{bs\\_unique}=\\{W^{r}_{bs\\_unique}\in\mathbb{R}^{k\times k\times
S},r\in[n{+}1..R]\\}$ are non-shared filter basis elements. For example,
Figure 2-(a) shows that the filter basis of ResNet34 is compose of both shared
and non-shared elements. One disadvantage of this hybrid scheme is that non-
shared filter basis elements require more parameters. The ratio of non-shared
basis elements can be varied to control the tradeoffs. But, our results in
Section A.1 show that providing only a few non-shared elements to a filter
basis is enough to achieve high performance.

##  4 Experiments

In this section, we perform a series of experiments on image classification
and object detection using several modern networks as base models. We also
analyze the effect of the orthogonality regularization.

For ResNets with conventional convolution filters (e.g, ResNet34), a filter is
replaced by the proposed factorized convolution block as shown in Figure
2-(a). A filter basis is shared within a residual block group having same
kernel dimensions. Throughout the experiments, we denote by ResNet$L$-S$s$U$u$
a ResNet with $L$ layers that has a filter basis with a $s$-to-$u$ ratio of
shared elements and non-shared elements.

ResNets with bottleneck blocks (e.g., ResNet50) and MobileNetV2 already have
decomposed block structures. Therefore, we designate one or two convolution
operators with the largest parameters in each block as filter bases sharing
weights across iterative blocks, as shown in Figure 2-(b) and -(c). During the
training, the proposed orthogonality regularization is applied to enforce the
elements of the shared filter bases orthonormal. In these models, we leave at
least one convolution operator as non-shared coefficients to improve
expression ability. In Section 4.1, we explore the effect of sharing both a
filter basis and coefficients in residual groups.

During experiments, all programs for training and evaluation run on PCs
equipped either with four RTX 2080Ti GPUs or two RTX 3090 GPUs and an Intel
i9-10900X CPU @3.7GHz.

###  4.1 Image Classification on ImageNet

We evaluate our method on the ILSVRC2012 dataset (Russakovsky et al., 2015)
that has 1000 classes. The dataset consists of 1.28M training and 50K
validation images. We use ResNets and MobileNetV2 as base models. We train the
ResNet-derived models for 150 epochs with SGD optimizer with a mini-batch size
of 256, a weight decay of 1e-4, and a momentum of 0.9. The learning rate
starts with 0.1 and decays by 0.1 at 60-th, 100-th, and 140-th epochs.
MobileNetV2 and our MobileNetV2-Shared models are trained for 300 epochs with
a weight decay of 1e-5. Its learning rate starts with 0.1 and decays by 0.1 at
150-th, 225-th, and 285-th epochs.

Table 1: Error (%) on ImageNet. In ResNet50/101-Shared‡ and MobileNetV2-Shared‡, first two convolution operators of each residual block are designated as recursively shared filter bases. In ResNet50-Shared-All, both a filter basis and coefficients are shared recursively across blocks of each residual group. In ResNet50-Shared-NoOrthoReg, orthogonality regulaization is not applied to ResNet50-Shared during training. Latency is measured on the Nvidia Jetson TX2 embedded board (GPU, batch size = 1).  Baseline | Model | Params | FLOPs | top-1 | top-5 | Latency
---|---|---|---|---|---|---
ResNet34 | ResNet34 (baseline) | 21.80M | 3.67G | 26.70 | 8.58 | 33.6ms
Filter Pruning (Li et al., 2017) | 19.30M | 2.76G | 27.83 | - | -
ResNet34-S48U1 (ours) | 11.79M | 3.26G | 26.67 | 8.54 | 38.6ms
ResNet50 | ResNet50 (baseline) | 25.56M | 4.11G | 23.85 | 7.13 | 43.8ms
Versatile-ResNet50 (Wang et al., 2018) | 19.0M | 3.2G | 24.5 | 7.6 | -
FSNet (Yang et al., 2020) | 13.9M | - | 26.89 | 8.63 | -
ResNet50-Shared (ours) | 20.51M | 4.11G | 23.64 | 6.98 | 43.3ms
ResNet50-Shared‡ (ours) | 18.26M | 4.11G | 23.95 | 7.14 | 43.3ms
ResNet50-Shared-All | 16.02M | 4.11G | 24.35 | 7.41 | -
ResNet50-Shared-NoOrthoReg | 20.51M | 4.11G | 24.19 | 7.34 | -
ResNet101 | ResNet101 (baseline) | 44.55M | 7.83G | 22.63 | 6.44 | 73.2ms
ResNet101-Shared (ours) | 29.47M | 7.83G | 22.31 | 6.47 | 72.9ms
MobileNetV2 | MobileNetV2 (baseline) | 3.50M | 0.33G | 28.0 | 9.71 | 18.4ms
DR-MobileNetV2 (Guo et al., 2019) | 2.96M | 0.27G | 28.2 | 9.72 | -
MobileNetV2-Shared (ours) | 3.24M | 0.33G | 27.61 | 9.34 | 17.9ms
MobileNetV2-Shared‡ (ours) | 2.98M | 0.33G | 28.21 | 9.85 | 17.8ms

In Table 1, we compare our results with competing techniques. Most parameter-
sharing and compression techniques save parameters at the cost of performance.
Unlikely, the results in Table 1 show that our approach consistently
outperforms the counterpart networks while saving a significant amount of
parameters. For example, ResNet34-S48U1 outperforms the counterpart ResNet34
while only using 54.0% parameters. Since our models derived from ResNet50/101
and MobileNetV2 already have factorized convolution blocks, their overall
parameter-saving is not as high as ResNet34-S48U1’s. However, they still save
about 19.8%, 33.9% and 7.5% parameters, respectively, while outperforming the
baselines. The benefit of our method is more pronounced in deeper networks. In
ResNet101, for example, a filter basis is shared by up to 23 recursive
convolution layers, resulting in 35.7% reduction of parameters. In Table 1,
ResNet50-Shared-NoOrhtoReg shows the effect of not applying the orthogonality
regularization on ResNet50-Shared. Our result shows that its top1 accuracy
drops to 75.81%, that is 0.34% lower than the counterpart ResNet50. In
contrast, when the orthogonality regularization is applied, the same model
achieves 76.35%, which is 0.21% higher than the counterpart ResNet50. This
0.55% improvement is obtained simply by applying orthogonality regularization.
In Table 1, ResNet50-Shared-All shows the effect of sharing not only a filter
basis but also coefficients recursively within residual groups. The result
shows that we can save further 18% parameters (reduced from 20.51M to 16.02M)
of ResNet50-Shared, but its top-1 (top-5) accuracy drops to 75.65% (92.59%),
which is 0.5% lower than the counterpart ResNet50. This result shows that
having unshared coefficients is important for high performance. Yang et al.
(2019) also demonstrates that sharing coefficients drops the accuracy by
additional 2%.

Although ResNet34-S48U1 requires lower FLOPs than the counterpart ResNet34, it
takes 14% longer latency on Jetson TX2. This overhead mostly comes from
separated convolution operations for shared and non-shared filter basis
elements. Current GPU-based deep learning libraries are not optimized to
process such separated operations efficiently. Unlikely, ResNet50/101-Shared
and MobileNetV2-Shared do not have such separated convolution operations, and,
hence, their latency is slightly lower on the device that is constrained by
limited cache and memory.

###  4.2 Image Classification on CIFAR

We evaluate the effectiveness of our method on various modern CNNs with the
CIFAR dataset that has 50,000 and 10,000 $32\times 32$ images for training and
testing, respectively. For training networks, we follow a similar training
scheme in He et al. (2016). Standardized data-augmentation and normalization
are applied to input data. Networks are trained for 300 epochs with SGD
optimizer with a weight decay of 5e-4 and a momentum of 0.9. The learning rate
is initialized to 0.1 and is decayed by 0.1 at 50% and 75% of the epochs.

Table 2: Error (%) on CIFAR-100. ‘⋆’ denotes that the orthogonality regularization is not applied. Baseline | Model | Params | FLOPs | Error
---|---|---|---|---
ResNet34 | ResNet34 (baseline) | 21.33M | 1.17G | 22.49
ResNet34-S32U1⋆ (ours) | 7.73M | 0.78G | 22.92
ResNet34-S32U1 (ours) | 7.73M | 0.78G | 21.84
DenseNet121 | DenseNet121 (baseline) | 7.05M | 0.91G | 21.95
DenseNet121-S16U1 (ours) | 5.08M | 0.72G | 22.15
ResNeXt50 | ResNeXt50 (baseline) | 23.17M | 1.36G | 20.71
ResNeXt50-S64U4 (ours) | 19.3M | 1.19G | 20.09
Table 3: Error (%) on CIFAR-10. ‘‡’ denotes having 2 shared bases in each convolution block group. ‘$\star$’ denotes that the orthogonality regularization is not applied. Baseline | Model | Params | FLOPs | Error
---|---|---|---|---
ResNet32 | ResNet32 (baseline) | 0.46M | 0.07G | 7.51
ResNet32-S16U1‡ (ours) | 0.24M | 0.08G | 6.95
ResNet56 | ResNet56 (baseline) | 0.85M | 0.16G | 6.97
Filter Pruning (Li et al., 2017) | 0.77M | 0.09G | 6.94
KSE (Li et al., 2019b) | 0.43M | 0.06G | 6.77
DR-Res 40 (Guo et al., 2019) | 0.50M | 0.11G | 6.51
ResNet56-S16U1⋆ (ours) | 0.27M | 0.15G | 7.70
| ResNet56-S16U1 (ours) | 0.27M | 0.15G | 7.46
|  ResNet56-S16U1‡ (ours) | 0.31M | 0.15G | 6.33

Table 2 shows the results on CIFAR-100. Networks trained with the proposed
method consistently outperform their counterparts. For instance,
ResNet34-S32U1 requires only 36.2% parameters and 66.6% FLOPs of the
counterpart ResNet34 while achieving lower test error (21.84%) than much
deeper ResNet50 (22.36%). To show the generality of our work, we apply the
proposed method to DenseNet (Huang et al., 2017), and ResNeXt (Xie et al.,
2017). Although the overall gain is not as pronounced as ResNets’, we still
observe reduction of resource usages in these networks. For instance,
ResNeXt50-S64U4 outperforms the counterpart ResNeXt50 while saving parameters
and FLOPs by 16.7% and 12.1%, respectively. In ResNeXt, the gain is limited
since they mainly exploit group convolutions; each group convolution is
decomposed for filter basis sharing in our network.

The result on CIFAR-10 with ResNets is presented in Table 3. Unlike networks
on CIFAR-100, networks on CIFAR-10 has much fewer channels (e.g. 16 channels
in the first residual block group) and, hence, projecting filters to such low
dimensional subspace might limit the performance of the networks. For
instance, in ResNet32-S8U1, filters are supposed to be projected onto 9
dimensional subspace consisting of 8 shared and 1 non-shared filter basis
elements. Further, for deeper networks such as ResNet56, a filter basis is
supposed to be shared by many residual blocks in the group, and it can damage
the performance. For example, every filter basis in ResNet56-S16U1 is shared
by 8 residual blocks, or 16 convolution operators. Due to this excessive
sharing, though ResNet56-S16U1 saves 41.3% parameters, its testing error
(7.46%) is higher than the counterpart ResNet56’s (6.97%).

To remedy this problem, we introduce a variant, in which each residual block
group of the networks uses 2 shared bases; one basis is shared by the first
convolution operators of recursive residual blocks, and the other basis is
shared by the second convolution operators. In Table 3, networks with a ‘‡’
mark denote this variant. Though this variant slightly increases the
parameters of the networks, it can prevent excessive sharing of parameters.
For example, although ResNet56-S16U1‡ needs 0.04M more parameters for
additional shared bases, it still saves 63% parameters of the counterpart
ResNet56 and achieves lower testing error of 6.33%.

In Table 3, we compare our results with similar state-of-the-art techniques.
Our method achieves better performance and parameter-saving than other
approaches such as filter pruning (Li et al., 2017), kernel clustering (Li et
al., 2019b), and recursive parameter sharing (Guo et al., 2019).

###  4.3 Object Detection on MS COCO

In order to explore the generalization ability of our approach, we next use
COCO 2017 dataset on object detection task using Faster R-CNN (Ren et al.,
2017), Mask R-CNN (He et al., 2017a), and RetinaNet (Lin et al., 2017) as
detectors. We compare ResNet50/101 and our ResNet50/101-Shared as backbone
networks of the detectors. These backbone networks are pre-trained on
ImageNet, then are transferred to MS COCO by fine-tuning. We use MMDetection
(Chen et al., 2019) toolbox and employ default settings for training and
evaluation. All networks are trained on train2017 for 12 epochs using SGD with
weight decay of 1e-4, momentum of 0.9 and mini-batch size of 8 (2 examples per
GPU). The learning rate is initialized to 0.01 and decays by 0.1 at 8-th and
11-th epochs.

Table 4: Object detection results on COCO 2017 validation set. ResNet50-Shared‡ uses first two convolution operators of each residual block as recursively shared filter bases.  Backbone | Detector | #Params | GFLOPs | $AP$ | $AP_{50}$ | $AP_{75}$ | $AP_{S}$ | $AP_{M}$ | $AP_{L}$
---|---|---|---|---|---|---|---|---|---
ResNet50 (baseline) | Faster R-CNN | 41.53M | 207.07 | 36.4 | 58.2 | 39.2 | 21.8 | 40.0 | 46.2
-Shared (ours) |  36.70M | 206.87 | 37.2 | 58.1 | 40.2 | 21.6 | 40.8 | 47.9
-Shared‡ (ours) |  34.46M | 206.87 | 36.6 | 57.4 | 39.8 | 21.4 | 40.1 | 47.3
ResNet101 (baseline) | 60.52M | 283.14 | 38.7 | 60.6 | 41.9 | 22.7 | 43.2 | 50.4
-Shared (ours) |  45.67M | 282.84 | 39.0 | 59.7 | 42.7 | 22.4 | 42.7 | 50.8
ResNet50 (baseline) | Mask R-CNN | 44.18M | 275.58 | 37.2 | 58.9 | 40.3 | 22.2 | 40.7 | 48.0
-Shared (ours) |  39.35M | 259.94 | 37.9 | 58.4 | 41.4 | 22.4 | 41.2 | 49.2
-Shared‡ (ours) |  37.10M | 259.94 | 37.3 | 57.6 | 40.7 | 21.5 | 40.3 | 48.5
ResNet101 (baseline) | 67.17M | 351.65 | 39.4 | 60.9 | 43.3 | 23.0 | 43.7 | 51.4
-Shared (ours) |  48.31M | 335.91 | 39.8 | 60.3 | 43.6 | 22.7 | 43.5 | 51.9
ResNet50 (baseline) | RetinaNet | 37.74M | 239.32 | 35.6 | 55.5 | 38.2 | 20.0 | 39.6 | 46.5
-Shared (ours) |  32.92M | 239.12 | 36.2 | 55.0 | 38.6 | 20.3 | 39.7 | 47.1
-Shared‡ (ours) |  30.67M | 239.12 | 35.6 | 54.1 | 38.2 | 19.8 | 39.2 | 46.9
ResNet101 (baseline) | 56.74M | 315.39 | 37.7 | 57.5 | 40.4 | 21.1 | 42.2 | 49.5
-Shared (ours) |  41.88M | 315.09 | 37.7 | 56.8 | 40.3 | 21.2 | 41.5 | 49.5
Table 5: Instance segmentation results using Mask R-CNN on COCO val2017. Backbone | $AP$ | $AP_{50}$ | $AP_{75}$ | $AP_{S}$ | $AP_{M}$ | $AP_{L}$
---|---|---|---|---|---|---
ResNet50 (baseline) | 34.1 | 55.5 | 36.2 | 16.1 | 36.7 | 50.0
-Shared (ours) | 34.5 | 55.4 | 36.9 | 18.9 | 37.7 | 46.5
ResNet101 (baseline) | 35.9 | 57.7 | 38.4 | 16.8 | 39.1 | 53.6
-Shared (ours) | 35.9 | 57.2 | 38.3 | 19.1 | 39.2 | 48.6

Table 4 shows the results on val2017 containing 5000 images. Our backbone
networks with shared filter bases consistently outperform baselines in all
detectors in terms of COCO’s standard metric AP while saving up to 28.1%
parameters. In Table 5, our models also achieve similar performance
improvement in instance segmentation using Mask R-CNN. This consistent
performance improvement is obtained simply by replacing the backbone networks
with ours and it demonstrates that the proposed recursive convolution block
design and the orthogonality regularization enables learning better feature
representations with a smaller amount of parameters.

###  4.4 Analysis: Effect of Orthogonality Regularization

To investigate the effect of the orthogonality regularization during training,
we track the flows of gradients while training ResNet34-S16U1. Figure 3 shows
the traces of the maximum and the mean gradients flowing in the filter bases
during an epoch. On every 10 iterations of a batch, the maximum and the mean
gradients are overlapped on top of the old plots. Therefore, the bars look
darker if bars are more overlapped. Jastrzebski et al. (2018) and Guo et al.
(2019) showed that unshared batch normalization (BN) strategy mitigates the
vanishing/exploding gradients problem, and our result in Figure 3-(b) shows
that unshared BNs following shared filter bases improve the flow of gradients.
When the proposed orthogonality regularization is applied to the shared filter
bases, in Figure 3-(c), similar effect is observed on gradients. When both
unshared BNs and the orthogonality regularization are applied together, in
Figure 3-(d), further stronger, but still bounded, flows of gradients are
observed. This trend is consistently observed during training. We conjecture
that this healthy flow of gradients enables learning better feature
representations during optimization process, resulting in superior performance
on various tasks.

Figure 3: The flows of gradients in 4 shared filter bases of ResNet34-S16U1 at
the same epoch. For comparison, the orthogonality regularization and the batch
normalization (BN) following the filter bases are turned on and off. In (b)
and (c), BNs and the orthogonality regularization, respectively, improve the
flow of gradients. In (d), when both BNs and orthogonality regularization are
applied together, the strongest flow of gradients is observed. This trend is
consistently observed during the training.  Figure 4: Cosine similarities of
coefficients of the 2nd and the 3rd residual block groups in ResNet34-S16U1. X
and Y axes are indexes to the coefficients of the residual block groups
sharing filter bases. Brighter colors indicate higher similarity. In (b), when
the orthogonality regularization is applied, the similarity is clearly
lowered, implying less redundancy in parameters.

To further analyze the effect of the orthogonality regularization, in Figure
4, we illustrate absolute cosine similarities of all coefficients of the 2nd
and the 3rd residual block groups of ResNet34-S16U1. The X and Y axes display
the indexes to the coefficients in the residual blocks. In Figure 4-(b), we
can clearly see that coefficients manifest lower similarities when the
orthogonality regularization is applied. In Figure 4-(a), when the
orthogonality regularization is not applied, an interesting grid pattern is
observed in Figure 4-(a). This repetitive grid pattern might be related to
ResNets’ nature of iterative feature refinement (Jastrzebski et al., 2018).
However, such high cosine similarity is directly related to the higher
redundancy in the networks. When our orthogonality regularization is applied,
such repetitive patterns are less evident, implying that the recursive
convolution layers learn better feature representations with less redundancy.

##  5 Conclusions

We introduce a recursive convolution block design and effective training
method for parameter-efficient CNNs. In this work, a common filter basis,
shared by repeating convolution layers, is learned while effectively avoiding
the vanishing/exploding gradients problem through the proposed orthogonality
regularization. Experimental results on image classification and image
detection show that our approach consistently outperforms over-parameterized
counterpart models while significantly saving parameters. This consistent
performance improvement demonstrates that the proposed approach enables
effective learning of better feature representations while a significant
amount of parameters are shared. We believe that the proposed recursive
convolution blocks and training method suggests important possibilities for
neural architecture search (NAS) to explore resource-efficient CNNs.

## Acknowledgments and Disclosure of Funding

We thank the anonymous reviewers for their constructive comments and
suggestions. This work was supported by the National Research Foundation of
Korea (NRF) Grants Funded by the Ministry of Science and ICT under Grant
NRF-2019R1F1A1060959.

## References

  * Bhalgat et al. (2020) Bhalgat, Y., Zhang, Y., Lin, J. M., and Porikli, F.  Structured convolutions for efficient neural network design.  _Advances in Neural Information Processing Systems_ , 33:5553–5564, 2020.
  * Chen et al. (2019) Chen, K., Wang, J., Pang, J., Cao, Y., Xiong, Y., Li, X., Sun, S., Feng, W., Liu, Z., Xu, J., Zhang, Z., Cheng, D., Zhu, C., Cheng, T., Zhao, Q., Li, B., Lu, X., Zhu, R., Wu, Y., Dai, J., Wang, J., Shi, J., Ouyang, W., Loy, C. C., and Lin, D.  MMDetection: Open mmlab detection toolbox and benchmark.  _arXiv preprint arXiv:1906.07155_ , 2019.
  * Denton et al. (2014) Denton, E. L., Zaremba, W., Bruna, J., LeCun, Y., and Fergus, R.  Exploiting linear structure within convolutional networks for efficient evaluation.  In _Advances in neural information processing systems_ , pp. 1269–1277, 2014.
  * Eigen et al. (2014) Eigen, D., Rolfe, J. T., Fergus, R., and LeCun, Y.  Understanding deep architectures using a recursive convolutional network.  In _2nd International Conference on Learning Representations, ICLR 2014_ , 2014.
  * Graves et al. (2013) Graves, A., Mohamed, A., and Hinton, G.  Speech recognition with deep recurrent neural networks.  In _2013 IEEE International Conference on Acoustics, Speech and Signal Processing_ , pp. 6645–6649, 2013.
  * Guo et al. (2019) Guo, Q., Yu, Z., Wu, Y., Liang, D., Qin, H., and Yan, J.  Dynamic recursive neural network.  In _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pp. 5142–5151, 2019.
  * Han et al. (2016) Han, S., Mao, H., and Dally, W. J.  Deep compression: Compressing deep neural network with pruning, trained quantization and huffman coding.  In Bengio, Y. and LeCun, Y. (eds.), _4th International Conference on Learning Representations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings_ , 2016.
  * He et al. (2016) He, K., Zhang, X., Ren, S., and Sun, J.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 770–778, 2016.
  * He et al. (2017a) He, K., Gkioxari, G., Dollar, P., and Girshick, R.  Mask r-cnn.  _2017 IEEE International Conference on Computer Vision (ICCV)_ , Oct 2017a.
  * He et al. (2017b) He, Y., Zhang, X., and Sun, J.  Channel pruning for accelerating very deep neural networks.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pp. 1389–1397, 2017b.
  * Howard et al. (2017) Howard, A. G., Zhu, M., Chen, B., Kalenichenko, D., Wang, W., Weyand, T., Andreetto, M., and Adam, H.  Mobilenets: Efficient convolutional neural networks for mobile vision applications.  _arXiv preprint arXiv:1704.04861_ , 2017.
  * Huang et al. (2017) Huang, G., Liu, Z., Van Der Maaten, L., and Weinberger, K. Q.  Densely connected convolutional networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 4700–4708, 2017.
  * Ioffe & Szegedy (2015) Ioffe, S. and Szegedy, C.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In _International Conference on Machine Learning_ , pp. 448–456, 2015.
  * Jaderberg et al. (2014) Jaderberg, M., Vedaldi, A., and Zisserman, A.  Speeding up convolutional neural networks with low rank expansions.  _arXiv preprint arXiv:1405.3866_ , 2014.
  * Jastrzebski et al. (2018) Jastrzebski, S., Arpit, D., Ballas, N., Verma, V., Che, T., and Bengio, Y.  Residual connections encourage iterative inference.  In _International Conference on Learning Representations_ , 2018.
  * Kim et al. (2016) Kim, J., Kwon Lee, J., and Mu Lee, K.  Deeply-recursive convolutional network for image super-resolution.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 1637–1645, 2016.
  * LeCun et al. (1990) LeCun, Y., Denker, J. S., and Solla, S. A.  Optimal brain damage.  In _Advances in neural information processing systems_ , pp. 598–605, 1990.
  * Li et al. (2017) Li, H., Kadav, A., Durdanovic, I., Samet, H., and Graf, H. P.  Pruning filters for efficient convnets.  In _International Conference on Learning Representations_ , 2017.
  * Li et al. (2019a) Li, Y., Gu, S., Van Gool, L., and Timofte, R.  Learning filter basis for convolutional neural network compression.  In _Proceedings of the IEEE International Conference on Computer Vision_ , 2019a.
  * Li et al. (2019b) Li, Y., Lin, S., Zhang, B., Liu, J., Doermann, D., Wu, Y., Huang, F., and Ji, R.  Exploiting kernel sparsity and entropy for interpretable cnn compression.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pp. 2800–2809, 2019b.
  * Liang & Hu (2015) Liang, M. and Hu, X.  Recurrent convolutional neural network for object recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 3367–3375, 2015.
  * Lin et al. (2017) Lin, T.-Y., Goyal, P., Girshick, R., He, K., and Dollár, P.  Focal loss for dense object detection.  In _Proceedings of the IEEE international conference on computer vision_ , 2017.
  * Pascanu et al. (2013) Pascanu, R., Mikolov, T., and Bengio, Y.  On the difficulty of training recurrent neural networks.  In _International conference on machine learning_ , pp. 1310–1318, 2013.
  * Polyak & Wolf (2015) Polyak, A. and Wolf, L.  Channel-level acceleration of deep face representations.  _IEEE Access_ , 3:2163–2175, 2015.
  * Qiu et al. (2018) Qiu, Q., Cheng, X., Sapiro, G., et al.  Dcfnet: Deep neural network with decomposed convolutional filters.  In _International Conference on Machine Learning_ , pp. 4198–4207. PMLR, 2018.
  * Ren et al. (2017) Ren, S., He, K., Girshick, R., and Sun, J.  Faster r-cnn: Towards real-time object detection with region proposal networks.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , Jun 2017.
  * Russakovsky et al. (2015) Russakovsky, O., Deng, J., Su, H., Krause, J., Satheesh, S., Ma, S., Huang, Z., Karpathy, A., Khosla, A., Bernstein, M., et al.  Imagenet large scale visual recognition challenge.  _International journal of computer vision_ , 115(3):211–252, 2015.
  * Savarese & Maire (2019) Savarese, P. and Maire, M.  Learning implicitly recurrent CNNs through parameter sharing.  In _International Conference on Learning Representations_ , 2019.
  * Socher et al. (2011) Socher, R., Lin, C. C., Manning, C., and Ng, A. Y.  Parsing natural scenes and natural language with recursive neural networks.  In _Proceedings of the 28th international conference on machine learning (ICML-11)_ , pp. 129–136, 2011.
  * Vorontsov et al. (2017) Vorontsov, E., Trabelsi, C., Kadoury, S., and Pal, C.  On orthogonality and learning recurrent networks with long term dependencies.  In _Proceedings of the 34th International Conference on Machine Learning_ , pp. 3570–3578, 2017.
  * Wang et al. (2018) Wang, Y., Xu, C., Chunjing, X., Xu, C., and Tao, D.  Learning versatile filters for efficient convolutional neural networks.  In _Advances in Neural Information Processing Systems_ , pp. 1608–1618, 2018.
  * Xie et al. (2017) Xie, S., Girshick, R., Dollár, P., Tu, Z., and He, K.  Aggregated residual transformations for deep neural networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 1492–1500, 2017.
  * Xingjian et al. (2015) Xingjian, S., Chen, Z., Wang, H., Yeung, D.-Y., Wong, W.-K., and Woo, W.-c.  Convolutional lstm network: A machine learning approach for precipitation nowcasting.  In _Advances in neural information processing systems_ , pp. 802–810, 2015.
  * Yang et al. (2020) Yang, Y., Yu, J., Jojic, N., Huan, J., and Huang, T. S.  Fsnet: Compression of deep convolutional neural networks by filter summary.  In _International Conference on Learning Representations_ , 2020.
  * Yang et al. (2019) Yang, Z., Wang, Y., Liu, C., Chen, H., Xu, C., Shi, B., Xu, C., and Xu, C.  LegoNet: Efficient convolutional neural networks with lego filters.  In Chaudhuri, K. and Salakhutdinov, R. (eds.), _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pp. 7005–7014. PMLR, 09–15 Jun 2019.
  * Zamir et al. (2017) Zamir, A. R., Wu, T.-L., Sun, L., Shen, W. B., Shi, B. E., Malik, J., and Savarese, S.  Feedback networks.  In _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pp. 1808–1817. IEEE, 2017.
  * Zhang et al. (2015) Zhang, X., Zou, J., He, K., and Sun, J.  Accelerating very deep convolutional networks for classification and detection.  _IEEE transactions on pattern analysis and machine intelligence_ , 38(10):1943–1955, 2015.

##  Appendix A Appendix

###  A.1 Effects of Ranks of Shared/Unshared Bases

Figure 5 shows test errors on CIFAR-100 as parameters and FLOPs are increased
by varying the number of shared/non-shared filter basis elements of networks.
In general, the higher performance is expected with the more parameters. We
observe that this presumption is true for shared basis elements. For instance,
when the number of shared basis elements $s$ is varied from 8 to 32, the test
error sharply decreases from 23.1% to 21.7%. However, non-shared basis
elements manifest counter-intuitive results. Although a small number of non-
shared basis elements (e.g., $u=1$) are clearly beneficial to the performance,
the higher $u$’s do not always lead to the higher performance. For instance,
when $u=4$, both ResNet34-S16U$u$ and ResNet34-S32U$u$ show the worst
performance. This result demonstrates the difficulty of training networks with
larger parameters. Further study is required for this problem.

Figure 5: Testing errors vs. the number of parameters and FLOPs on CIFAR-100.
The number of shared basis elements (s), and non-shared basis elements (u) are
varied. Using more shared basis elements results in better performance. In
contrast, using more non-shared elements does not always improve performance,
implying the difficulty of training networks with larger parameters.