# Using Shape Metrics to Describe 2D Data Points

William Franz Lamberti  Center for Public Health Genomics
Department of Biomedical Engineering
School of Data Science
University of Virginia
Charlottesville, VA, United States
<EMAIL_ADDRESS>
###### Abstract

Traditional machine learning (ML) algorithms, such as multiple regression,
require human analysts to make decisions on how to treat the data. These
decisions can make the model building process subjective and difficult to
replicate for those who did not build the model. Deep learning approaches
benefit by allowing the model to learn what features are important once the
human analyst builds the architecture. Thus, a method for automating certain
human decisions for traditional ML modeling would help to improve the
reproducibility and remove subjective aspects of the model building process.
To that end, we propose to use shape metrics to describe 2D data to help make
analyses more explainable and interpretable. The proposed approach provides a
foundation to help automate various aspects of model building in an
interpretable and explainable fashion. This is particularly important in
applications in the medical community where the ‘right to explainability’ is
crucial. We provide various simulated data sets ranging from probability
distributions, functions, and model quality control checks (such as QQ-Plots
and residual analyses from ordinary least squares) to showcase the breadth of
this approach.

###### Index Terms:

Shape Analysis, Image Measures, Explainability, Interpretability, Emerging
Applications and Systems

##  I Introduction

In the age of big data, 2D representations of multivariate data are still
commonplace and a requirement for many experiments[1, 2, 3, 4, 5]. For
example, analysts continue to simplify results to 2D representations by
analyzing various underlying properties and implicit assumptions. However,
many of these analyses are subjective and not numeric. For instance, analyzing
the residuals of multiple regression is vital to ensure that the assumptions
of the model are met[6, 7]. However, these checks are usually checked by human
analysts and are not verified using numeric measures. There are many other
aspects of multiple regression that require human decisions to be made such as
the inclusion of interaction terms. These human decisions can be a downside
when compared to deep learning approaches since deep learning models are able
to learn what features are important for a given analysis without the need for
a human analyst[8]. Having these human elements makes experiments using these
approaches less explainable and interpretable for those not involved in the
experiment.

Ensuring that a computational analysis is as explainable and interpretable as
possible is key for explainable artificial intelligence (XAI) applications.
XAI is used in various fields such as medicine and national security[9]. For
example, the European Union has passed laws ensuring a patient’s “right to
explainability”[10]. In short, this law states that if a computational model
is used to help make a diagnosis, the aspects of the computational model must
be able to be described in layman’s terms. However, XAI is not limited to the
intersection of science and policy. XAI is key for making scientific
inferences since scientists need to understand how AI models use features to
better understand the scientific phenomena. Thus, there is a fundamental need
to make computational analyses as intuitive and clear as possible. A shape
metric like area is a prime example of a metric which has a clear definition
and is a concept that is understandable to the average person[9]. Thus, we
desire to use explainable and interpretable shape metrics to describe data to
help quantify the shapes of 2D data. This paper posits that all data that
resides in 2D feature spaces can be represented as images. These images can be
analyzed by extracting various useful shape metrics.

The work that provides the foundation for the idea of analyzing 2D data using
shape metrics is eigenvale decomposition. Eigenvalues correspond to the
relative length of the axes of the data[9]. For example, if we observe a 2D
scatterplot, the major and minor axis lengths will be captured by the first
and second eigenvalues, respectively [9]. Thus, eigenvalue decomposition is
the first idea that measures the shape of data. However, this idea has not
been significantly expanded. To that end, we are providing foundational
evidence and an approach to use a variety of shape metrics to describe data in
2D feature spaces.

Figure 1: Overview of analyzing 2D data as images using shape metrics. Once
the 2D data is collected, step 1) converts the data to binary images. Step 2)
collects various shape metrics of interest to describe the images. Step 3)
analyzes the shape metrics. In this example, the differently shaped data
points are classified. However, other analyses could be performed.

We provide a new manner to analyze 2D data as images. We first convert the
data in 2D space to 2D images of the shapes. We then collect shape metrics
from the images. We lastly analyze the images for a given analysis. By
quantifying 2D data using tangible shape metrics, we make analyses of 2D data
more explainable and interpretable. An overview of our contribution is
provided in Figure 1. The code for our experiments are provided at our GitHub
link: https://github.com/billyl320/2d_shape_points.

##  II Methods and Materials

There are various simulated scenerios we will analyze: the discrimination
between different Normal distributions, the detection of outliers in QQ-Plots,
the discrimination of different 2D functions, and the analysis of multiple
regression or ordinary least squares (OLS) residual analysis of variance.

Normal Distributions: The discrimination of different Normal distributions
will help to provide a foundation by which we are able to understand how shape
metrics are useful on conceptually concrete examples. We will denote Normal or
Gaussian distributions with a mean $\mu$ and variance $\sigma^{2}$ as
$N(\mu,\sigma^{2})$. The first experiment will have two simulated Normal
distributions with a common variance, but different means. This experiment
will show that the means of the distribution are not important for
discriminating the distributions. The following two experiments will have
pairs of Normals with the same mean but differing variances. These two
experiments provide evidence that the shape metrics are useful for
discriminating Normals with differing variances.

The followup to these experiments uses various mixtures of Gaussian
distributions. This will show that different mixtures of Normals can be
discriminated using shape metrics. Examples are provided in Figure 2. Thus,
this set of four experiments using a variety of different Normal distributions
shows that the variance and the number of Gaussians mixtures present are key
for providing differently shaped data.

(a) 1 Gaussian

(b) 2 Gaussians

(c) 3 Gaussians

(d) 4 Gaussians

Figure 2: Examples of mixtures of Gaussians.

QQ-Plots and Outlier Detection: The second experiment aims to show that shape
metrics can be used to classify QQ-Plots that have outliers. QQ-Plots are a
visualization technique used to ensure that data follows a particular
distribution[11]. For OLS, we want the residuals to follow a Normal
distribution [12, 6, 7]. We simulated various QQ-Plots with no outliers (1000
random $N(0,1)$), minor outliers (990 random $N(0,1)$ and 10 random $N(3,1)$),
medium outliers (990 random $N(0,1)$ and 10 random $N(5,1)$), and major
outliers (990 random $N(0,1)$ and 10 random $N(10,1)$) to showcase the ability
of our approach to identify outliers. Examples of these QQ-Plots as images are
provided in Figures 3(a) \- 3(d), respectively.

(a) No outliers

(b) Minor outliers

(c) Medium outliers

(d) Major outliers

Figure 3: QQ-Plot examples with different levels of outliers.

Functions: The third experiment aims to show that different shaped 2D
functions can be discriminated using shape metrics. Each function had 1000
random simulated observations per resulting image. The first function was

$Y=3X+\epsilon$ (1)

where $X\sim N(0,100)$ and $\epsilon\sim N(0,1)$. This is referred to as the
“Linear” function. Note that $\sim$ represents “distributed as”. The second
function was

$Y=4\sin{X}+\epsilon$ (2)

where $X\sim N(0,100)$ and $\epsilon\sim N(0,0.25)$. This is referred to as
the “Sine” model. The third function was

$Y=X^{2}+\epsilon$ (3)

where $X\sim N(0,100)$ and $\epsilon\sim N(0,1)$. This is referred to as the
“Parabola” function. The fourth function was

$Y=X^{4}+10X^{3}-7X^{2}+\epsilon$ (4)

where $X\sim N(0,100)$ and $\epsilon\sim N(0,1)$. The is referred to as the
“polynomial” or “Poly.” model. Examples of these functions are provided in
Figures 4(a) \- 4(d), respectively.

(a) Linear

(b) Sine

(c) Parabola

(d) Poly

Figure 4: Examples of the various simulated functional models resulting
images.

OLS Residual Analysis: The fourth and last experiment corresponds to
evaluating different types of variance plots for multiple regression. The
plots should have random scatter and not display any obvious patterns. An
example of this is provided in Figure 5(a). Cone-like shape could indicates an
increase in variance as the response increases. This is common for Poisson
phenomena. Another concerning pattern looks like an almond or an eye. This
occurs when the response is a proportion from a Binomial random variable
(shorted to “Binom”). The last error discussed here is multiplicative errors
(shorted to “Multi.”), where the plot will look like a bowtie. An example of
multiplicative errors is provided in Figure 5(b). Once these patterns are
identified, transformations can be applied to the response to correct for
these errors[6]. Automating these transformations would help to make multiple
regression more standardized and reproducible.

(a) Random normal scatter

(b) Multiplicative scatter

Figure 5: Examples of a subset of the OLS residual simulated resulting images.

Converting Data to Images: The crucial step is converting the raw data into
images. This is done by converting the images using 2D histograms and then
accepting all positive signals. We will describe this more precisely using
image operator notation[13]. Thus, 2D raw data are converted into 2D images
using

$\textbf{b}[\vec{x}]=\Gamma_{>0}\mathbb{H}_{2,2}X,$ (5)

where $X$ is the input data, $\mathbb{H}_{2,2}$ converts the data into a 2D
histogram[14], and $\Gamma_{>0}$ is the threshold image operator. From here,
we would collect various shape metrics of interest.

Shape Metrics: Due to space limitations, we do not provide the exact manner in
which these shape metrics are collected. However, we will provide a brief
description of the shape metrics used in this analysis. Much of the
description of these metrics is borrowed from Lamberti[15]. Lamberti’s
extended descriptions on the shape metrics used in this analysis can be found
at various sources[16, 15, 9]. The first metrics were the shape proportions
(SP) and encircled image-histograms (EI), which were collected from the shape
proportion and encircled image-histogram (SPEI) algorithm [16]. The EI is the
black and white pixel counts of the shape after the shape is placed in the
minimum encompassing circle and then the minimum encompassing square [16]. In
other words, this is the area and the surrounding area of the shape. The SP
value is the proportion of the area of the shape relative to the sum of the
EI. The other shape metrics collected that were used in the model were the
eigenvalues of the shapes, eccentricity [13], and circularity [13, 17]. The
eigenvalues measure the major and minor axes of the shape. Eccentricity is the
ratio of the major axis over the minor one. These are calculated using the
$1^{\text{st}}$ and $2^{\text{nd}}$ eigenvalues of the shapes, respectively
[13]. Circularity measures how circular a given shape is. This results in a
total of 7 total metrics used during our analyses. The used metrics are
summarized in Table I.

TABLE I:  Table provides the metrics used in this analysis on a given image, $i$. The first column is the $q^{\text{th}}$ metric, where $q\in\\{1,2,...,7\\}$. The last column provides the number of times the classification trees used each metric. These variables make our models interpretable and explainable[9].  $\vec{m}_{q,i}$ | Metric | Counts
---|---|---
1 | White EI | 6
2 | Black EI | 0
3 | SP value | 3
4 | Eccentricity | 4
5 |  $1^{\text{st}}$ Eigenvalue | 2
6 |  $2^{\text{nd}}$ Eigenvalue | 0
7 | Circularity | 1

Modeling: Once we collected the 7 metrics, collected 100 simulated cases
(which results in 100 images) for each class for a given experiment. We used a
classification tree to discriminate between the different classes[18]. We used
80% of the data as training and 20% as validation. We used stratified random
sampling to ensure that the proportions were evenly split between the
different classes. On the training data, we used 5-fold cross-validation to
select the complexity parameter[19, 3].

##  III Results

A summary of the results are provided in Table II. Each set of experiments is
analyzed in more detail in the following sections. As a whole, these
experiments show that raw 2D data can be analyzed as images using their shape
metrics.

TABLE II: Validation accuracy measures and confidence intervals (CIs) of classification trees for different experiments. The first 4 rows correspond to the comparison of multivariate Normal distributions under different conditions. The fifth row corresponds to the QQ-Plot experiment. The sixth row corresponds to the Function experiment. The seventh row corresponds to the OLS residual analysis experiment. Note that $\mu=\begin{bmatrix}0\\\ 0\end{bmatrix}$ and $\Sigma=\begin{bmatrix}1&0\\\ 0&1\end{bmatrix}$.  Experiment | Accuracy | Accuracy 95% CI
---|---|---
$N_{1}(\begin{bmatrix}0\\\ 0\end{bmatrix},\Sigma),N_{2}(\begin{bmatrix}10\\\ 10\end{bmatrix},\Sigma)$ | 0.525 | (0.3613, 0.6849)
$N_{1}(\mu,\Sigma),N_{2}(\mu,\begin{bmatrix}1&0.9\\\ 0.9&1\end{bmatrix})$ | 1.00 | (0.9119, 1.00)
$N_{1}(\mu,\Sigma),N_{2}(\mu,\begin{bmatrix}0.001&0\\\ 0&0.001\end{bmatrix})$ | 1.00 | (0.9119, 1.00)
4 Gaussian Mixtures | 0.975 | (0.9126, 0.997)
QQ-Plots and Outlier Detection | 0.95 | (0.8769, 0.9862)
Functions | 0.9375 | (0.8601, 0.9794)
OLS Residual Analysis | 0.9375 | (0.8601, 0.9794)

Normal Distributions: This analysis shows that shape metrics are useful for
only discriminating the shape of data, not the location of the data. The first
experiment shows where shape metrics will be unhelpful as only the means
differ between the distributions. The shape of the distribution is primarily
described by the variance of the distribution.

QQ-Plots and Outlier Detection: Despite having varying levels of outliers, we
were able to discriminate between different levels of outliers on our QQ-
Plots. This shows that the shape metrics are useful in quantifying different
levels of outliers present in 2D data. This would help to quantify and detect
outliers automatically. Thus, we have evidence that shape metrics can be used
to detect the presence of outliers automatically.

Functions: Shape metrics are useful for classifying different kinds of 2D
functions from one another. This is useful for evaluating a variety of
different functional shapes. This can be used to help guide analysts on the
kind of model to utilize for a given analysis. This experiment provides
evidence for automatic function determination.

OLS Residual Analysis: Shape metrics are useful for determining the kind of
variance pattern observed in multiple regression residual plots. This helps to
remove the more subjective aspects of model evaluation and help to provide a
standard and reproducible process for modeling choices. This provides evidence
for automated OLS residual analysis using shape metrics.

Useful Features: Table I shows the number of times each of the metrics
collected were used across all of the experiments in the classification trees.
The White EI (or area) was used the most out of all of the chosen metrics.
Eccentricity was the second most used metric. Thus, White EI and Eccentricity
should be used in future analyses and applications at a minimum.

##  IV Conclusions

Shape metrics are useful for describing a variety of different 2D data
scenarios. This provides strong evidence that all 2D raw data are images and
can be analyzed as such. Thus, we can standardize traditional analyses like
multiple regression to help automate modeling decisions and make them more
explainable and interpretable to those not involved in the analysis. This
paper provides a strong foundation for analysts to analyze 2D data as images
using shape metrics.

##  V Acknowledgements

UVA Engineering Graduate Writing Lab Peer Review Group provided valuable
feedback during initial drafts of this paper. We would also like to thank the
Zang Lab for Computational Biology at the University of Virginia for their
support.

## References

  * [1] L. McInnes, J. Healy, and J. Melville, “UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction,” _arXiv:1802.03426 [cs, stat]_ , Sep. 2020, arXiv: 1802.03426. [Online]. Available: http://arxiv.org/abs/1802.03426
  * [2] L. J. Mateo, S. E. Murphy, A. Hafner, I. S. Cinquini, C. A. Walker, and A. N. Boettiger, “Visualizing DNA folding and RNA in embryos at single-cell resolution,” _Nature_ , vol. 568, no. 7750, pp. 49–54, Apr. 2019. [Online]. Available: https://www.nature.com/articles/s41586-019-1035-4
  * [3] Y. Takei, J. Yun, S. Zheng, N. Ollikainen, N. Pierson, J. White, S. Shah, J. Thomassie, S. Suo, C.-H. L. Eng, M. Guttman, G.-C. Yuan, and L. Cai, “Integrated spatial genomics reveals global architecture of single nuclei,” _Nature_ , vol. 590, no. 7845, pp. 344–350, Feb. 2021, number: 7845 Publisher: Nature Publishing Group. [Online]. Available: https://www.nature.com/articles/s41586-020-03126-2
  * [4] T. Alexandrov, “Spatial Metabolomics and Imaging Mass Spectrometry in the Age of Artificial Intelligence,” _Annual Review of Biomedical Data Science,_ , vol. 3, Jul. 2020. [Online]. Available: https://papers.ssrn.com/abstract=3658948
  * [5] H. Hotelling, “Analysis of a complex of statistical variables into principal components,” _Journal of Educational Psychology_ , vol. 24, no. 7, pp. 498–520, Oct. 1933, publisher: Warwick & York. [Online]. Available: http://mutex.gmu.edu/login?url=http://search.ebscohost.com/login.aspx?direct=true&db=pdh&AN=1934-01354-001&site=ehost-live
  * [6] W. Mendenhall and T. T. Sincich, _A Second Course in Statistics: Regression Analysis_ , 7th ed. Boston, MA: Pearson, Jan. 2011.
  * [7] N. R. Draper and H. Smith, _Applied Regression Analysis_ , third edition ed. New York: Wiley-Interscience, Apr. 1998.
  * [8] J. Gu, Z. Wang, J. Kuen, L. Ma, A. Shahroudy, B. Shuai, T. Liu, X. Wang, G. Wang, J. Cai, and T. Chen, “Recent advances in convolutional neural networks,” _Pattern Recognition_ , vol. 77, pp. 354–377, May 2018. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0031320317304120
  * [9] W. F. Lamberti, “An Overview of Explainable and Interpretable Artificial Intelligence,” in _AI Assurance: Towards Valid, Explainable, Fair, and Ethical AI_. Elsevier, 2022.
  * [10] E. Union, “Regulation (EU) 2016/679 of the European Parliament and of the Council of 27 April 2016 on the protection of natural persons with regard to the processing of personal data and on the free movement of such data, and repealing Directive 95/46/EC (General Data Protection Regulation) (Text with EEA relevance),” May 2016, legislative Body: EP, CONSIL. [Online]. Available: http://data.europa.eu/eli/reg/2016/679/oj/eng
  * [11] M. B. Wilk and R. Gnanadesikan, “Probability Plotting Methods for the Analysis of Data,” _Biometrika_ , vol. 55, no. 1, pp. 1–17, 1968. [Online]. Available: https://www.jstor.org/stable/2334448
  * [12] G. K. Bhattacharyya and R. A. Johnson, _Statistical Concepts and Methods_ , 1st ed. Wiley, 1977. [Online]. Available: https://www.wiley.com/en-us/Statistical+Concepts+and+Methods-p-9780471072041
  * [13] J. M. Kinser, _Image Operators: Image Processing in Python_ , 1st ed. Boca Raton, FL: CRC Press, Oct. 2018.
  * [14] “numpy.histogram2d — NumPy v1.21 Manual.” [Online]. Available: https://numpy.org/doc/stable/reference/generated/numpy.histogram2d.html
  * [15] W. F. Lamberti, “Classification of Synthetic Aperture Radar Images of Icebergs and Ships Using Random Forests Outperforms Convolutional Neural Networks,” in _2020 IEEE Radar Conference (RadarConf20)_ , Sep. 2020, pp. 1–6, iSSN: 2375-5318.
  * [16] ——, “Algorithms to Improve Analysis and Classification for Small Data,” Ph.D., George Mason University, United States – Virginia, 2020, iSBN: 9798557033350. [Online]. Available: http://search.proquest.com/docview/2476825035/abstract/90CC4207B46B4068PQ/1
  * [17] A. Rosenfeld, “Compact Figures in Digital Pictures,” _IEEE Transactions on Systems, Man, and Cybernetics_ , vol. SMC-4, no. 2, pp. 221–223, Mar. 1974.
  * [18] T. Therneau, B. Atkinson, B. R. p. o. t. i. R. port, and maintainer 1999-2017), “rpart: Recursive Partitioning and Regression Trees,” Apr. 2019. [Online]. Available: https://CRAN.R-project.org/package=rpart
  * [19] G. James, D. Witten, T. Hastie, and R. Tibshirani, Eds., _An introduction to statistical learning: with applications in R_ , ser. Springer texts in statistics. New York: Springer, 2013, no. 103, oCLC: ocn828488009.


homogeneous photometry extracted by us from the DESI-LIS. Our scaling
relations are constructed from all possible combinations of nine structural
parameters for late-type galaxies. The resulting scaling relations are
homogeneously fit using a BCES bisector. Both Bayesian and classical intrinsic
scatter values are then computed for each relation. The agreement between our
slope and scatter values and the literature is generally good (Table 6).

Our analysis has yielded the following data products and main results:

  * •

A robust set of homogeneously fit scaling relations from a large sample of
late-type galaxies covering a broad range of structural properties (Section
4). This includes tables of observed, classical, and Bayesian intrinsic
scatters (Table 5), and a detailed comparison with literature values (Section
5);

  * •

A code to compute Bayesian intrinsic scatters (Appendix C);

  * •

An appreciation that Bayesian intrinsic scatters are typically larger than
classical intrinsic scatters, as they account for parameter covariances and
error propagation through robust fitting techniques (Section 2). Most
relations have intrinsic scatters ranging from 60 to 90 percent of the total
observed value, with a few relations spanning greater extremes.

  * •

An appreciation that Bayesian intrinsic scatters are more robust to biases in
the data than classical methods (Appendix B);

  * •

An identification that the scatter of $R_{23.5}-L_{23.5}$, $V_{23.5}-(g-z)$,
and $R_{23.5}-j_{\rm dyn}$ scaling relations is mostly intrinsic, largely
because most observational errors slide along the scaling relations, making
some of these relations ideally suited for comparisons with galaxy formation
models;

  * •

An identification of the tightest scaling relations by intrinsic orthogonal
scatter (Section 4.2). For each structural parameter, these relations are:
$V_{23.5}-M_{*}$; $R_{23.5}-L_{23.5},j_{*},j_{\rm dyn}$;
$L_{23.5}-V_{23.5},R_{23.5},j_{\rm dyn}$; $\Sigma_{1}-V_{23.5}$;
$(g-z)-V_{23.5}$; $M_{*}-V_{23.5}$; $M_{\rm dyn}-j_{\rm dyn}$;
$j_{*}-R_{23.5}$; and $j_{\rm dyn}-V_{23.5},R_{23.5},L_{23.5},M_{\rm dyn}$,
with some parameters having multiple equally tight relations.

  * •

A revised discussion on the nature and impact of projection corrections for
galaxy structural parameters (Section 3.4)

Genuine care is required in order to achieve meaningful analyses of galaxy
scaling relations. As discussed in Section 5, the parameters associated with
any given galaxy scaling relation depend greatly on sample selection and
analysis methods. Accurate and representative structural parameters of
galaxies, the simultaneous study of multiple scaling relations, and the
derivation of robust measures of intrinsic scatter will facilitate comparisons
with galaxy formation and evolutionary models. The proper characterization of
intrinsic scaling relation scatters ought to help bridge the gap between
observations and simulations by removing one of the layers that stands between
them.

We wish to acknowledge our Queen’s University colleague Aaron Vincent for
illuminating discussions about Bayesian statistics and data analysis, as well
as the referee for insightful comments. Special thanks go to Arjun Dey and
John Moustakas for discussions about the DESI Legacy Imaging Survey, and to
the DESI team for the extensive database they provided. This research has made
use of the NASA/IPAC Extragalactic Database (NED), which is operated by the
Jet Propulsion Laboratory, California Institute of Technology, under contract
with the National Aeronautics and Space Administration. We are grateful to the
Natural Sciences and Engineering Research Council of Canada, Ontario
Government, and Queen’s University for support through various scholarships
and grants.

##  Appendix A Classical Uncertainty Functions

Here, we report the classical uncertainty functions for all of the extracted
parameters in Section 3.3. These expressions are determined via Equation 1 and
include all sources of uncertainty from Section 3 for which we could assign a
value and calculate derivatives. Table 7 includes expressions for uncertainty
for all parameters as used in our analysis.

Table 7: Classical Uncertainty Equations

Variable | Classical uncertainty
---|---
(1) | (2)
$\sigma_{\log_{10}(R_{23.5})}^{2}$ | $=\left(\frac{\sigma_{SB}}{R_{\rm obs}\Delta_{SB,R}\ln(10)}\right)^{2}+\left(\frac{\sigma_{D}}{D\ln(10)}\right)^{2}+\sigma_{C_{R_{23.5}}}^{2}$
$\sigma_{\log_{10}(L_{23.5,z})}^{2}$ | $=\left(\frac{\sigma_{m_{z}}}{2.5}\right)^{2}+\left(2\frac{\sigma_{D}}{D\ln(10)}\right)^{2}\left(\frac{\sigma_{R_{\rm obs}}\Delta_{m,R}}{2.5}\right)^{2}+\sigma_{C_{L_{23.5}}}^{2}$
$\sigma_{\log_{10}(V_{23.5})}^{2}$ | $=\frac{1}{\ln^{2}(10)}\left[\left(\frac{\sigma_{V_{\rm obs}}}{V_{\rm obs}}\right)^{2}+\left(\frac{\sigma_{V_{\rm sys}}}{V_{\rm sys}}\right)^{2}+\left(\frac{\sigma_{i}}{\tan(i)}\right)^{2}+\left(\frac{\sigma_{z}}{1+z}\right)^{2}\right]$
$\sigma_{g-z}^{2}$ | $=\left(\sigma_{m_{g}}\right)^{2}+\left(\sigma_{m_{z}}\right)^{2}+\left(\sigma_{R_{\rm obs}}\Delta_{m_{g},R}\right)^{2}+\left(\sigma_{R_{\rm obs}}\Delta_{m_{z},R}\right)^{2}+\sigma_{C_{g-z}}^{2}$
$\sigma_{\log_{10}(M_{*})}^{2}$ | $=\frac{1}{\ln^{2}(10)}\left[\left(\frac{\sigma_{L_{23.5}}}{L_{23.5}}\right)^{2}+\left(\frac{\sigma_{g-z}\Delta_{\Upsilon,g-z}}{\Upsilon}\right)+\left(\frac{\sigma_{\Upsilon}}{\Upsilon}\right)^{2}\right]$
$\sigma_{\log_{10}(\Sigma_{1})}^{2}$ | $=\left(\frac{\sigma_{m_{1,z}}}{2.5}\right)^{2}+\left(2\frac{\sigma_{D}}{D\ln(10)}\right)^{2}+\sigma_{C_{\Sigma_{1}}}^{2}+\left(\frac{\sigma_{\Upsilon_{1}}}{\Upsilon_{1}\ln(10)}\right)^{2}$
$\sigma_{\log_{10}(M_{\rm dyn})}^{2}$ | $=\left(\frac{2\sigma_{V_{23.5}}}{V_{23.5}\ln(10)}\right)^{2}+\left(\frac{\sigma_{R_{23.5}}}{R_{23.5}\ln(10)}\right)^{2}$
$\sigma_{\log_{10}(j_{*})}^{2}$ | $=\sigma_{C_{j_{*}}}^{2}+\sigma_{\log_{10}(M_{*})}^{2}+\left(\frac{2\pi}{j_{*}\ln(10)}\right)^{2}\left[\left(\sigma_{I}\int_{0}^{R_{23.5}}\Upsilon(r)V(r)r^{2}dr\right)^{2}+\left(\sigma_{\Upsilon}\int_{0}^{R_{23.5}}I(r)V(r)r^{2}dr\right)^{2}\right.$
|
$\left.+\left(\sigma_{V_{23.5}}\int_{0}^{R_{23.5}}I(r)\Upsilon(r)r^{2})dr\right)^{2}+\left(\sigma_{R_{23.5}}I(R_{23.5})\Upsilon(R_{23.5})V(R_{23.5})R_{23.5}^{2}\right)^{2}\right]$
$\sigma_{\log_{10}(j_{\rm dyn})}^{2}$ | $=\sigma_{\log_{10}(M_{\rm dyn})}^{2}+\left(\frac{\sigma_{V_{23.5}}}{j_{\rm dyn}}\int_{0}^{R_{23.5}}\frac{3V(r)^{2}r}{G\ln(10)}dr\right)^{2}+\left(\frac{\sigma_{R_{23.5}}V(R_{23.5})^{3}R_{23.5}}{j_{\rm dyn}G\ln(10)}\right)^{2}$

Note. — Expressions for classical uncertainty propagation as described in
Equation 1. These are computed for all points in each scaling relation and
used for the classical uncertainty propagation in Section 2.2. The variables
are described in the text and Section 3.3.

Some variable names appear in Table 7 for the first time and so we describe
them here: $\sigma_{D}$ is the uncertainty on the distance to a galaxy and
appears in many of the relations.
$\sigma_{C_{R_{23.5}}},\sigma_{C_{L_{23.5}}},\sigma_{C_{g-z}},\sigma_{C_{g-z}},\sigma_{C_{\Sigma_{1}}},$
and $\sigma_{C_{j_{*}}}$ are the uncertainty on inclination correction factors
for the corresponding variables. $\sigma_{SB}$ is the uncertainty in the
surface brightness at the selected radius ($R_{\rm obs}$) and $\Delta_{SB,R}$
is the local slope of the SB profile at $R_{\rm obs}$. $\sigma_{m_{z}}$ is the
uncertainty on the apparent $z$ band magnitude, and $\Delta_{m,R}$ is the
slope of the curve of growth at the isophotal radius $R_{\rm obs}$.
$\sigma_{i}$ is the inclination uncertainty (propagation through Equation 5
not shown for clarity), and $\sigma_{z}$ is the redshift uncertainty.
$\sigma_{V_{\rm obs}}$ is the uncertainty on the individual velocity
measurement at $R_{23.5}$, or the closest point to it. The uncertainty on the
systematic velocity $\sigma_{V_{\rm sys}}$ was determined by fitting rotation
curves with alternative methods (such as an arctan model) and finding typical
agreement to within $3\,\text{km s${}^{-1}$}$. This value is used as the
systematic velocity uncertainty for all galaxies. $\Delta_{m_{g},R}$ and
$\Delta_{m_{z},R}$ are the slopes of the $g-$ and $z-$band growth curves,
respectively. $\Delta_{\Upsilon,g-z}$ is the slope of the mass-to-light ratio
as a function of colour, and $\sigma_{\Upsilon}$ is the uncertainty on the
mass-to-light ratio within $R_{\rm obs}$ set to $0.05\,\text{dex}$ (Sec 4.1 of
Roediger & Courteau 2015). Similarly, $\sigma_{\Upsilon_{1}}$ is the
uncertainty on the mass-to-light ratio within $1\,\text{kpc}$.
$\sigma_{m_{1,z}}$ is the uncertainty on the magnitude within $1\,\text{kpc}$.
$\sigma_{I}$ is the uncertainty on the intensity at $R_{23.5}$.

##  Appendix B Testing Data Truncations

Various alterations of the toy model can be made to test the effects of
realistic biases in the data. Truncating the data on an axis can represent
sampling biases. We test the effect of truncating extreme values on the
$x$-axis by setting a lower bound to the axis that removes a given percentage
of the data points. Truncating based on uncertainty can represent data quality
cuts, and we test the effect of setting an uncertainty threshold that removes
some percentage of the data. Truncating based on the residual after fitting a
relation can represent sigma clipping, and we test the effect of removing data
points with large residuals as well. Shown in Figure 6 are tests for the
stability of each method against five types of bias in the data; each test is
performed using 1000 trials of the toy model. We compare the relative error
between the predicted intrinsic scatter value and the true value as a function
of the degree of data bias.

Figure 6: Tests on intrinsic scatter estimates, showing the relative error
between each method and the true mock intrinsic scatter value averaged over
1000 trials. Shown are the relative errors for four levels of data
manipulation as given in Table 8; zero indicates no modification to the data.
The Bayesian method consistently performs better than the classical method and
is more robust to data manipulation. Table 8: Toy Model Modifications

Bias | 1 | 2 | 3 | 4
---|---|---|---|---
(1) | (2) | (3) | (4) | (5)
quality cuts | 0.05 | 0.10 | 0.25 | 0.50
sampling bias | 0.05 | 0.10 | 0.25 | 0.50
sigma clipping | 0.01 | 0.02 | 0.03 | 0.05
scaled | 0.02 | 0.05 | 0.10 | 0.25
covariance | 0.01 | 0.02 | 0.05 | 0.10

Note. — Table of data modifications from Figure 6 to test the Bayesian and
classical intrinsic scatter methods. Column (1) indicates the type of bias
applied to the data. Columns (2) - (5) indicate the degree of bias. The
numbers in the table have a different meaning for each bias type. Quality
cuts, sampling bias, and sigma clipping values refer to the fraction of
removed data points. Scaled refers to the factor by which the mock observation
uncertainties are scaled relative to their true value. Covariance refers to
the degree of extra covariant scatter introduced to the mock observations
(note that the toy model intrinsic scatter is approximately 0.13 in the same
arbitrary units).

Figure 6 makes clear that scatter analyses are adversely impacted by sigma-
clipping methods, since their accuracy and uncertainty bounds are invalidated.
This effect becomes significant for the both classical and Bayesian algorithms
once as little as $1-2\,\text{\%}$ of the data are removed. Indeed, sigma
clipping eliminates some of the most informative data points from the sample;
points close to the relation are consistent with any intrinsic scatter, while
those far from the relation that have small uncertainty can only be explained
with a nonzero intrinsic scatter. Figure 6 also shows that truncations based
on an axis or on uncertainty (as are used in Section 3.6) have essentially no
effect on either method. We tested truncating up to $50\,\text{\%}$ of the
data, and the relative error only increased slightly as would be expected with
less data to analyze. Finally, we see that the Bayesian method performs better
than the classical method in the ideal toy model (zero bias) by about
$25\,\text{\%}$. The difference is even more pronounced if there is some bias
to the data such as correlated errors.

More sinister is the possibility that the observational uncertainties are
themselves incorrect. Perhaps a given uncertainty is overestimated as the
observer wished to present conservative values, or the uncertainty could be
underestimated due to an unknown factor influencing the data. Whatever the
cause, biased uncertainties can potentially present a significant barrier to
measuring accurate intrinsic scatter values. We test the effect of modifying
the uncertainty quality. We also test the impact of correlated uncertainties
and verify the statement in Section 2 that the Bayesian algorithm properly
handles correlated uncertainties.

To test for incorrect observational uncertainties, we randomly selected
$50\,\text{\%}$ of the data points and scaled their uncertainty by some
factor, thus considering both inconsistent data (only half the data is
modified) and incorrect uncertainty values. This scaling is done after each
data point is perturbed by its uncertainty, so the algorithms will only have
access to the scaled (incorrect) values.

To test for correlated uncertainties, an extra uncertainty $\sigma_{z}$ is
added to simultaneously scatter the data in both axes. The classical method
considers each axis separately and so the $\sigma_{z}$ error will be added in
quadrature to both axes, while the Bayesian technique includes the new
uncertainty in its sampling algorithm as a covariant term.

Figure 6 also shows that both algorithms can handle a small-scale manipulation
to the uncertainty values; however the Bayesian algorithm estimate and its
uncertainty are consistently more robust. At $10\,\text{\%}$ scaling, the
classical algorithm produces an essentially meaningless uncertainty range. The
test for covariance in Figure 6 demonstrates, as expected, that the classical
method is negatively impacted by this type of uncertainty while the Bayesian
method handles it smoothly.

##  Appendix C Bayesian Intrinsic Scatter Code

Here we present the Bayesian intrinsic scatter code. The Python code below
relies on the public numpy (Harris et al. 2020) and scipy (Jones et al. 2001)
packages. The code defines the function BayesianIntrinsicScatter which has
several arguments described below.

The argument phi_list is a list object where each element contains all of the
measurement information for a galaxy (rotation curve, surface brightness
profile, distance, etc.) and the associated uncertainties. The galaxy object
(elements of the phi_list list) may be formatted in any way that is easiest
for the user as the code does not directly interact with them, instead only
passing them to other functions. sigma_max is the maximum possible intrinsic
scatter value for the prior; typically the total scatter of the scaling
relation is used. The two X and Y arguments are functions that take a phi_list
element and evaluate the $x$/$y$ axes of the scaling relation in question.
These functions can return None to indicate that a point has exceeded a data
quality cut and should be ignored. Other than returning None, no other
assumptions are made about the output of X and Y as they are only passed to
user-defined functions. sample_phi_params is a function that takes no
arguments and returns any parameters that are needed for the sample_phi
function. sample_phi takes parameters from sample_phi_params and an element
from phi_list and returns the element resampled about its uncertainty values.
For example, the distance measurement would be resampled by a normal
distribution about the original measurement with the uncertainty as the
standard deviation. relation_f is the function for the scaling relation; it
takes as arguments a tuple of parameters, an $x$-axis value, and a $y$-axis
value, then returns the residual. relation_f_fit takes a list of $x$-axis
values and a list of $y$-axis values and returns the tuple of parameters
needed for relation_f. N_samples is the number of times to sample the
sample_phi function for each galaxy and should be at least 500 for the
posterior to converge. N_sigma is the number of points to evaluate the
intrinsic scatter pdf. sigma_min is the minimum value at which to evaluate the
intrinsic scatter pdf. nprocs is the number of processors used for the
calculations; this can accelerate the Bayesian intrinsic scatter measurement.
If the speed up is not needed, one can replace the instances of pool.map with
map for the same functionality. min_pass is the minimum number of evaluations
out of N_samples that do not return None required for a galaxy to be included
in the intrinsic scatter calculation.

While the Bayesian intrinsic scatter calculation is somewhat more complicated
than a classical intrinsic scatter analysis, there are some distinct
advantages. From a coding perspective, there are only a few functions that
need to be defined. The X and Y functions are effectively already needed to
construct the scaling relation in the first place. The sample_phi function is
relatively simple to construct as it typically only involves sampling normal
distributions. Again the relation_f and relation_f_fit functions must already
be constructed in order to fit the scaling relation. Contrast this with the
classical error analysis where derivatives must be computed for each axis as a
function of each input value that has an uncertainty. Once implemented for the
first time, the Bayesian intrinsic scatter calculation is actually easier than
the classical analysis, though it does take longer to compute.

    import numpy as np
    from scipy.integrate import trapz
    from scipy.stats import norm
    from functools import partial
    from multiprocessing import Pool

    def _pdf(r,s):
        return np.sum(norm.pdf(r, loc = 0, scale = s))

    def BayesianIntrinsicScatter(phi_list, sigma_max, X, Y,
                                 sample_phi_params, sample_phi,
                                 relation_f, relation_f_fit,
                                 N_samples = 500, N_sigma = 100,
                                 sigma_min = 0.01, nprocs = 4,
                                 min_pass = 100):
        pool = Pool(nprocs)
        # array of values at which to evaluate the intrinsic scatter pdf
        S = np.linspace(sigma_min, sigma_max, N_sigma)

        residuals = [[] for n in range(len(phi_list))]
        for r in range(N_samples):
            params = sample_phi_params()
            # Resample galaxy list
            sample = pool.map(partial(sample_phi, params), phi_list)
            # Evaluate scaling relation axes
            XY = zip(pool.map(X,sample), pool.map(Y,sample),
                     range(len(phi_list)))
            XY = list(filter(lambda xy: not None in xy, XY))
            # Fit scaling relation
            fit = relation_f_fit(list(xy[0] for xy in XY),
                                 list(xy[1] for xy in XY))
            for xy in XY:
                # Store scaling relation residual
                residuals[xy[2]].append(relation_f(fit, xy[0], xy[1]))

        posteriors = []
        for r in filter(lambda r: len(r) > min_pass, residuals):
            # evaluate the intrinsic scatter pdf for this galaxy
            pdf = np.array(pool.map(partial(_pdf, r), S))

            # normalize pdf to integral 1
            posteriors.append(np.log10(pdf/trapz(pdf, S)))

        # Take product of single galaxy posteriors and normalize
        P_sigmai = np.sum(posteriors, axis = 0)
        P_sigmai -= np.max(P_sigmai)
        P_sigmai = (10**P_sigmai)/trapz(10**P_sigmai, S)

        return S, P_sigmai

## References

  * Abell et al. (1989) Abell, G. O., Corwin, Harold G., J., & Olowin, R. P. 1989, ApJS, 70, 1
  * Akritas & Bershady (1996) Akritas, M. G. & Bershady, M. A. 1996, ApJ, 470, 706
  * Andreon & Hurn (2013) Andreon, S. & Hurn, M. 2013, Statistical Analysis and Data Mining: The ASA Data Science Journal, 9, 15
  * Aquino-Ortíz et al. (2020) Aquino-Ortíz, E., Sánchez, S. F., et al. 2020, ApJ, 900, 109
  * Arora et al. (in prep.) Arora, N., Stone, C., et al. in prep., in prep.
  * Avila-Reese et al. (2008) Avila-Reese, V., Zavala, J., et al. 2008, AJ, 136, 1340
  * Barro et al. (2017) Barro, G., Faber, S. M., et al. 2017, ApJ, 840, 47
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57
  * Beifiori et al. (2012) Beifiori, A., Courteau, S., et al. 2012, MNRAS, 419, 2497
  * Bender et al. (1993) Bender, R., Burstein, D., & Faber, S. M. 1993, ApJ, 411, 153
  * Bernstein et al. (1994) Bernstein, G. M., Guhathakurta, P., et al. 1994, AJ, 107, 1962
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
  * Bottrell et al. (2017) Bottrell, C., Torrey, P., et al. 2017, MNRAS, 467, 2879
  * Bradford et al. (2016) Bradford, J. D., Geha, M. C., & van den Bosch, F. C. 2016, ApJ, 832, 11
  * Bradley et al. (2020) Bradley, L., Sipőcz, B., et al. 2020
  * Brook et al. (2012a) Brook, C. B., Stinson, G., et al. 2012a, MNRAS, 419, 771
  * Brook et al. (2012b) Brook, C. B., Stinson, G., et al. 2012b, MNRAS, 424, 1275
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
  * Bundy et al. (2015) Bundy, K., Bershady, M. A., et al. 2015, ApJ, 798, 7
  * Burstein et al. (1991) Burstein, D., Haynes, M. P., & Faber, M. 1991, Nature, 353, 515
  * Byun et al. (1994) Byun, Y. I., Freeman, K. C., & Kylafis, N. D. 1994, ApJ, 432, 114
  * Chen et al. (2020a) Chen, G., Zhang, H.-X., et al. 2020a, ApJ, 895, 146
  * Chen et al. (2020b) Chen, Z., Faber, S. M., et al. 2020b, ApJ, 897, 102
  * Cheung et al. (2012) Cheung, E., Faber, S. M., et al. 2012, ApJ, 760, 131
  * Cho & Park (2009) Cho, J. & Park, C. 2009, ApJ, 693, 1045
  * Conroy (2013) Conroy, C. 2013, ARA&A, 51, 393
  * Cooley & Tukey (1965) Cooley, J. W. & Tukey, J. W. 1965, Mathematics of computation, 19, 19
  * Cortese et al. (2016) Cortese, L., Fogarty, L. M. R., et al. 2016, MNRAS, 463, 170
  * Courteau (1997) Courteau, S. 1997, AJ, 114, 2402
  * Courteau et al. (2007) Courteau, S., Dutton, A. A., et al. 2007, ApJ, 671, 203
  * Courteau & Rix (1999) Courteau, S. & Rix, H.-W. 1999, ApJ, 513, 561
  * Courteau et al. (2000) Courteau, S., Willick, J. A., et al. 2000, ApJ, 544, 636
  * Crain et al. (2015) Crain, R. A., Schaye, J., et al. 2015, MNRAS, 450, 1937
  * Crampin & Hoyle (1964) Crampin, D. J. & Hoyle, F. 1964, ApJ, 140, 99
  * Dale et al. (1999) Dale, D. A., Giovanelli, R., et al. 1999, AJ, 118, 1489
  * Devour & Bell (2016) Devour, B. M. & Bell, E. F. 2016, MNRAS, 459, 2054
  * Devour & Bell (2017) Devour, B. M. & Bell, E. F. 2017, MNRAS, 468, L31
  * Devour & Bell (2019) Devour, B. M. & Bell, E. F. 2019, ApJS, 244, 3
  * Dey et al. (2019) Dey, A., Schlegel, D. J., et al. 2019, AJ, 157, 168
  * Driver et al. (2008) Driver, S. P., Popescu, C. C., et al. 2008, ApJ, 678, L101
  * Driver et al. (2007) Driver, S. P., Popescu, C. C., et al. 2007, MNRAS, 379, 1022
  * Dutton et al. (2010) Dutton, A. A., Conroy, C., et al. 2010, MNRAS, 407, 407
  * Dutton et al. (2011) Dutton, A. A., van den Bosch, F. C., et al. 2011, MNRAS, 410, 1660
  * Efron (1992) Efron, B. 1992, in Breakthroughs in statistics (Springer), 569–593
  * Ellison et al. (2020) Ellison, S. L., Thorp, M. D., et al. 2020, MNRAS, 493, L39
  * Fall & Romanowsky (2018) Fall, S. M. & Romanowsky, A. J. 2018, ApJ, 868, 133
  * Fang et al. (2013) Fang, J. J., Faber, S. M., et al. 2013, ApJ, 776, 63
  * Fernández Lorenzo et al. (2013) Fernández Lorenzo, M., Sulentic, J., et al. 2013, MNRAS, 434, 325
  * Ferrero et al. (2020) Ferrero, I., Navarro, J. F., et al. 2020, arXiv e-prints, arXiv:2009.03916
  * Ferrero et al. (2017) Ferrero, I., Navarro, J. F., et al. 2017, MNRAS, 464, 4736
  * Friedman et al. (2010) Friedman, J., Hastie, T., & Tibshirani, R. 2010, Journal of statistical software, 33, 33
  * García-Benito et al. (2019) García-Benito, R., González Delgado, R. M., et al. 2019, A&A, 621, A120
  * Giovanelli (2013) Giovanelli, R. 2013, Advancing the Physics of Cosmic Distances, 289, 296
  * Giovanelli et al. (1994) Giovanelli, R., Haynes, M. P., et al. 1994, AJ, 107, 2036
  * Gnedin et al. (2007) Gnedin, O. Y., Weinberg, D. H., et al. 2007, ApJ, 671, 1115
  * Graham & Worley (2008) Graham, A. W. & Worley, C. C. 2008, MNRAS, 388, 1708
  * Hall et al. (2012) Hall, M., Courteau, S., et al. 2012, MNRAS, 425, 2741
  * Han (1992) Han, M. 1992, ApJ, 391, 617
  * Harris et al. (2020) Harris, C. R., Millman, K. J., et al. 2020, Nature, 585, 585
  * Haynes & Giovanelli (1984) Haynes, M. P. & Giovanelli, R. 1984, AJ, 89, 758
  * Holmberg (1958) Holmberg, E. 1958, Meddelanden fran Lunds Astronomiska Observatorium Serie II, 136, 1
  * Hopkins et al. (2009) Hopkins, P. F., Hernquist, L., et al. 2009, ApJ, 691, 1424
  * Hubble (1926) Hubble, E. P. 1926, ApJ, 64
  * Isobe et al. (1990) Isobe, T., Feigelson, E. D., et al. 1990, ApJ, 364, 104
  * Jacoby et al. (1992) Jacoby, G. H., Branch, D., et al. 1992, PASP, 104, 599
  * Jones et al. (2001) Jones, E., Oliphant, T., et al. 2001, SciPy: Open source scientific tools for Python, [Online; accessed Jan 2021]
  * Jonsson et al. (2010) Jonsson, P., Groves, B. A., & Cox, T. J. 2010, MNRAS, 403, 17
  * Kelly (2007) Kelly, B. C. 2007, ApJ, 665, 1489
  * Kim et al. (2014) Kim, J.-h., Abel, T., et al. 2014, ApJS, 210, 14
  * Knebe et al. (2018) Knebe, A., Pearce, F. R., et al. 2018, MNRAS, 475, 2936
  * Kourkchi et al. (2020) Kourkchi, E., Tully, R. B., et al. 2020, ApJ, 902, 145
  * Kourkchi et al. (2019) Kourkchi, E., Tully, R. B., et al. 2019, ApJ, 884, 82
  * Lagos et al. (2018) Lagos, C. d. P., Tobar, R. J., et al. 2018, MNRAS, 481, 3573
  * Lambas et al. (1992) Lambas, D. G., Maddox, S. J., & Loveday, J. 1992, MNRAS, 258, 404
  * Lange et al. (2015) Lange, R., Driver, S. P., et al. 2015, MNRAS, 447, 2603
  * Lapi et al. (2018) Lapi, A., Salucci, P., & Danese, L. 2018, ApJ, 859, 2
  * Lelli et al. (2016) Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016, AJ, 152, 157
  * Lelli et al. (2017) Lelli, F., McGaugh, S. S., et al. 2017, ApJ, 836, 152
  * Maller et al. (2009) Maller, A. H., Berlind, A. A., et al. 2009, ApJ, 691, 394
  * Mancera Piña et al. (2021) Mancera Piña, P. E., Posti, L., et al. 2021, A&A, 647, A76
  * Masters et al. (2010) Masters, K. L., Nichol, R., et al. 2010, MNRAS, 404, 792
  * Mathewson & Ford (1996) Mathewson, D. S. & Ford, V. L. 1996, ApJS, 107, 97
  * Mathewson et al. (1992) Mathewson, D. S., Ford, V. L., & Buchhorn, M. 1992, ApJS, 81, 413
  * McGaugh (2004) McGaugh, S. S. 2004, ApJ, 609, 652
  * McGaugh (2005) McGaugh, S. S. 2005, ApJ, 632, 859
  * McGaugh et al. (2016) McGaugh, S. S., Lelli, F., & Schombert, J. M. 2016, Phys. Rev. Lett., 117, 201101
  * McGaugh et al. (2000) McGaugh, S. S., Schombert, J. M., et al. 2000, ApJ, 533, L99
  * Mo et al. (1998) Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319
  * Möllenhoff et al. (2006) Möllenhoff, C., Popescu, C. C., & Tuffs, R. J. 2006, A&A, 456, 941
  * Mosenkov et al. (2015) Mosenkov, A. V., Sotnikova, N. Y., et al. 2015, MNRAS, 451, 2376
  * Moster et al. (2010) Moster, B. P., Somerville, R. S., et al. 2010, ApJ, 710, 903
  * Mowla et al. (2019) Mowla, L. A., van Dokkum, P., et al. 2019, ApJ, 880, 57
  * Obreschkow & Glazebrook (2014) Obreschkow, D. & Glazebrook, K. 2014, ApJ, 784, 26
  * Ouellette et al. (2017) Ouellette, N. N.-Q., Courteau, S., et al. 2017, ApJ, 843, 74
  * Peebles (1969) Peebles, P. J. E. 1969, ApJ, 155, 393
  * Peng et al. (2010) Peng, Y.-j., Lilly, S. J., et al. 2010, ApJ, 721, 193
  * Pierce & Tully (1988) Pierce, M. J. & Tully, R. B. 1988, ApJ, 330, 579
  * Pierce & Tully (1992) Pierce, M. J. & Tully, R. B. 1992, ApJ, 387, 47
  * Pizagno et al. (2005) Pizagno, J., Prada, F., et al. 2005, ApJ, 633, 844
  * Pizagno et al. (2007) Pizagno, J., Prada, F., et al. 2007, AJ, 134, 945
  * Ponomareva et al. (2017) Ponomareva, A. A., Verheijen, M. A. W., et al. 2017, MNRAS, 469, 2387
  * Posti et al. (2018) Posti, L., Fraternali, F., et al. 2018, A&A, 612, L6
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Reyes et al. (2011) Reyes, R., Mandelbaum, R., et al. 2011, MNRAS, 417, 2347
  * Rhee (2000) Rhee, M. H. 2000, in Astronomical Society of the Pacific Conference Series, Vol. 201, Cosmic Flows Workshop, ed. S. Courteau & J. Willick, 134
  * Robotham & Obreschkow (2015) Robotham, A. S. G. & Obreschkow, D. 2015, PASA, 32, e033
  * Roediger & Courteau (2015) Roediger, J. C. & Courteau, S. 2015, MNRAS, 452, 3209
  * Romanowsky & Fall (2012) Romanowsky, A. J. & Fall, S. M. 2012, ApJS, 203, 17
  * Rong et al. (2018) Rong, Y., Li, H., et al. 2018, MNRAS, 477, 230
  * Saintonge & Spekkens (2011) Saintonge, A. & Spekkens, K. 2011, ApJ, 726, 77
  * Sakai et al. (2000) Sakai, S., Mould, J. R., et al. 2000, ApJ, 529, 698
  * Scannapieco et al. (2012) Scannapieco, C., Wadepuhl, M., et al. 2012, MNRAS, 423, 1726
  * Schaller et al. (2015) Schaller, M., Dalla Vecchia, C., et al. 2015, MNRAS, 454, 2277
  * Schaye (2004) Schaye, J. 2004, ApJ, 609, 667
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525
  * Sereno & Ettori (2015) Sereno, M. & Ettori, S. 2015, MNRAS, 450, 3633
  * Shai Shalev-Shwartz (2014) Shai Shalev-Shwartz, S. B.-D. 2014, Understanding Machine Learning: From Theory to Algorithms (Cambrige University Press)
  * Shao et al. (2007) Shao, Z., Xiao, Q., et al. 2007, ApJ, 659, 1159
  * Shen et al. (2003) Shen, S., Mo, H. J., et al. 2003, MNRAS, 343, 978
  * Snyder et al. (2015) Snyder, G. F., Torrey, P., et al. 2015, MNRAS, 454, 1886
  * Springob et al. (2007) Springob, C. M., Masters, K. L., et al. 2007, ApJS, 172, 599
  * Steinmetz & Navarro (1999) Steinmetz, M. & Navarro, J. F. 1999, ApJ, 513, 555
  * Stetson (1987) Stetson, P. B. 1987, PASP, 99, 191
  * Stone & Courteau (2019) Stone, C. & Courteau, S. 2019, ApJ, 882, 6
  * Strauss & Willick (1995) Strauss, M. A. & Willick, J. A. 1995, Phys. Rep., 261, 271
  * Sweet et al. (2018) Sweet, S. M., Fisher, D., et al. 2018, ApJ, 860, 37
  * Szomoru et al. (2012) Szomoru, D., Franx, M., & van Dokkum, P. G. 2012, ApJ, 749, 121
  * Takase & Kinoshita (1967) Takase, B. & Kinoshita, H. 1967, PASJ, 19, 409
  * Taylor et al. (2011) Taylor, E. N., Hopkins, A. M., et al. 2011, MNRAS, 418, 1587
  * Teimoorinia et al. (2016) Teimoorinia, H., Bluck, A. F. L., & Ellison, S. L. 2016, MNRAS, 457, 2086
  * Tempel et al. (2010) Tempel, E., Tamm, A., & Tenjes, P. 2010, A&A, 509, A91
  * Tibshirani (1996) Tibshirani, R. 1996, Journal of the Royal Statistical Society: Series B (Methodological), 58, 58
  * Tody (1986) Tody, D. 1986, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 627, Instrumentation in astronomy VI, ed. D. L. Crawford, 733
  * Torrey et al. (2015) Torrey, P., Snyder, G. F., et al. 2015, MNRAS, 447, 2753
  * Trujillo et al. (2020a) Trujillo, I., Chamba, N., & Knapen, J. H. 2020a, MNRAS, 493, 87
  * Trujillo et al. (2020b) Trujillo, I., Chamba, N., & Knapen, J. H. 2020b, MNRAS, 495, 3777
  * Trujillo-Gomez et al. (2011) Trujillo-Gomez, S., Klypin, A., et al. 2011, ApJ, 742, 16
  * Tully & Fisher (1977) Tully, R. B. & Fisher, J. R. 1977, A&A, 54, 661
  * Tully & Fouqué (1985) Tully, R. B. & Fouqué, P. 1985, ApJS, 58, 67
  * Tully et al. (1998) Tully, R. B., Pierce, M. J., et al. 1998, AJ, 115, 2264
  * Unterborn & Ryden (2008) Unterborn, C. T. & Ryden, B. S. 2008, ApJ, 687, 976
  * van der Wel et al. (2014) van der Wel, A., Franx, M., et al. 2014, ApJ, 788, 28
  * Verheijen (2001) Verheijen, M. A. W. 2001, ApJ, 563, 694
  * Willick (1999) Willick, J. A. 1999, in Formation of Structure in the Universe, 213
  * Willick et al. (1995) Willick, J. A., Courteau, S., et al. 1995, ApJ, 446, 12
  * Willick et al. (1997) Willick, J. A., Strauss, M. A., et al. 1997, ApJ, 486, 629
  * Woo et al. (2008) Woo, J., Courteau, S., & Dekel, A. 2008, MNRAS, 390, 1453
  * Woo & Ellison (2019) Woo, J. & Ellison, S. L. 2019, MNRAS, 487, 1927
  * Wu et al. (2020) Wu, P.-F., van der Wel, A., et al. 2020, ApJ, 888, 77
  * Xiao et al. (2012) Xiao, T., Wang, T., et al. 2012, MNRAS, 421, 486
  * Yip et al. (2010) Yip, C.-W., Szalay, A. S., et al. 2010, ApJ, 709, 780
  * Yung et al. (2019) Yung, L. Y. A., Somerville, R. S., et al. 2019, MNRAS, 490, 2855
  * Zasov & Rubtsova (1989) Zasov, A. V. & Rubtsova, T. V. 1989, Soviet Astronomy Letters, 15, 51
  * Zhang et al. (2017) Zhang, H.-X., Puzia, T. H., & Weisz, D. R. 2017, ApJS, 233, 13
  * Zolotov et al. (2015) Zolotov, A., Dekel, A., et al. 2015, MNRAS, 450, 2327
  * Zou & Hastie (2005) Zou, H. & Hastie, T. 2005, Journal of the royal statistical society: series B (statistical methodology), 67, 67
  * Zwaan et al. (1995) Zwaan, M. A., van der Hulst, J. M., et al. 1995, MNRAS, 273, L35


# Continual Learning with Neuron Activation Importance

Sohee Kim Seungkyu Lee
Kyunghee University, Department of Computer Engineering,
Yongin, Republic of Korea
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Continual learning is a concept of online learning with multiple sequential
tasks. One of the critical barriers of continual learning is that a network
should learn a new task keeping the knowledge of old tasks without access to
any data of the old tasks. In this paper, we propose a neuron activation
importance-based regularization method for stable continual learning
regardless of the order of tasks. We conduct comprehensive experiments on
existing benchmark data sets to evaluate not just the stability and plasticity
of our method with improved classification accuracy also the robustness of the
performance along the changes of task order.

##  1 Introduction

Continual learning is a sequential learning scheme on multiple different
tasks. New tasks do not necessarily consist of only existing classes of
previous tasks nor statistically similar instances of existing classes. In
challenging situations, new tasks may consist of mutually disjoint classes or
existing classes with unseen types of instances in previous tasks. One of the
main challenges is learning such new tasks without catastrophic forgetting
existing knowledge of previous tasks. Researchers have proposed diverse
continual learning approaches to achieve both stability (remembering past
tasks) and plasticity (adapting to new tasks) of their deep neural networks
from sequential tasks of irregular composition of classes and varying
characteristics of training instances. Since the training of a neural network
is influenced more by recently and frequently observed data, the neural
network forgets what it has learned in prior tasks without continuing access
to them in the following tasks. A rigorous approach that maintains the
knowledge of entire tasks may solve the problem while sacrificing
computational cost, however, it is impractical with an undefined number of
tasks in real applications of continual learning. Continual learning model has
to adapt to a new task without access to some or entire classes of past tasks
while it maintains acquired knowledge from the past tasks [19]. In addition,
the continual learning model has to be evaluated with arbitrary order of tasks
since the order of tasks is not able to be fixed nor predicted in real
applications. The continual learning model is required to function
consistently regardless of the order of tasks.

Figure 1: Three different (Importance-Measurement) ways (a) WI-WM: Weight
Importance ($\Omega_{1}\sim\Omega_{l}$) by respective Weight Measurement (b)
NI-WM: Neuron Importance ($\Omega_{k}$) by Weight Measurements. The maximum
value of weight importance out of ($\Omega_{1}\sim\Omega_{l}$) is assigned to
neuron importance ($\Omega_{k}$). And then all weights connected to the neuron
get the same importance of the neuron. (c) NI-NM: Neuron Importance
($\Omega_{k}$) by Neuron Measurement, where $l$ and $\Omega$ indicate weight
index and importance of either weight or neuron respectively. The proposed
method belongs to (c) NI-NM.

(a) Task order: 1 $\rightarrow$ 2 $\rightarrow$ 3 $\rightarrow$ 4
$\rightarrow$ 5

(b) Task order: 2 $\rightarrow$ 4 $\rightarrow$ 3 $\rightarrow$ 5
$\rightarrow$ 1

Figure 2: Classification accuracy of continual learning on Split Cifar10. SI
[22], MAS [2] and UCL [1] show critical changes in their performance as the
order of tasks changes.

| Absolute Task Order
---|---
Method | T1 | T2 | T3 | T4 | T5
SI [22] | 11.9 | 5.9 | 2.4 | 1.5 | 9.48
MAS [2] | 11.4 | 5.2 | 7.0 | 3.3 | 4.5
UCL [1] | 10.2 | 3.1 | 2.0 | 5.0 | 4.7
OURS | 11.4 | 1.5 | 0.5 | 1.7 | 1.6

Table 1: Performance disparity(%) between Figure 2(a) and 2(b) on Split CIFAR
10. ”Absolute task order” represents the sequence of tasks that a model
learns. (Additional explanation is discussed in Section 3.)

There are three major categories in prior continual learning approaches; 1)
architecture modification of neural networks [21, 18, 17], 2) rehearsal using
sampled data from previous tasks [16, 3], and 3) regularization freezing
significant weights of a model calculating the importance of weights or
neurons [14, 10, 22, 15, 2, 4, 23, 1, 7, 8]. Most recent methods have tackled
the problem with fundamental regularization approaches that utilize the
weights of given networks to the fullest. The basic idea of regularization
approaches is to constrain essential weights of prior tasks not to change. In
general, they alleviate catastrophic interference with a new task by imposing
a penalty on the difference of weights between the prior tasks and the new
task. The extent of the penalty is controlled by the significance of weights
or neurons in solving a certain task using respective measurements. As Figure
1 illustrates, weight importance can be decided by three different
(Importance-Measurement) ways.

WI-WM (Weight Importance by Weight Measurement) [10, 22, 2, 15, 23] calculates
weight importance based on the measurement of the corresponding weight as
described in Figure 1a. Elastic weight consolidation (EWC) [10] estimates
parameter importance using the diagonal of the Fisher information matrix
equivalent to the second derivative of the loss. Synaptic intelligence (SI)
[22] measures the importance of weights in an online manner by calculating
each parameter’s sensitivity to the loss change while it trains a network.
When a certain parameter changes slightly during training batches but its
contribution to the loss is high (i.e., rapid change of its gradient), the
parameter is considered to be crucial and restricted not to be updated in
future tasks. Unlike SI [22], Memory aware synapses (MAS) [2] assesses the
contribution of each weight to the change of a learned function. It considers
the gradient of outputs of a model with a mean square error loss. Gradient
itself represents a change of outputs concerning the weights. Variational
Continual Learning (VCL) [15], a Bayesian neural network-based method, decides
weight importance through variational inference. Bayesian Gradient Descent
(BGD) [23] finds posterior parameters (e.g., mean and variance) assuming that
the posterior and the prior distributions are Gaussian.

To mitigate the interference across multiple tasks in continual learning,
weight importance-based approaches let each weight have its weight importance.
However, in the case of convolutional neural networks, since a convolutional
filter makes one feature map that can be regarded as one neuron, those weights
should have the same importance. Furthermore, those methods that consider the
amount of change of weights [10, 22, 2] are impossible to reinitialize weights
at each training of a new task, which possibly decreases the plasticity of the
network. (Additional explanation of weight re-initialization is discussed in
section 3.)

NI-WM (Neuron Importance by Weight Measurement) calculates neuron importance
based on the measurement of all weights. Weight importance is redefined as the
importance of its connected neuron [1]. Uncertainty-regularized Continual
Learning (UCL) [1] measures weight importance by its uncertainty indicating
the variance of weight distribution. It claims that the distribution of
essential weights for past tasks has low variance, and such stable weights
during training a task are regarded as important weights not to forget. As
illustrated in Figure 1b, it suggests neuron-based importance in neural
networks. The smallest variance value (maximum importance) among the weights
incoming to and outgoing from a corresponding neuron decides the importance of
the neuron, and then the importance of all those weights is updated as the
neuron importance.

Figure 3: Normalized Weight importance distribution of each convolution layer.
To show the proportion of the average value of weight importance among layers,
we normalize the values to sum 1. This is based on the first task of Split
CIFAR 10 (task order: 3-1-2-4-5).

NI-NM (Neuron Importance by Neuron Measurement) calculates neuron importance
based on the measurement of the corresponding neuron [8]. Weight importance is
defined as the importance of its connected neuron.

[8] exploits proximal gradient descents using a neuron importance. Its neuron
importance depends on the average activation value. Activation value itself is
a measurement of neuron importance, and weights connected to the neuron get
identical weight importance.

One critical observation in prior experimental evaluations of existing
continual learning methods is that the accuracy of each task significantly
changes when the order of tasks is changed. As discussed in [20], proposing a
continual learning method robust to the order of tasks is another critical
aspect. Therefore, performance evaluation with fixed task order does not
coincide with the fundamental aim of continual learning where no dedicated
order of tasks is given in reality. Figure 2 shows sample test results of
state-of-the-art continual learning methods compared to our proposed method.
As summarized in Table 1, classification accuracy values of prior methods
fluctuate as the order of tasks changes(from Figure 2(a) to Figure 2(b)).

In this work, we propose a regularization approach for continual learning
assigning neuron importance by the measurement of average neuron activation.
As Figure 3 describes, we balance neuron importance distribution among layers
based on the average neuron activation divided by standard deviation, which is
critical to performance consistency along the changes of task order. We assign
calculated neuron importance to all weights of incoming edges connected to the
neuron.

A Neuron with high activation to the majority of instances is defined as an
essential neuron. We freeze essential neurons by freezing the weights of all
connected incoming edges (essential weights) during the learning of a new task
so that our model remembers past tasks. We propose to evaluate the robustness
to the order of tasks in a comprehensive manner in which we evaluate the
average and standard deviation of classification accuracy with multiple sets
of randomly shuffled tasks.

Our approach remembers past tasks robustly compared to recent regularization
methods [22, 15, 2, 1]. To measure performance fluctuation along the change of
task order, we evaluate our method with numerous shuffled orders. We
quantitatively evaluate our classification performance based on a measure of
interference from past tasks on MNIST [13, 6], CIFAR10, CIFAR100 [12] and Tiny
ImageNet [5] data sets.

Key contributions of our work include 1) a simple but intuitive and effective
continual learning method introducing activation based neuron importance, 2) a
comprehensive experimental evaluation framework on existing benchmark data
sets to evaluate not just the final accuracy of continual learning also the
robustness of the accuracy along the changes of the order of tasks. Based on
the evaluation framework, existing state-of-the-art methods and the proposed
method are evaluated.

##  2 Proposed Method

###  2.1 Neuron Importance by Average Neuron Activation

The proposed method extracts neuron importance based on the average activation
value of all instances. And then the neuron importance is assigned to all
weights of incoming edges connected to the neuron. In convolutional neural
networks, activation value of a neuron corresponds to the average value of one
feature map (i.e., global average pooling value). The average activation value
of neuron corresponds to the average of global average pooling value. The
average activation values at each layer are independently calculated but are
considered together. In other words, the individual average activation values
represent the importance of each neuron of a whole model. However, encoded
features at each layer describe different aspects of an input image and, as a
result, the average activation values at each layer should not be evaluated
together. Therefore, the average activation value is not able to fully
represent the characteristics of the essential neuron. Besides, in convolution
neural networks, the absolute magnitude of average activation value (i.e., the
average of global average pooling value) varies along the location of layer:
in high-level feature maps, the portion of activated area decreases. Due to
the difference in absolute average activation values across the layers,
weights of earlier layers tend to be considered more essential as Figure 4
shows. If the average activation value is used as neuron importance, networks
will prefer to keep the weights of earlier layers.

Figure 4: Normalized Weight importance distribution of each convolution layer.
To show the proportion of the average value of weight importance among layers,
we normalize the values to sum 1. Our method relaxes the tendency to
excessively consolidate weights of earlier layers. This is based on the first
task of Split CIFAR 10 (task order: 3-1-2-4-5).

Instead, we propose to use layer-wise average activation divided by the
respective standard deviation for neuron importance measurement. Compared to
the average activation-based neuron importance [8], ours prevents earlier
layers from getting excessive importance compared to other layers, which, in
turn, prevents a network from vulnerable to changing the order of tasks in
terms of forgetting past tasks. Figure 4 shows normalized average weight
importance of each layer(total 6 layers). Prior average activation based
regularization term assigns around 57% of total importance to layer 1(57%,
12%, 10%, 6%, 8%, 8%, respectively for the 6 layers). On the other hand, our
proposed regularization loss term assigns 26% of total importance to layer 1.
Furthermore, our method avoids assigning excessive importance to certain
layer(26%, 16%, 16%, 15%, 15%, 12%).

Then, why this improves the continual learning performance regardless of task
order? In prior works, more weights of lower layers tend to be frozen in
earlier tasks that eliminate the chance of upcoming tasks to build new low-
level feature sets. Only a new task that is fortunately able to rebuild
higher-layer features based on the frozen lower layer weights from previous
tasks could survive. On the other hand, ours keeps the balance of frozen
weights in all layers securing more freedom of feature descriptions for new
tasks in both lower and higher layers. Indeed, lower layer features such as
edges are not class (task) dependent features. Therefore, excessively freezing
lower layer features is not preferable in continual learning. Even though
tasks change, a new task may find alternative low-level features that have
high similarity with them of past tasks, as discussed in [11]. In order to
encode such relation, we propose to use the average and standard deviation of
neuron activation values at each layer. Our loss function is described as
follows.

$\displaystyle
L_{t}=\tilde{L_{t}}+\alpha\sum\limits_{l}\Omega_{k}^{t}(w_{l}^{t-1}-w_{l}^{t})^{2},$
(1)

where $\tilde{L_{t}}$ is loss of current task (e.g., cross entropy loss), $t$
is task index, $l$ is weight index, and $\Omega_{k}^{t}$ indicates $k^{th}$
neuron importance. $\alpha$ is a strength parameter to control the amount of
weights consolidation. Neuron importance is defined as follows.

$\begin{gathered}\Omega^{t}_{k}=\frac{\frac{1}{N_{t}}\sum\limits_{{i=1}}^{N_{t}}{f_{k}(x^{(t)}_{i})}}{\sigma+\epsilon},\\\
\sigma=\sqrt{\frac{\sum\limits_{{i=1}}^{N_{t}}{\\{{f_{k}(x^{(t)}_{i})-\frac{1}{N_{t}}\sum\limits_{{i=1}}^{N_{t}}{f_{k}(x^{(t)}_{i})}}\\}^{2}}}{N_{t}}},\end{gathered}$
(2)

where $N_{t}$ is the number of instances, $x$ is input, $k$ is neuron index,
$f_{k}(\cdot)$ is activation value (global average value, in the case of
convolution neural network), and $i$ is instance index. We introduce
$\epsilon$ to prevent the numerator from being zero when the standard
deviation becomes zero. Proposed method considers the variation of average
activation value among instances and the differences of average activation
value among different layers. It encourages freezing more weights of later
layers than earlier layers which are more likely to describe given task-
specific features.

Our experiments(Table 2 in Section 3.2) show that prior methods tend to forget
past tasks in learning new tasks. In the prior methods, weights of later
layers are more likely to change than weights of earlier layers during
learning a new task.

In general, if the essential weights of later layers of previous tasks change,
the network forgets past tasks and hardly recovers previous task-specific
features. On the other hand, even though weights of earlier layers of previous
tasks change, there are other chances to recover general low-level features
which are shared with following new tasks. Since our method puts relatively
more constraints on the weights of task-specific features not to change than
the prior methods(Figure 3), our method forgets past tasks less showing stable
performance along the change in the order of tasks.

###  2.2 Weight Re-initialization for Better Plasticity

In continual learning, networks have to not only avoid catastrophic forgetting
but also learn new tasks. According to the extent of difference in optimal
classification feature space of different tasks, optimized feature space in
the previous task might be significantly changed with a new task. In the
learning of a new task, we can let the model start either from random weights
or from optimized weights with previous tasks. Even though the optimized
weights on previous tasks can be considered as a set of random weights for a
new task, we avoid a situation where the optimized weights for one task work
as a local optimal for another similar task that may hinder new training from
obtaining new optimal weights through weight re-initialization. The situation
can be explained with $\Omega_{k}(w_{k}^{t-1}-w_{k}^{t})^{2}$ term in the loss
function of our network. During the learning of a new task, the network is
informed of past tasks by $\Omega_{k}(w_{k}^{t-1}-w_{k}^{t})^{2}$ term which
lets the network maintain essential weights of the past tasks assigning high
$\Omega_{k}$ values. In other words, $\Omega_{k}(w_{k}^{t-1}-w_{k}^{t})^{2}$
delivers the knowledge of previous tasks. Whatever the magnitude of
$\Omega_{k}$ is, however, $\Omega_{k}(w_{k}^{t-1}-w_{k}^{t})^{2}$ term is
ignored if $w_{k}^{t-1}$ almost equals to $w_{k}{t}$ already in the initial
epoch of the training of a new task, which prevents the network from learning
a new task. This situation is alleviated by weight re-initialization that
allows the value of $\Omega_{k}(w_{k}^{t-1}-w_{k}^{t})^{2}$ to be high enough
regardless of the magnitude of $\Omega_{k}$ in the training of a new task. In
this case, still the knowledge of previous tasks will be delivered by
$\Omega_{k}$ and affect the training of a new task.

##  3 Experimental Evaluations

We perform experimental evaluations of our method compared to existing state-
of-the-art methods for continual learning on several benchmark data sets;
Split and permuted MNIST [13, 6], and incrementally learning classes of
CIFAR10, CIFAR100 [12] and Tiny ImageNet [5]. We set hyper-parameters of other
existing approaches based on the description in [1] which has tested existing
approaches with different hyper-parameters to find their best performance. We
train all different tasks with a batch size of 256 and Adam [9] using the same
learning rate (0.001). For the Split CIFAR tasks and Split Tiny ImageNet, as
aforementioned, we perform the evaluation multiple times shuffling the order
of tasks randomly to evaluate the robustness to task orders. We test with all
120, 200, and 50 random orders for Split CIFAR10, Split CIFAR10-100 and Split
Tiny ImageNet respectively. To minimize statistical fluctuations of accuracy,
each combination of task sequences is repeated three times.

Figure 5: Evaluation Metrics: DOI(Degree of Interference) and LA Accuracy of
task. T, L and n stands for task, learning step, and the number of tasks
respectively.

As described in Figure 5, we define several evaluation metrics. ”Absolute task
order” indicates the sequence of tasks that a model learns. For instance, task
1 stands for the first task that a model learns no matter which classes
comprise the task. ”Learning step-wise average accuracy(LA Accuracy)”
represents the accuracy of each learning step averaged through the whole tasks
involved. (i.e., $LA_{k}=Average(L_{k})$). ”Degree of interference(DOI)”
indicates the decreased extent of accuracy of each task after all learning
steps are conducted. It is calculated by $(T_{k},L_{k})-(T_{k},L_{n})$. When
we report the performance of randomly shuffled order experiment, we
respectively average LA accuracy and DOI of randomly shuffled ordered test.

###  3.1 MNIST

We first evaluate our algorithm on a Split MNIST benchmark. In this
experiment, two sequential classes compose each task (total 5 tasks). We use
multi-headed and multi-layer perceptrons with two hidden layers with 400 ReLU
activations. Each task has its output layer with two outputs and Softmax. We
train our network for 40 epochs with $\alpha$ = 0.0045. In Figure 6, we
compare the accuracy of each task for at every learning step (column-wise
comparison in Figure 5) and LA accuracy. MAS [2] outperforms all other
baselines reaching 99.81% while ours achieves 99.7%. However, the accuracy is
almost saturated due to the low complexity of the data.

Figure 6: Results on Split MNIST benchmark. Here, VCL indicates VCL(without
coreset)[15].

We also evaluate methods on permuted MNIST data set. Our model used in this
evaluation is MLP which consists of two hidden layers with 400 ReLUs each and
one output layer with Softmax. The network is trained for 20 epochs with
$\lambda$ = 0.005. Also, to normalize the range of activation value, ReLU is
applied to the output layer additionally when computing neuron importance
$\Omega_{k}$. Our algorithm (95.21%) outperforms MAS [2] (94.70%), EWC [10]
(82.45%) and VCL(without coreset) [15] (89.76%) and on the other hand, UCL [1]
(96.72%), SI [22] (96.39%) and BGD [23] (96.168%) show better results.
However, most results on this data set achieve almost saturated accuracy.

###  3.2 Split CIFAR10

We test our method on a Split CIFAR10 benchmark. In this experiment, two
sequential classes compose each task (total 5 tasks). Evaluation on Split
CIFAR10 data set is based on the multi-headed network with six convolution
layers and two fully connected layers where the output layer is different for
each task. We train our network for 100 epochs with $\alpha=0.7$. The order of
5 tasks that comprise CIFAR10 is randomly shuffled (total 120 random orders).

Figure 7: Average LA Accuracy and its std. of Split CIFAR10 benchmark.

| Average DOI of Absolute Task Order
---|---
Method | D1 | D2 | D3 | D4
SI [22] | 28.05(±11.7) | 20.00(±7.4) | 15.51(±8.2) | 8.68(±5.8)
MAS [2] | 33.59(±11.7) | 27.37(±11.3) | 19.15(±10.6) | 11.45(±6.7)
UCL [1] | 11.36(±5.8) | 8.56(±3.6) | 5.94(±3.0) | 3.55(±6.5)
BGD [23] | 39.06(±10.1) | 34.83(±8.5) | 29.19(±8.8) | 19.71(±2.1)
OURS | 1.44(±1.1) | 1.59(±1.2) | 1.18(±0.8) | 0.70(±0.7)

Table 2: Average DOI(Degree of interference) and its std.(%) on Split CIFAR10.
Note that proposed method forgets past tasks less regardless of the order of
tasks.

As Figure 7 describes, our method overall outperforms all other methods with
large margins. Also, the standard deviation graph shows that our algorithm is
more robust to the order of tasks.

As Table 2 shows, proposed method shows better stability in the order of tasks
and also has a low degree of forgetting. In our method, average degraded
degree of performance is lowest as $1.23\%$, whereas SI [22] is $18.06\%$, UCL
[1] is $7.35\%$, MAS [2] is $22.89\%$, and BGD [23] is $30.7\%$.

####  3.2.1 Ablation study

To verify the effect of weight re-initialization for the learning of new
tasks, we compare performance of ours and UCL [1] with those without weight
re-initialization.

| Task Order
---|---
Method | T5 | T4 | T3 | T2 | T1
UCL [1] | 0 | -0.425 | -0.9 | 4.93 | 6.38
OURS | -1 | 15.7 | 21.44 | 18.62 | 21.44

Table 3: Performance difference(%) = (accuracy with weight re-initialization)
- (accuracy without weight re-initialization). Note that the task order is
fixed.

As Table 3 indicates, accuracy increases in both methods when weight re-
initialization is applied. It suggests that weight re-initialization
encourages better plasticity. Note that several weight importance based
methods [10, 22, 2] cannot employ weight re-initialization since they consider
the amount of weight changes in the methods.

###  3.3 Split CIFAR10-100

Figure 8: Average LA Accuracy and its std. of CIFAR10-100 benchmark.

| Average DOI of Absolute Task Order
---|---
Method | D1 | D2 | D3 | D4 | D5 | D6 | D7 | D8 | D9 | D10
SI [22] | 5.85(±0.9) | 7.37(±2.3) | 6.58(±2.1) | 5.87(±2.0) | 5.57(±1.7) | 5.09(±1.7) | 4.45(±1.5) | 3.97(±1.3) | 3.28(±1.2) | 2.17(±1.0)
MAS [2] | 9.32(±1.4) | 9.18(±2.8) | 8.19(±2.2) | 7.39(±2.1) | 6.65(±1.9) | 6.17(±1.9) | 5.30(±1.6) | 4.64(±1.4) | 3.70(±1.2) | 2.50(±1.0)
UCL [1] | 3.74(±0.4) | 1.20(±1.2) | 1.31(±1.1) | 0.98(±1.0) | 0.78(±0.9) | 0.68(±0.8) | 0.64(±0.7) | 0.59(±0.7) | 0.37(±0.5) | 0.27(±0.4)
OURS | 2.03(±0.4) | 2.08(±0.9) | 1.26(±0.9) | 0.94(±0 8) | 0.84(±0.8) | 0.77(±0.7) | 0.61(±0.7) | 0.59(±0.6) | 0.47(±0.5) | 0.24(±0.4)

Table 4: Average DOI(Degree of interference) and its std.(%) on Split
CIFAR10-100. Note that proposed method forgets past tasks less regardless of
the order of tasks.

We evaluate our method on Split CIFAR10-100 benchmark where each task has 10
consecutive classes (total 11 tasks). We use the same multi-headed setup as in
the case of Split CIFAR10. We train our network for 100 epochs with
$\alpha=0.5$. We fix task 1 as CIFAR10 due to the difference in the size of
data set between CIFAR10 and CIFAR100. The order of remaining tasks that
consist of CIFAR100 is randomly shuffled (total 200 random orders).

Our method shows better stability showing the best accuracy values in old
tasks. On the other hand, previous methods seem to prefer to be better with
recent new tasks proving that our importance based continual learning is
working appropriately. Indeed, as Figure 8 and Table 4 represent, SI [22] and
MAS [2] seem that they learn new tasks very well forgetting what they have
learned before.

(a) LA accuracy of Split Cifar10

(b) LA accuracy of Split Cifar10-100

Figure 9: The performance on Split CIFAR10 and CIFAR10-100 with doubled
channel. Accuracy increases when we use a doubled channel network. Note that
the task order is fixed.

Since all incoming weights are tied to the neuron in our method, the higher
number of weights to be consolidated during training new tasks causes lower
accuracy of final task. In practice, the decrease of plasticity in our method
can be addressed by using a larger network (e.g., the larger number of
channels). We test the performance with a network of a doubled number of
channels (256 to 512). Figure 9 shows that our network with doubled number of
channels has improved accuracy keeping its stability and better plasticity.

Table 4 shows that our method obtains lowest average degraded degree of
performance $0.98\%$ compared to SI [22], MAS [2], UCL [1] achieving $5.02\%$,
$6.3\%$, $1.06\%$ respectively. Also, the proposed method shows the lowest
standard deviation of DOI, which indicates that our method is robust to the
interference from various combinations of tasks.

Figure 10: Average LA Accuracy and its std. of Tiny ImageNet data set

| Average DOI of Absolute Task Order
---|---
Method | D1 | D2 | D3 | D4 | D5 | D6 | D7 | D8 | D9
SI [22] | 27.35(±3.7) | 29.02(±4.4) | 27.09(±4.9) | 22.67(±4.2) | 20.24(±3.7) | 18.25(±4.4) | 14.86(±3.8) | 11.70(±3.2) | 7.78(±3.0)
MAS [2] | 20.49(±4.2) | 16.32(±3.6) | 14.40(±4.1) | 11.02(±3.6) | 9.14(±3.6) | 8.27(±2.7) | 6.20(±2.3) | 5.41(±2.3) | 3.86(±1.8)
OURS | 11.47(±3.6) | 6.75(±2.8) | 5.78(±2.3) | 3.98(±1.5) | 3.16(±1.2) | 2.76(±1.5) | 2.15(±1.1) | 1.70(±0.9) | 0.98(±0.7)

Table 5: Average DOI(Degree of interference) and its std.(%) on Split Tiny
ImageNet. Note that proposed method forgets past tasks less regardless of the
order of tasks.

###  3.4 Split Tiny ImageNet

We evaluate our method on Split Tiny ImageNet data set where each task has 20
consecutive classes (total 10 tasks). We use the same multi-headed setup as in
the case of Split CIFAR10 and Split CIFAR10-100. We train our network for 100
epochs with $\alpha=0.5$. The order of tasks is randomly shuffled (total 50
random orders). Only convolution neural networks based methods are tested for
a fair comparison.

In Figure 10, our method outperforms all other methods with large margins. The
standard deviation graph shows that our method algorithm shows the least
performance disparity under the change in the order of tasks. Table 5 presents
that our method acquires lowest average degraded degree of performance among
SI [22], MAS [2] and ours, achieving $19.08\%$, $10.5\%$, and $4.3\%$
respectively. Also, ours has the lowest standard deviation of DOI. This
implies that our method is robust to the interference from various
combinations of tasks.

##  4 Conclusion

We have proposed an activation importance-based continual learning method that
consolidates important neurons of past tasks. Comprehensive evaluation has
proved that the proposed method has implemented regularization-based continual
learning achieving the fundamental aim of continual learning tasks not only
balancing between stability and plasticity but also keeping robustness of the
performance to the changes in the order of tasks.

## References

  * [1] Hongjoon Ahn, Sungmin Cha, Donggyu Lee, and Taesup Moon.  Uncertainty-based continual learning with adaptive regularization.  In Advances in Neural Information Processing Systems, pages 4392–4402, 2019.
  * [2] Rahaf Aljundi, Francesca Babiloni, Mohamed Elhoseiny, Marcus Rohrbach, and Tinne Tuytelaars.  Memory aware synapses: Learning what (not) to forget.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 139–154, 2018.
  * [3] Rahaf Aljundi, Min Lin, Baptiste Goujaud, and Yoshua Bengio.  Gradient based sample selection for online continual learning.  In Advances in Neural Information Processing Systems, pages 11816–11825, 2019.
  * [4] Rahaf Aljundi, Marcus Rohrbach, and Tinne Tuytelaars.  Selfless sequential learning.  arXiv preprint arXiv:1806.05421, 2018.
  * [5] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database, 2009.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 248–255, 2020.
  * [6] Ian J Goodfellow, Mehdi Mirza, Da Xiao, Aaron Courville, and Yoshua Bengio.  An empirical investigation of catastrophic forgetting in gradient-based neural networks.  arXiv preprint arXiv:1312.6211, 2013.
  * [7] Khurram Javed and Martha White.  Meta-learning representations for continual learning.  In Advances in Neural Information Processing Systems, pages 1820–1830, 2019.
  * [8] Sangwon Jung, Hongjoon Ahn, Sungmin Cha, and Taesup Moon.  Adaptive group sparse regularization for continual learning.  arXiv preprint arXiv:2003.13726, 2020.
  * [9] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [10] James Kirkpatrick, Razvan Pascanu, Neil Rabinowitz, Joel Veness, Guillaume Desjardins, Andrei A Rusu, Kieran Milan, John Quan, Tiago Ramalho, Agnieszka Grabska-Barwinska, et al.  Overcoming catastrophic forgetting in neural networks.  Proceedings of the national academy of sciences, 114(13):3521–3526, 2017.
  * [11] Simon Kornblith, Mohammad Norouzi, Honglak Lee, and Geoffrey Hinton.  Similarity of neural network representations revisited.  In International Conference on Machine Learning, pages 3519–3529. PMLR, 2019.
  * [12] Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * [13] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [14] Zhizhong Li and Derek Hoiem.  Learning without forgetting.  IEEE transactions on pattern analysis and machine intelligence, 40(12):2935–2947, 2017.
  * [15] Cuong V Nguyen, Yingzhen Li, Thang D Bui, and Richard E Turner.  Variational continual learning.  arXiv preprint arXiv:1710.10628, 2017.
  * [16] Matthew Riemer, Ignacio Cases, Robert Ajemian, Miao Liu, Irina Rish, Yuhai Tu, and Gerald Tesauro.  Learning to learn without forgetting by maximizing transfer and minimizing interference.  arXiv preprint arXiv:1810.11910, 2018.
  * [17] Andrei A Rusu, Neil C Rabinowitz, Guillaume Desjardins, Hubert Soyer, James Kirkpatrick, Koray Kavukcuoglu, Razvan Pascanu, and Raia Hadsell.  Progressive neural networks.  arXiv preprint arXiv:1606.04671, 2016.
  * [18] Ali Sharif Razavian, Hossein Azizpour, Josephine Sullivan, and Stefan Carlsson.  Cnn features off-the-shelf: an astounding baseline for recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition workshops, pages 806–813, 2014.
  * [19] Sebastian Thrun.  Is learning the n-th thing any easier than learning the first?  In Advances in neural information processing systems, pages 640–646, 1996.
  * [20] Jaehong Yoon, Saehoon Kim, Eunho Yang, and Sung Ju Hwang.  Scalable and order-robust continual learning with additive parameter decomposition.  arXiv preprint arXiv:1902.09432, 2019.
  * [21] Jaehong Yoon, Eunho Yang, Jeongtae Lee, and Sung Ju Hwang.  Lifelong learning with dynamically expandable networks.  arXiv preprint arXiv:1708.01547, 2017.
  * [22] Friedemann Zenke, Ben Poole, and Surya Ganguli.  Continual learning through synaptic intelligence.  Proceedings of machine learning research, 70:3987, 2017.
  * [23] Chen Zeno, Itay Golan, Elad Hoffer, and Daniel Soudry.  Task agnostic continual learning using online variational bayes.  arXiv preprint arXiv:1803.10123, 2018.